Index: palm/trunk/INSTALL/test_oceanml_p3d
===================================================================
--- palm/trunk/INSTALL/test_oceanml_p3d	(revision 3291)
+++ palm/trunk/INSTALL/test_oceanml_p3d	(revision 3294)
@@ -19,5 +19,4 @@
     initializing_actions       = 'set_constant_profiles', ! initial conditions
 
-    ocean                      = .TRUE., ! enable ocean modeling
     latitude                   = 55.6,
 
@@ -37,5 +36,4 @@
     bc_s_b                     = 'initial_gradient',
     bc_s_t                     = 'neumann',
-    bc_sa_t                    = 'neumann',
 
     use_top_fluxes             = .T.,
@@ -47,5 +45,4 @@
 
     top_heatflux               = 0.0,
-    top_salinityflux           = 0.0,
 !
 !-- numerics
@@ -75,6 +72,4 @@
 !-- data output
 !-------------------------------------------------------------------------------
-    netcdf_data_format         = 5, ! use NetCDF3
-
     dt_run_control             = 0.0,    ! output interval for run control
 
@@ -83,3 +78,15 @@
 / ! end of runtime parameter namelist
  
+!-------------------------------------------------------------------------------
+!-- OCEAN PARAMETER NAMELIST
+!   Documentation: https://palm.muk.uni-hannover.de/trac/wiki/doc/app/oceanpar
+!-------------------------------------------------------------------------------
+&ocean_parameters
+!
+!-- top boundary condition for salinity
+    bc_sa_t                    = 'neumann',
 
+!   zero salinityflux at ocean surface
+    top_salinityflux           = 0.0,
+
+/ ! end of ocean parameters
Index: palm/trunk/INSTALL/test_oceanml_rc
===================================================================
--- palm/trunk/INSTALL/test_oceanml_rc	(revision 3291)
+++ palm/trunk/INSTALL/test_oceanml_rc	(revision 3294)
@@ -1,10 +1,10 @@
 
  ******************************    --------------------------------------------
- * PALM 4.0  Rev: 2507M       *    ocean - run without 1D - prerun
+ * PALM 5.0  Rev: 3291M       *    ocean - run without 1D - prerun
  ******************************    --------------------------------------------
 
- Date:                 02-10-17    Run:       test_oceanml        
- Time:                 10:50:00    Run-No.:   00
- Run on host:             lcmuk
+ Date:                 01-10-18    Run:       test_oceanml                      
+ Time:                 04:15:16    Run-No.:   00
+ Run on host:              imuk
  Number of PEs:               4    Processor grid (x,y): (   2,   2) calculated
  ------------------------------------------------------------------------------
@@ -13,4 +13,5 @@
  -----------------
 
+ --> Use the Moeng_Wyngaard turbulence closure (LES mode).
  --> Use the boussinesq approximation for the model equations.
  --> Solve perturbation pressure via FFT using temperton-algorithm routines
@@ -34,6 +35,10 @@
  ----------------------------------
 
- Grid length:      dx =      1.250 m    dy =      1.250 m    dz =      1.250 m
+ Grid length:      dx =       1.250 m    dy =       1.250 m
+                dz(1) =       1.250 m
+
  Domain size:       x =     80.000 m     y =     80.000 m  z(u) =    -80.000 m
+
+
 
  Number of gridpoints (x,y,z):  (0:  63, 0:  63, 0:  65)
@@ -114,5 +119,5 @@
 
     1D-Profiles:
-       Output every             ******** s
+       Output every             9999999.90 s
 
 
@@ -141,8 +146,12 @@
  -------------------
 
-    Geograph. latitude  :   phi    = 55.6 degr
+    Geograph. latitude  :   latitude  = 55.6 degr
+    Geograph. longitude :   longitude =  0.0 degr
     Angular velocity    :   omega  = 0.729E-04 rad/s
     Coriolis parameter  :   f      =  0.000120 1/s
                             f*     =  0.000082 1/s
+
+    Day of the year at model start :   day_init      =     172
+    UTC time at model start        :   time_utc_init = 43200.0 s
 
     Gravity             :   g      =  9.8 m/s**2
Index: palm/trunk/SCRIPTS/palmrun
===================================================================
--- palm/trunk/SCRIPTS/palmrun	(revision 3291)
+++ palm/trunk/SCRIPTS/palmrun	(revision 3294)
@@ -27,4 +27,7 @@
 # -----------------
 # $Id$
+# more detailed error message for outdated fileconnection file
+# 
+# 3236 2018-09-10 09:07:06Z raasch
 # warning issued in case of file connection file version mismatch,
 # host_configuration renamed configuration_identifier,
@@ -409,5 +412,6 @@
           printf "\n               working directory, but its revision does not match"
           printf "\n               the revision of the default (trunk) version."
-          printf "\n               You may need to update your connection file!"
+          printf "\n               You may need to update your connection file"
+          printf "\n               \"${working_directory}/.palm.iofiles\" !"
        fi
     fi
Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 3291)
+++ palm/trunk/SOURCE/Makefile	(revision 3294)
@@ -25,4 +25,8 @@
 # -----------------
 # $Id$
+# changes related to modularization of the ocean mode,
+# bugfix: dependency to advec_ws was missed in chemistry_model_mod
+# 
+# 3274 2018-09-24 15:42:55Z knoop
 # Added palm dependency of multi_agent_system_mod, because of mas_last_actions
 # call at the end of palm run
@@ -494,5 +498,4 @@
 	disturb_field.f90 \
 	disturb_heatflux.f90 \
-	eqn_state_seawater.f90 \
 	exchange_horiz_2d.f90 \
 	exchange_horiz.f90 \
@@ -509,9 +512,9 @@
 	init_grid.f90 \
 	init_masks.f90 \
-	init_ocean.f90 \
 	init_pegrid.f90 \
 	init_pt_anomaly.f90 \
 	init_rankine.f90 \
 	init_slope.f90 \
+	init_vertical_profiles.f90 \
 	interaction_droplets_ptq.f90 \
 	land_surface_model_mod.f90 \
@@ -548,4 +551,5 @@
 	netcdf_data_input_mod.f90 \
 	netcdf_interface_mod.f90 \
+	ocean_mod.f90 \
 	outflow_turbulence.f90 \
 	package_parin.f90 \
@@ -707,4 +711,5 @@
 	mod_kinds.o \
 	modules.o \
+	ocean_mod.o \
 	radiation_model_mod.o \
 	turbulence_closure_mod.o \
@@ -753,4 +758,5 @@
 	chemistry_model_mod.o \
 	gust_mod.o \
+	init_vertical_profiles.o \
 	land_surface_model_mod.o \
 	large_scale_forcing_nudging_mod.o \
@@ -760,4 +766,5 @@
 	netcdf_data_input_mod.o \
 	netcdf_interface_mod.o \
+	ocean_mod.o \
 	plant_canopy_model_mod.o \
 	pmc_interface_mod.o \
@@ -772,4 +779,5 @@
 	wind_turbine_model_mod.o
 chemistry_model_mod.o: \
+	advec_ws.o \
 	chem_gasphase_mod.o \
 	chem_modules.o \
@@ -867,4 +875,5 @@
 	modules.o \
 	netcdf_interface_mod.o \
+	ocean_mod.o \
 	radiation_model_mod.o \
 	surface_mod.o \
@@ -882,4 +891,5 @@
 	modules.o \
 	netcdf_interface_mod.o \
+	ocean_mod.o \
 	turbulence_closure_mod.o \
 	urban_surface_mod.o
@@ -914,8 +924,4 @@
 	modules.o \
 	random_generator_parallel_mod.o \
-	surface_mod.o
-eqn_state_seawater.o: \
-	mod_kinds.o \
-	modules.o \
 	surface_mod.o
 exchange_horiz.o: \
@@ -963,4 +969,5 @@
 	modules.o \
 	netcdf_interface_mod.o \
+	ocean_mod.o \
 	plant_canopy_model_mod.o \
 	pmc_handle_communicator_mod.o \
@@ -993,4 +1000,5 @@
 	netcdf_data_input_mod.o \
 	netcdf_interface_mod.o \
+	ocean_mod.o \
 	plant_canopy_model_mod.o \
 	pmc_interface_mod.o \
@@ -1030,9 +1038,4 @@
 	modules.o \
 	netcdf_interface_mod.o
-init_ocean.o: \
-	basic_constants_and_equations_mod.o \
-	eqn_state_seawater.o \
-	mod_kinds.o \
-	modules.o
 init_pegrid.o: \
 	mod_kinds.o \
@@ -1052,4 +1055,7 @@
 	mod_kinds.o \
 	modules.o
+init_vertical_profiles.o: \
+	modules.o \
+	mod_kinds.o
 interaction_droplets_ptq.o: \
 	basic_constants_and_equations_mod.o \
@@ -1238,4 +1244,5 @@
 	mod_kinds.o \
 	modules.o \
+	ocean_mod.o \
 	radiation_model_mod.o \
 	spectra_mod.o \
@@ -1243,4 +1250,15 @@
 	urban_surface_mod.o \
 	uv_exposure_model_mod.o
+ocean_mod.o: \
+ advec_s_pw.o \
+	advec_s_up.o \
+	advec_ws.o \
+	basic_constants_and_equations_mod.o \
+	diffusion_s.o \
+	init_vertical_profiles.o \
+	mod_kinds.o \
+	modules.o \
+	pmc_interface_mod.o \
+	user_actions.o
 outflow_turbulence.o: \
 	cpulog_mod.o \
@@ -1278,4 +1296,5 @@
 	multi_agent_system_mod.o \
 	netcdf_interface_mod.o \
+	ocean_mod.o \
 	plant_canopy_model_mod.o \
 	pmc_interface_mod.o \
@@ -1387,5 +1406,4 @@
 	diffusion_v.o \
 	diffusion_w.o \
-	eqn_state_seawater.o \
 	large_scale_forcing_nudging_mod.o \
 	mod_kinds.o \
@@ -1393,4 +1411,5 @@
 	modules.o \
 	modules.o \
+	ocean_mod.o \
 	plant_canopy_model_mod.o \
 	radiation_model_mod.o \
@@ -1508,4 +1527,5 @@
 	mod_kinds.o \
 	modules.o \
+ ocean_mod.o \
 	pmc_interface_mod.o \
 	turbulence_closure_mod.o \
@@ -1538,4 +1558,5 @@
 	modules.o \
 	multi_agent_system_mod.o \
+ ocean_mod.o \
 	pmc_interface_mod.o \
 	prognostic_equations.o \
@@ -1595,4 +1616,5 @@
 	mod_particle_attributes.o \
 	modules.o \
+ ocean_mod.o \
 	plant_canopy_model_mod.o \
 	pmc_interface_mod.o \
@@ -1785,4 +1807,5 @@
 	modules.o \
 	netcdf_interface_mod.o \
+	ocean_mod.o \
 	radiation_model_mod.o \
 	random_function_mod.o \
Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 3291)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! ocean renamed ocean_mode
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -293,5 +296,5 @@
 
        USE control_parameters,                                                 &
-           ONLY:  humidity, loop_optimization, passive_scalar, ocean,          &
+           ONLY:  humidity, loop_optimization, passive_scalar, ocean_mode,     &
                   rans_tke_e, ws_scheme_mom, ws_scheme_sca
 
@@ -369,5 +372,5 @@
           ENDIF
 
-          IF ( ocean )  THEN
+          IF ( ocean_mode )  THEN
              ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) )
              sums_wssas_ws_l = 0.0_wp
@@ -454,5 +457,5 @@
              ENDIF
 
-             IF ( ocean )  THEN
+             IF ( ocean_mode )  THEN
                 ALLOCATE( flux_s_sa(nzb+1:nzt,0:threads_per_task-1),           &
                           diss_s_sa(nzb+1:nzt,0:threads_per_task-1) )
@@ -1084,5 +1087,6 @@
     
        USE control_parameters,                                                 &
-           ONLY:  humidity, passive_scalar, ocean, ws_scheme_mom, ws_scheme_sca
+           ONLY:  humidity, passive_scalar, ocean_mode, ws_scheme_mom,         &
+                  ws_scheme_sca
 
        USE kinds
@@ -1123,5 +1127,5 @@
              sums_wsnrs_ws_l = 0.0_wp
           ENDIF
-          IF ( ocean )  sums_wssas_ws_l = 0.0_wp
+          IF ( ocean_mode )  sums_wssas_ws_l = 0.0_wp
 
        ENDIF
Index: palm/trunk/SOURCE/average_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/average_3d_data.f90	(revision 3291)
+++ palm/trunk/SOURCE/average_3d_data.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! changes concerning modularization of ocean option
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -149,5 +152,5 @@
     USE control_parameters,                                                    &
         ONLY:  air_chemistry, average_count_3d, doav, doav_n, land_surface,    &
-               urban_surface, varnamelength
+               ocean_mode, urban_surface, varnamelength
 
     USE cpulog,                                                                &
@@ -165,4 +168,7 @@
         ONLY:  lsm_3d_data_averaging
 
+    USE ocean_mod,                                                             &
+        ONLY:  ocean_3d_data_averaging
+
     USE radiation_model_mod,                                                   &
         ONLY:  radiation, radiation_3d_data_averaging
@@ -179,8 +185,8 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i  !< running index
-    INTEGER(iwp) ::  ii !< running index
-    INTEGER(iwp) ::  j  !< running index
-    INTEGER(iwp) ::  k  !< running index
+    INTEGER(iwp) ::  i   !< loop index
+    INTEGER(iwp) ::  ii  !< loop index
+    INTEGER(iwp) ::  j   !< loop index
+    INTEGER(iwp) ::  k   !< loop index
 
     CHARACTER (LEN=varnamelength) ::  trimvar  !< TRIM of output-variable string
@@ -392,15 +398,4 @@
              ENDIF
 
-          CASE ( 'rho_ocean' )
-             IF ( ALLOCATED( rho_ocean_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) / REAL( average_count_3d, KIND=wp )
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
           CASE ( 's' )
              IF ( ALLOCATED( s_av ) ) THEN
@@ -414,15 +409,4 @@
              ENDIF
 
-          CASE ( 'sa' )
-             IF ( ALLOCATED( sa_av ) ) THEN
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         sa_av(k,j,i) = sa_av(k,j,i) / REAL( average_count_3d, KIND=wp )
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
          CASE ( 'shf*' )
              IF ( ALLOCATED( shf_av ) ) THEN
@@ -549,40 +533,39 @@
              ENDIF
 
+          CASE ( 'usm_output' )
 !             
-!--       Block of urban surface model outputs   
-          CASE ( 'usm_output' )
+!--          Block of urban surface model outputs
              CALL usm_average_3d_data( 'average', doav(ii) )
 
           CASE DEFAULT
 !
-!--          Turbulence closure module
-             CALL tcm_3d_data_averaging( 'average', doav(ii) )
-!
-!--          microphysics module quantities
+!--          Averaging of data from other modules
+             IF ( air_chemistry )  THEN
+                CALL chem_3d_data_averaging( 'average', doav(ii) )
+             ENDIF
+
              IF ( bulk_cloud_model )  THEN
                 CALL bcm_3d_data_averaging( 'average', doav(ii) )
              ENDIF
-!
-!--          Land surface quantity
+
+             IF ( gust_module_enabled )  THEN
+                CALL gust_3d_data_averaging( 'average', doav(ii) )
+             ENDIF
+
              IF ( land_surface )  THEN
                 CALL lsm_3d_data_averaging( 'average', doav(ii) )
              ENDIF
-!
-!--          Radiation quantity
+
+             IF ( ocean_mode )  THEN
+                CALL ocean_3d_data_averaging( 'average', doav(ii) )
+             ENDIF
+
              IF ( radiation )  THEN
                 CALL radiation_3d_data_averaging( 'average', doav(ii) )
              ENDIF
-!
-!--          Gust module quantities
-             IF ( gust_module_enabled )  THEN
-                CALL gust_3d_data_averaging( 'average', doav(ii) )
-             ENDIF
-!
-!--          Chemistry quantity
-             IF ( air_chemistry )  THEN
-                CALL chem_3d_data_averaging( 'average', doav(ii) )
-             ENDIF
-!
-!--          User-defined quantity
+
+             CALL tcm_3d_data_averaging( 'average', doav(ii) )
+!
+!--          User-defined quantities
              CALL user_3d_data_averaging( 'average', doav(ii) )
 
Index: palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- palm/trunk/SOURCE/boundary_conds.f90	(revision 3291)
+++ palm/trunk/SOURCE/boundary_conds.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! changes concerning modularization of ocean option
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -203,4 +206,7 @@
         ONLY:  kappa
 
+    USE bulk_cloud_model_mod,                                                  &
+        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
+
     USE chemistry_model_mod,                                                   &
         ONLY:  chem_boundary_conds 
@@ -212,8 +218,7 @@
                child_domain, constant_diffusion, coupling_mode, &
                dt_3d, humidity, ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, &
-               ibc_s_t,ibc_sa_t, ibc_uv_b, ibc_uv_t,                           &
-               intermediate_timestep_count,                                    &
-               nesting_offline, nudging,                                       &
-               ocean, passive_scalar, rans_mode, rans_tke_e, tsc, use_cmax
+               ibc_s_t, ibc_uv_b, ibc_uv_t,                           &
+               intermediate_timestep_count, nesting_offline, nudging,          &
+               ocean_mode, passive_scalar, rans_mode, rans_tke_e, tsc, use_cmax
 
     USE grid_variables,                                                        &
@@ -221,11 +226,10 @@
 
     USE indices,                                                               &
-        ONLY:  nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,             &
-               nzb, nzt
+        ONLY:  nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nzt
 
     USE kinds
 
-    USE bulk_cloud_model_mod,                                                  &
-        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
+    USE ocean_mod,                                                             &
+        ONLY:  ibc_sa_t
 
     USE pegrid
@@ -529,5 +533,5 @@
 !
 !-- Boundary conditions for salinity
-    IF ( ocean )  THEN
+    IF ( ocean_mode )  THEN
 !
 !--    Bottom boundary: Neumann condition because salinity flux is always
Index: palm/trunk/SOURCE/buoyancy.f90
===================================================================
--- palm/trunk/SOURCE/buoyancy.f90	(revision 3291)
+++ palm/trunk/SOURCE/buoyancy.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! module use statements moved to global declaration section
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -53,5 +56,4 @@
 ! 1850 2016-04-08 13:29:27Z maronga
 ! Module renamed
-! 
 ! 
 ! 1682 2015-10-07 23:56:08Z knoop
@@ -121,4 +123,13 @@
         ONLY:  g
 
+    USE arrays_3d,                                                             &
+        ONLY:  pt, pt_slope_ref, ref_state, tend
+
+    USE control_parameters,                                                    &
+        ONLY:  atmos_ocean_sign, cos_alpha_surface, message_string, pt_surface,&
+               sin_alpha_surface, sloping_surface
+
+    USE kinds
+
     PRIVATE
     PUBLIC buoyancy
@@ -139,18 +150,7 @@
     SUBROUTINE buoyancy( var, wind_component )
 
-       USE arrays_3d,                                                          &
-           ONLY:  pt, pt_slope_ref, ref_state, tend
-
-       USE control_parameters,                                                 &
-           ONLY:  atmos_ocean_sign, cos_alpha_surface, message_string,         &
-                  pt_surface, sin_alpha_surface, sloping_surface
-
        USE indices,                                                            &
            ONLY:  nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nzb,       &
                   nzt, wall_flags_0
-
-       USE kinds
-
-       USE pegrid
 
 
@@ -245,18 +245,7 @@
     SUBROUTINE buoyancy_ij( i, j, var, wind_component )
 
-       USE arrays_3d,                                                          &
-           ONLY:  pt, pt_slope_ref, ref_state, tend
-
-       USE control_parameters,                                                 &
-           ONLY:  atmos_ocean_sign, cos_alpha_surface, message_string,         &
-                  pt_surface, sin_alpha_surface, sloping_surface
 
        USE indices,                                                            &
            ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
-
-       USE kinds
-
-       USE pegrid
-
 
        IMPLICIT NONE
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 3291)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 3294)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id: check_parameters.f90 2520 2017-10-05 13:50:26Z gronemeier &
+! changes concerning modularization of ocean option,
+! init_vertical_profiles moved to separate file to avoid circular dependency
+! 
+! 2520 2017-10-05 13:50:26Z gronemeier
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -669,4 +673,6 @@
     USE arrays_3d
 
+    USE basic_constants_and_equations_mod
+
     USE bulk_cloud_model_mod,                                                  &
         ONLY:  bulk_cloud_model, bcm_check_parameters, bcm_check_data_output,  &
@@ -676,13 +682,21 @@
         ONLY:  chem_boundary_conds, chem_check_data_output,                    &
                chem_check_data_output_pr, chem_species
+
     USE chem_modules
-    USE basic_constants_and_equations_mod
+
     USE control_parameters
+
     USE dvrp_variables
+
     USE grid_variables
+
     USE gust_mod,                                                              &
         ONLY: gust_check_data_output, gust_check_data_output_pr,               &
               gust_check_parameters, gust_module_enabled
+
+    USE kinds
+
     USE indices
+
     USE land_surface_model_mod,                                                &
         ONLY: lsm_check_data_output, lsm_check_data_output_pr,                 &
@@ -691,6 +705,4 @@
     USE lsf_nudging_mod,                                                       &
         ONLY:  lsf_nudging_check_parameters, lsf_nudging_check_data_output_pr
-
-    USE kinds
 
     USE model_1d_mod,                                                          &
@@ -706,6 +718,12 @@
                waterflux_output_unit, momentumflux_output_unit
 
+    USE ocean_mod,                                                             &
+        ONLY:  ocean_check_data_output, ocean_check_data_output_pr,            &
+               ocean_check_parameters
+
     USE particle_attributes
+
     USE pegrid
+
     USE plant_canopy_model_mod,                                                &
         ONLY:  pcm_check_data_output, pcm_check_parameters
@@ -715,23 +733,35 @@
 
     USE profil_parameter
+
     USE radiation_model_mod,                                                   &
         ONLY:  radiation, radiation_check_data_output,                         &
                radiation_check_data_output_pr, radiation_check_parameters
+
     USE spectra_mod,                                                           &
         ONLY:  calculate_spectra, spectra_check_parameters
+
     USE statistics
+
     USE subsidence_mod
+
     USE statistics
+
     USE synthetic_turbulence_generator_mod,                                    &
         ONLY:  stg_check_parameters
+
     USE transpose_indices
+
     USE turbulence_closure_mod,                                                &
         ONLY:  tcm_check_data_output, tcm_check_parameters
+
     USE urban_surface_mod,                                                     &
         ONLY:  usm_check_data_output, usm_check_parameters
+
     USE uv_exposure_model_mod,                                                 &
         ONLY:  uvem_check_data_output
+
     USE wind_turbine_model_mod,                                                &
         ONLY:  wtm_check_parameters
+
     USE vertical_nesting_mod,                                                  &
         ONLY:  vnested, vnest_check_parameters
@@ -1083,8 +1113,4 @@
 
 !
-!-- Check turbulence closure setup
-    CALL tcm_check_parameters
-
-!
 !-- Check approximation
     IF ( TRIM( approximation ) /= 'boussinesq'   .AND.                         &
@@ -1127,4 +1153,5 @@
        CALL message( 'check_parameters', 'PA0450', 1, 2, 0, 6, 0 )
     ENDIF
+!
 !-- Set flux input mode according to approximation if applicable
     IF ( TRIM( flux_input_mode ) == 'approximation-specific' )  THEN
@@ -1145,4 +1172,5 @@
        CALL message( 'check_parameters', 'PA0451', 1, 2, 0, 6, 0 )
     ENDIF
+!
 !-- Set flux output mode according to approximation if applicable
     IF ( TRIM( flux_output_mode ) == 'approximation-specific' )  THEN
@@ -1163,5 +1191,5 @@
 
 !
-!-- set the flux output units according to flux_output_mode
+!-- Set the flux output units according to flux_output_mode
     IF ( TRIM( flux_output_mode ) == 'kinematic' ) THEN
         heatflux_output_unit              = 'K m/s'
@@ -1174,23 +1202,9 @@
     ENDIF
 
+!
 !-- set time series output units for fluxes
     dots_unit(14:16) = heatflux_output_unit
     dots_unit(21)    = waterflux_output_unit
     dots_unit(19:20) = momentumflux_output_unit
-
-!
-!-- Check ocean setting
-    IF ( TRIM( coupling_mode ) == 'uncoupled'  .AND.                           &
-         TRIM( coupling_char ) == '_O' .AND.                                   &
-         .NOT. ocean )  THEN
-
-!
-!--    Check whether an (uncoupled) atmospheric run has been declared as an
-!--    ocean run (this setting is done via palmrun-option -y)
-       message_string = 'ocean = .F. does not allow coupling_char = "' //      &
-                        TRIM( coupling_char ) // '" set by palmrun-option "-y"'
-       CALL message( 'check_parameters', 'PA0317', 1, 2, 0, 6, 0 )
-
-    ENDIF
 
 !
@@ -1430,40 +1444,21 @@
 
 !
-!-- Check for synthetic turbulence generator settings
+!-- Check the module settings
+    IF ( ocean_mode )  CALL ocean_check_parameters
+    CALL tcm_check_parameters
     CALL stg_check_parameters
-!
-!-- When plant canopy model is used, peform addtional checks
-    IF ( plant_canopy )  CALL pcm_check_parameters
-
-!
-!-- Additional checks for spectra
-    IF ( calculate_spectra )  CALL spectra_check_parameters
-!
-!-- When land surface model is used, perform additional checks
-    IF ( land_surface )  CALL lsm_check_parameters
-!
-!-- When microphysics module is used, perform additional checks
-    IF ( bulk_cloud_model )  CALL bcm_check_parameters
-
-!
-!-- When urban surface model is used, perform additional checks
-    IF ( urban_surface )  CALL usm_check_parameters
-
-!
-!-- If wind turbine model is used, perform additional checks
-    IF ( wind_turbine )  CALL wtm_check_parameters
-!
-!-- When radiation model is used, peform addtional checks
-    IF ( radiation )  CALL radiation_check_parameters
-!
-!-- When gust module is used, perform additional checks
+    IF ( plant_canopy )         CALL pcm_check_parameters
+    IF ( calculate_spectra )    CALL spectra_check_parameters
+    IF ( land_surface )         CALL lsm_check_parameters
+    IF ( bulk_cloud_model )     CALL bcm_check_parameters
+    IF ( urban_surface )        CALL usm_check_parameters
+    IF ( wind_turbine )         CALL wtm_check_parameters
+    IF ( radiation )            CALL radiation_check_parameters
     IF ( gust_module_enabled )  CALL gust_check_parameters
-!
-!-- When large-scale forcing or nudging is used, peform addtional checks
     IF ( large_scale_forcing  .OR.  nudging )  CALL lsf_nudging_check_parameters
 
 !
-!-- In case of no model continuation run, check initialising parameters and
-!-- deduce further quantities
+!-- In case of no restart run, check initialising parameters and calculate
+!-- further quantities
     IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
 
@@ -1472,5 +1467,4 @@
        pt_init = pt_surface
        IF ( humidity       )  q_init  = q_surface
-       IF ( ocean          )  sa_init = sa_surface
        IF ( passive_scalar )  s_init  = s_surface
        IF ( air_chemistry )  THEN
@@ -1486,5 +1480,5 @@
        gradient = 0.0_wp
 
-       IF (  .NOT.  ocean )  THEN
+       IF ( .NOT. ocean_mode )  THEN
 
           ug_vertical_gradient_level_ind(1) = 0
@@ -1550,5 +1544,5 @@
        gradient = 0.0_wp
 
-       IF (  .NOT.  ocean )  THEN
+       IF ( .NOT. ocean_mode )  THEN
 
           vg_vertical_gradient_level_ind(1) = 0
@@ -1691,23 +1685,12 @@
 !--    Compute initial chemistry profile using the given chemical species gradients
        IF ( air_chemistry )  THEN                                                          
-         DO lsp = 1, nvar
-          CALL init_vertical_profiles( cs_vertical_gradient_level_ind(lsp,:),  &
-                                       cs_vertical_gradient_level(lsp,:),      &
-                                       cs_vertical_gradient(lsp,:),            &
-                                       chem_species(lsp)%conc_pr_init,         &
-                                       cs_surface(lsp), bc_cs_t_val(lsp) )
-         ENDDO
-       ENDIF
-!
-!
-!--    If required, compute initial salinity profile using the given salinity
-!--    gradients
-       IF ( ocean )  THEN
-          CALL init_vertical_profiles( sa_vertical_gradient_level_ind,          &
-                                       sa_vertical_gradient_level,              &
-                                       sa_vertical_gradient, sa_init,           &
-                                       sa_surface, dum )
-       ENDIF
-
+          DO  lsp = 1, nvar
+            CALL init_vertical_profiles( cs_vertical_gradient_level_ind(lsp,:),&
+                                         cs_vertical_gradient_level(lsp,:),    &
+                                         cs_vertical_gradient(lsp,:),          &
+                                         chem_species(lsp)%conc_pr_init,       &
+                                         cs_surface(lsp), bc_cs_t_val(lsp) )
+          ENDDO
+       ENDIF
 
     ENDIF
@@ -1777,15 +1760,4 @@
                         TRIM( reference_state ) // '"'
        CALL message( 'check_parameters', 'PA0056', 1, 2, 0, 6, 0 )
-    ENDIF
-
-!
-!-- Sign of buoyancy/stability terms
-    IF ( ocean )  atmos_ocean_sign = -1.0_wp
-
-!
-!-- Ocean version must use flux boundary conditions at the top
-    IF ( ocean .AND. .NOT. use_top_fluxes )  THEN
-       message_string = 'use_top_fluxes must be .TRUE. in ocean mode'
-       CALL message( 'check_parameters', 'PA0042', 1, 2, 0, 6, 0 )
     ENDIF
 
@@ -2067,39 +2039,4 @@
 
 !
-!-- 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
-          message_string = 'unknown boundary condition: bc_sa_t = "' //        &
-                           TRIM( bc_sa_t ) // '"'
-          CALL message( 'check_parameters', 'PA0068', 1, 2, 0, 6, 0 )
-       ENDIF
-
-       IF ( top_salinityflux == 9999999.9_wp )  constant_top_salinityflux = .FALSE.
-       IF ( ibc_sa_t == 1  .AND.  top_salinityflux == 9999999.9_wp )  THEN
-          message_string = 'boundary condition: bc_sa_t = "' //                &
-                           TRIM( bc_sa_t ) // '" requires to set ' //          &
-                           'top_salinityflux'
-          CALL message( 'check_parameters', 'PA0069', 1, 2, 0, 6, 0 )
-       ENDIF
-
-!
-!--    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_wp )  THEN
-          message_string = 'boundary condition: bc_sa_t = "' //                &
-                           TRIM( bc_sa_t ) // '" is not allowed with ' //      &
-                           'top_salinityflux /= 0.0'
-          CALL message( 'check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
-       ENDIF
-
-    ENDIF
-
-!
 !-- Set boundary conditions for total water content
     IF ( humidity )  THEN
@@ -2158,5 +2095,5 @@
                .AND.  top_scalarflux /= 0.0_wp )  THEN
           message_string = 'boundary condition: bc_s_t = "' //                 &
-                           TRIM( bc_sa_t ) // '" is not allowed with ' //      &
+                           TRIM( bc_s_t ) // '" is not allowed with ' //       &
                            'top_scalarflux /= 0.0'
           CALL message( 'check_parameters', 'PA0441', 1, 2, 0, 6, 0 )
@@ -2230,5 +2167,5 @@
 
     IF ( rayleigh_damping_height == -1.0_wp )  THEN
-       IF (  .NOT.  ocean )  THEN
+       IF (  .NOT.  ocean_mode )  THEN
           rayleigh_damping_height = 0.66666666666_wp * zu(nzt)
        ELSE
@@ -2236,5 +2173,5 @@
        ENDIF
     ELSE
-       IF (  .NOT.  ocean )  THEN
+       IF (  .NOT.  ocean_mode )  THEN
           IF ( rayleigh_damping_height < 0.0_wp  .OR.                          &
                rayleigh_damping_height > zu(nzt) )  THEN
@@ -2519,22 +2456,4 @@
              hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
 
-          CASE ( 'sa', '#sa' )
-             IF ( .NOT. ocean )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for ocean = .FALSE.'
-                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
-             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_wp    ! because zu(nzb) is negative
-                   data_output_pr(i)     = data_output_pr(i)(2:)
-                ENDIF
-             ENDIF
-
           CASE ( 'u*2' )
              dopr_index(i) = 30
@@ -2846,5 +2765,4 @@
              hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
 
-
           CASE ( 'u"pt"' )
              dopr_index(i) = 58
@@ -2876,58 +2794,4 @@
              dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
-
-          CASE ( 'rho_ocean' )
-             IF (  .NOT.  ocean ) THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for ocean = .FALSE.'
-                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
-             ELSE
-                dopr_index(i) = 64
-                dopr_unit(i)  = 'kg/m3'
-                hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
-                IF ( data_output_pr(i)(1:1) == '#' )  THEN
-                   dopr_initial_index(i) = 77
-                   hom(:,2,77,:)         = SPREAD( zu, 2, statistic_regions+1 )
-                   hom(nzb,2,77,:)       = 0.0_wp    ! because zu(nzb) is negative
-                   data_output_pr(i)     = data_output_pr(i)(2:)
-                ENDIF
-             ENDIF
-
-          CASE ( 'w"sa"' )
-             IF (  .NOT.  ocean ) THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for ocean = .FALSE.'
-                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
-             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
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for ocean = .FALSE.'
-                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
-             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
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for ocean = .FALSE.'
-                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
-             ELSE
-                dopr_index(i) = 67
-                dopr_unit(i)  = 'psu m/s'
-                hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 )
-             ENDIF
 
           CASE ( 'w*p*' )
@@ -2953,16 +2817,4 @@
              ENDIF
 
-          CASE ( 'prho' )
-             IF (  .NOT.  ocean ) THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for ocean = .FALSE.'
-                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
-             ELSE
-                dopr_index(i) = 71
-                dopr_unit(i)  = 'kg/m3'
-                hom(:,2,71,:) = SPREAD( zu, 2, statistic_regions+1 )
-             ENDIF
-
           CASE ( 'hyp' )
              dopr_index(i) = 72
@@ -3014,20 +2866,34 @@
              ENDIF
 
-
-
           CASE DEFAULT
+
+!
+!--          Check for other modules
+             IF ( unit == 'illegal'  .AND.  bulk_cloud_model  )  THEN
+                CALL bcm_check_data_output_pr( data_output_pr(i), i, unit,     &
+                                               dopr_unit(i) )
+             ENDIF
+
+             IF ( unit == 'illegal' )  THEN
+                CALL chem_check_data_output_pr( data_output_pr(i), i, unit,    &
+                                                dopr_unit(i) )
+             ENDIF
+
+             IF ( unit == 'illegal'  .AND.  gust_module_enabled  )  THEN
+                CALL gust_check_data_output_pr( data_output_pr(i), i, unit,    &
+                                                dopr_unit(i) )
+             ENDIF
+
+             IF ( unit == 'illegal' )  THEN
+                CALL lsf_nudging_check_data_output_pr( data_output_pr(i), i,   &
+                                                       unit, dopr_unit(i) )
+             ENDIF
 
              CALL lsm_check_data_output_pr( data_output_pr(i), i, unit,        &
                                             dopr_unit(i) )
-!
-!--          Block of microphysics module profile outputs
-             IF ( unit == 'illegal'  .AND.  bulk_cloud_model  )  THEN
-                   CALL bcm_check_data_output_pr(data_output_pr(i), i, unit,   &
+
+             IF ( unit == 'illegal' )  THEN
+                CALL ocean_check_data_output_pr( data_output_pr(i), i, unit,   &
                                                  dopr_unit(i) )
-             ENDIF
-
-             IF ( unit == 'illegal' )  THEN
-                CALL lsf_nudging_check_data_output_pr( data_output_pr(i), i,   &
-                                                    unit, dopr_unit(i) )
              ENDIF
 
@@ -3036,16 +2902,7 @@
                                                      unit, dopr_unit(i) )
              ENDIF
-!
-!--          Block of gust module profile outputs
-             IF ( unit == 'illegal'  .AND.  gust_module_enabled  )  THEN
-                   CALL gust_check_data_output_pr( data_output_pr(i), i, unit, &
-                                                   dopr_unit(i) )
-             ENDIF
-
-             IF ( unit == 'illegal' )  THEN                                        
-                CALL chem_check_data_output_pr( data_output_pr(i), i,          &
-                                                unit, dopr_unit(i) )
-             ENDIF 
-
+
+!
+!--          Finally, check for user defined quantities
              IF ( unit == 'illegal' )  THEN
                 unit = ''
@@ -3053,4 +2910,6 @@
              ENDIF
 
+!
+!--          No valid quantity found
              IF ( unit == 'illegal' )  THEN
                 IF ( data_output_pr_user(1) /= ' ' )  THEN
@@ -3184,12 +3043,4 @@
              unit = 'kg/kg'
 
-          CASE ( 'rho_ocean' )
-             IF (  .NOT.  ocean )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res ocean = .TRUE.'
-                CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
-             ENDIF
-             unit = 'kg/m3'
-
           CASE ( 's' )
              IF (  .NOT.  passive_scalar )  THEN
@@ -3199,12 +3050,4 @@
              ENDIF
              unit = 'kg/m3'
-
-          CASE ( 'sa' )
-             IF (  .NOT.  ocean )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res ocean = .TRUE.'
-                CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
-             ENDIF
-             unit = 'psu'
 
           CASE ( 'p', 'pt', 'u', 'v', 'w' )
@@ -3279,28 +3122,12 @@
           CASE DEFAULT
 
-             CALL tcm_check_data_output ( var, unit )
-
-!
-!--          Block of microphysics module outputs
+             CALL tcm_check_data_output( var, unit )
+
+!
+!--          Check for other modules
              IF ( unit == 'illegal'  .AND.  bulk_cloud_model  )  THEN
-                CALL bcm_check_data_output ( var, unit )
-             ENDIF
-
-             IF ( unit == 'illegal' )  THEN
-                CALL lsm_check_data_output ( var, unit, i, ilen, k )
-             ENDIF
-
-             IF ( unit == 'illegal' )  THEN
-                CALL radiation_check_data_output( var, unit, i, ilen, k )
-             ENDIF
-
-!
-!--          Block of gust module outputs
-             IF ( unit == 'illegal'  .AND.  gust_module_enabled  )  THEN
-                CALL gust_check_data_output ( var, unit )
-             ENDIF
-
-!
-!--          Block of chemistry model outputs
+                CALL bcm_check_data_output( var, unit )
+             ENDIF
+
              IF ( unit == 'illegal'  .AND.  air_chemistry                      &
                                      .AND.  var(1:3) == 'kc_' )  THEN
@@ -3308,22 +3135,37 @@
              ENDIF
 
-!
-!--          Block of urban surface model outputs
-             IF ( unit == 'illegal' .AND. urban_surface .AND. var(1:4) == 'usm_' ) THEN
+             IF ( unit == 'illegal' )  THEN
+                CALL lsm_check_data_output ( var, unit, i, ilen, k )
+             ENDIF
+
+             IF ( unit == 'illegal'  .AND.  gust_module_enabled  )  THEN
+                CALL gust_check_data_output ( var, unit )
+             ENDIF
+
+             IF ( unit == 'illegal'  .AND.  ocean_mode )  THEN
+                CALL ocean_check_data_output( var, unit )
+             ENDIF
+
+             IF ( unit == 'illegal'  .AND.  plant_canopy                       &
+                                     .AND.  var(1:4) == 'pcm_' )  THEN
+                CALL pcm_check_data_output( var, unit )
+             ENDIF
+
+             IF ( unit == 'illegal' )  THEN
+                CALL radiation_check_data_output( var, unit, i, ilen, k )
+             ENDIF
+
+             IF ( unit == 'illegal'  .AND.  urban_surface                      &
+                                     .AND.  var(1:4) == 'usm_' )  THEN
                  CALL usm_check_data_output( var, unit )
              ENDIF
 
-!
-!--          Block of plant canopy model outputs
-             IF ( unit == 'illegal' .AND. plant_canopy .AND. var(1:4) == 'pcm_' ) THEN
-                 CALL pcm_check_data_output( var, unit )
-             ENDIF
-
-!
-!--          Block of uv exposure model outputs
-             IF ( unit == 'illegal' .AND. uv_exposure .AND. var(1:5) == 'uvem_' )  THEN
+             IF ( unit == 'illegal'  .AND.  uv_exposure                        &
+                                     .AND.  var(1:5) == 'uvem_' )  THEN
                 CALL uvem_check_data_output( var, unit, i, ilen, k )
              ENDIF
 
+!
+!--          Finally, check for user-defined quantities
              IF ( unit == 'illegal' )  THEN
                 unit = ''
@@ -3628,5 +3470,5 @@
 !-- Determine upper and lower hight level indices for random perturbations
     IF ( disturbance_level_b == -9999999.9_wp )  THEN
-       IF ( ocean )  THEN
+       IF ( ocean_mode )  THEN
           disturbance_level_b     = zu((nzt*2)/3)
           disturbance_level_ind_b = ( nzt * 2 ) / 3
@@ -3653,5 +3495,5 @@
 
     IF ( disturbance_level_t == -9999999.9_wp )  THEN
-       IF ( ocean )  THEN
+       IF ( ocean_mode )  THEN
           disturbance_level_t     = zu(nzt-3)
           disturbance_level_ind_t = nzt - 3
@@ -4074,106 +3916,4 @@
 
 
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Inititalizes the vertical profiles of scalar quantities.
-!------------------------------------------------------------------------------!
-
-    SUBROUTINE init_vertical_profiles( vertical_gradient_level_ind,            &
-                                       vertical_gradient_level,                &
-                                       vertical_gradient,                      &
-                                       pr_init, surface_value, bc_t_val )
-
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i     !< counter
-       INTEGER(iwp), DIMENSION(1:10) ::  vertical_gradient_level_ind !< vertical grid indices for gradient levels
-
-       REAL(wp)     ::  bc_t_val      !< model top gradient
-       REAL(wp)     ::  gradient      !< vertica gradient of the respective quantity
-       REAL(wp)     ::  surface_value !< surface value of the respecitve quantity
-
-
-       REAL(wp), DIMENSION(0:nz+1) ::  pr_init                 !< initialisation profile
-       REAL(wp), DIMENSION(1:10)   ::  vertical_gradient       !< given vertical gradient
-       REAL(wp), DIMENSION(1:10)   ::  vertical_gradient_level !< given vertical gradient level
-
-       i = 1
-       gradient = 0.0_wp
-
-       IF (  .NOT.  ocean )  THEN
-
-          vertical_gradient_level_ind(1) = 0
-          DO  k = 1, nzt+1
-             IF ( i < 11 )  THEN
-                IF ( vertical_gradient_level(i) < zu(k)  .AND.            &
-                     vertical_gradient_level(i) >= 0.0_wp )  THEN
-                   gradient = vertical_gradient(i) / 100.0_wp
-                   vertical_gradient_level_ind(i) = k - 1
-                   i = i + 1
-                ENDIF
-             ENDIF
-             IF ( gradient /= 0.0_wp )  THEN
-                IF ( k /= 1 )  THEN
-                   pr_init(k) = pr_init(k-1) + dzu(k) * gradient
-                ELSE
-                   pr_init(k) = pr_init(k-1) + dzu(k) * gradient
-                ENDIF
-             ELSE
-                pr_init(k) = pr_init(k-1)
-             ENDIF
-   !
-   !--       Avoid negative values
-             IF ( pr_init(k) < 0.0_wp )  THEN
-                pr_init(k) = 0.0_wp
-             ENDIF
-          ENDDO
-
-       ELSE
-
-          vertical_gradient_level_ind(1) = nzt+1
-          DO  k = nzt, 0, -1
-             IF ( i < 11 )  THEN
-                IF ( vertical_gradient_level(i) > zu(k)  .AND.            &
-                     vertical_gradient_level(i) <= 0.0_wp )  THEN
-                   gradient = vertical_gradient(i) / 100.0_wp
-                   vertical_gradient_level_ind(i) = k + 1
-                   i = i + 1
-                ENDIF
-             ENDIF
-             IF ( gradient /= 0.0_wp )  THEN
-                IF ( k /= nzt )  THEN
-                   pr_init(k) = pr_init(k+1) - dzu(k+1) * gradient
-                ELSE
-                   pr_init(k)   = surface_value - 0.5_wp * dzu(k+1) * gradient
-                   pr_init(k+1) = surface_value + 0.5_wp * dzu(k+1) * gradient
-                ENDIF
-             ELSE
-                pr_init(k) = pr_init(k+1)
-             ENDIF
-   !
-   !--       Avoid negative humidities
-             IF ( pr_init(k) < 0.0_wp )  THEN
-                pr_init(k) = 0.0_wp
-             ENDIF
-          ENDDO
-
-       ENDIF
-
-!
-!--    In case of no given gradients, choose zero gradient conditions
-       IF ( vertical_gradient_level(1) == -999999.9_wp )  THEN
-          vertical_gradient_level(1) = 0.0_wp
-       ENDIF
-!
-!--    Store gradient at the top boundary for possible Neumann boundary condition
-       bc_t_val  = ( pr_init(nzt+1) - pr_init(nzt) ) / dzu(nzt+1)
-
-    END SUBROUTINE init_vertical_profiles
-
-
-
 !------------------------------------------------------------------------------!
 ! Description:
Index: palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_2d.f90	(revision 3291)
+++ palm/trunk/SOURCE/data_output_2d.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! changes concerning modularization of ocean option
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -238,7 +241,6 @@
     USE arrays_3d,                                                             &
         ONLY:  dzw, e, heatflux_output_conversion, nc, nr, p, pt,              &
-               precipitation_amount, prr, q, qc, ql, ql_c, ql_v, qr,           &
-               rho_ocean, s, sa, tend, u, v, vpt, w, zu, zw,                   &
-               waterflux_output_conversion, hyrho, d_exner
+               precipitation_amount, prr, q, qc, ql, ql_c, ql_v, qr, s, tend,  &
+               u, v, vpt, w, zu, zw, waterflux_output_conversion, hyrho, d_exner
 
     USE averaging
@@ -258,5 +260,5 @@
                ibc_uv_b, io_blocks, io_group, land_surface, message_string,    &
                ntdim_2d_xy, ntdim_2d_xz, ntdim_2d_yz,                          &
-               psolver, section, simulated_time,                               &
+               ocean_mode, psolver, section, simulated_time,                   &
                time_since_reference_point, uv_exposure
 
@@ -284,4 +286,7 @@
                id_var_time_xy, id_var_time_xz, id_var_time_yz, nc_stat,        &
                netcdf_data_format, netcdf_handle_error
+
+    USE ocean_mod,                                                             &
+        ONLY:  ocean_data_output_2d
 
     USE particle_attributes,                                                   &
@@ -940,15 +945,4 @@
                 level_z(nzb+1) = zu(nzb+1)
 
-             CASE ( 'rho_ocean_xy', 'rho_ocean_xz', 'rho_ocean_yz' )
-                IF ( av == 0 )  THEN
-                   to_be_resorted => rho_ocean
-                ELSE
-                   IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN
-                      ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                      rho_ocean_av = REAL( fill_value, KIND = wp )
-                   ENDIF
-                   to_be_resorted => rho_ocean_av
-                ENDIF
-
              CASE ( 's_xy', 's_xz', 's_yz' )
                 IF ( av == 0 )  THEN
@@ -960,15 +954,4 @@
                    ENDIF
                    to_be_resorted => s_av
-                ENDIF
-
-             CASE ( 'sa_xy', 'sa_xz', 'sa_yz' )
-                IF ( av == 0 )  THEN
-                   to_be_resorted => sa
-                ELSE
-                   IF ( .NOT. ALLOCATED( sa_av ) ) THEN
-                      ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                      sa_av = REAL( fill_value, KIND = wp )
-                   ENDIF
-                   to_be_resorted => sa_av
                 ENDIF
 
@@ -1306,12 +1289,5 @@
 
 !
-!--             Turbulence closure variables
-                IF ( .NOT. found )  THEN
-                   CALL tcm_data_output_2d( av, do2d(av,if), found, grid, mode,&
-                                             local_pf, nzb_do, nzt_do )
-                ENDIF
-
-!
-!--             Microphysics module quantities
+!--             Quantities of other modules
                 IF ( .NOT. found  .AND.  bulk_cloud_model )  THEN
                    CALL bcm_data_output_2d( av, do2d(av,if), found, grid, mode,&
@@ -1319,6 +1295,9 @@
                 ENDIF
 
-!
-!--             Land surface model quantity
+                IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
+                   CALL gust_data_output_2d( av, do2d(av,if), found, grid,     &
+                                             local_pf, two_d, nzb_do, nzt_do )
+                ENDIF
+
                 IF ( .NOT. found  .AND.  land_surface )  THEN
                    CALL lsm_data_output_2d( av, do2d(av,if), found, grid, mode,&
@@ -1326,6 +1305,9 @@
                 ENDIF
 
-!
-!--             Radiation quantity
+                IF ( .NOT. found  .AND.  ocean_mode )  THEN
+                   CALL ocean_data_output_2d( av, do2d(av,if), found, grid,    &
+                                              mode, local_pf, nzb_do, nzt_do )
+                ENDIF
+
                 IF ( .NOT. found  .AND.  radiation )  THEN
                    CALL radiation_data_output_2d( av, do2d(av,if), found, grid,&
@@ -1334,13 +1316,9 @@
                 ENDIF
 
-!
-!--             Gust module quantities
-                IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
-                   CALL gust_data_output_2d( av, do2d(av,if), found, grid,     &
-                                             local_pf, two_d, nzb_do, nzt_do )
-                ENDIF
-
-!
-!--             UV exposure model quantity
+                IF ( .NOT. found )  THEN
+                   CALL tcm_data_output_2d( av, do2d(av,if), found, grid, mode,&
+                                             local_pf, nzb_do, nzt_do )
+                ENDIF
+
                 IF ( .NOT. found  .AND.  uv_exposure )  THEN
                    CALL uvem_data_output_2d( av, do2d(av,if), found, grid,     &
@@ -1349,5 +1327,5 @@
 
 !
-!--             User defined quantity
+!--             User defined quantities
                 IF ( .NOT. found )  THEN
                    CALL user_data_output_2d( av, do2d(av,if), found, grid,     &
Index: palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_3d.f90	(revision 3291)
+++ palm/trunk/SOURCE/data_output_3d.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! changes concerning modularization of ocean option
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -208,6 +211,6 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  e, nc, nr, p, pt, prr, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, &
-               sa, tend, u, v, vpt, w, d_exner
+        ONLY:  d_exner, e, nc, nr, p, pt, prr, q, qc, ql, ql_c, ql_v, qr, s,   &
+               tend, u, v, vpt, w
 
     USE averaging
@@ -225,5 +228,5 @@
         ONLY:  air_chemistry, do3d, do3d_no, do3d_time_count,                  &
                io_blocks, io_group, land_surface, message_string,              &
-               ntdim_3d, nz_do3d,  plant_canopy,                               &
+               ntdim_3d, nz_do3d,  ocean_mode, plant_canopy,                   &
                psolver, simulated_time, time_since_reference_point,            &
                urban_surface, varnamelength
@@ -252,4 +255,7 @@
                netcdf_data_format, netcdf_handle_error
         
+    USE ocean_mod,                                                             &
+        ONLY:  ocean_data_output_3d
+
     USE particle_attributes,                                                   &
         ONLY:  grid_particles, number_of_particles, particles,                 &
@@ -614,15 +620,4 @@
              ENDIF
 
-          CASE ( 'rho_ocean' )
-             IF ( av == 0 )  THEN
-                to_be_resorted => rho_ocean
-             ELSE
-                IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN
-                   ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                   u_av = REAL( fill_value, KIND = wp )
-                ENDIF
-                to_be_resorted => rho_ocean_av
-             ENDIF
-
           CASE ( 's' )
              IF ( av == 0 )  THEN
@@ -634,15 +629,4 @@
                 ENDIF
                 to_be_resorted => s_av
-             ENDIF
-
-          CASE ( 'sa' )
-             IF ( av == 0 )  THEN
-                to_be_resorted => sa
-             ELSE
-                IF ( .NOT. ALLOCATED( sa_av ) ) THEN
-                   ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                   sa_av = REAL( fill_value, KIND = wp )
-                ENDIF
-                to_be_resorted => sa_av
              ENDIF
 
@@ -701,12 +685,6 @@
           CASE DEFAULT
 
-             IF ( .NOT. found )  THEN
-                CALL tcm_data_output_3d( av, do3d(av,if), found, local_pf,     &
-                                         nzb_do, nzt_do )
-                resorted = .TRUE.
-             ENDIF
-
-!
-!--          Microphysic module quantities
+!
+!--          Quantities of other modules
              IF ( .NOT. found  .AND.  bulk_cloud_model )  THEN
                 CALL bcm_data_output_3d( av, do3d(av,if), found, local_pf,     &
@@ -715,6 +693,16 @@
              ENDIF
 
-!
-!--          Land surface quantity
+             IF ( .NOT. found  .AND.  air_chemistry )  THEN
+                CALL chem_data_output_3d( av, do3d(av,if), found,              &
+                                          local_pf, fill_value, nzb_do, nzt_do )
+                resorted = .TRUE.
+             ENDIF
+
+             IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
+                CALL gust_data_output_3d( av, do3d(av,if), found, local_pf,    &
+                                          nzb_do, nzt_do )
+                resorted = .TRUE.
+             ENDIF
+
              IF ( .NOT. found  .AND.  land_surface )  THEN
 !
@@ -742,6 +730,16 @@
              ENDIF
 
-!
-!--          Radiation quantity
+             IF ( .NOT. found  .AND.  ocean_mode )  THEN
+                CALL ocean_data_output_3d( av, do3d(av,if), found, local_pf,   &
+                                           nzb_do, nzt_do )
+                resorted = .TRUE.
+             ENDIF
+
+             IF ( .NOT. found  .AND.  plant_canopy )  THEN
+                CALL pcm_data_output_3d( av, do3d(av,if), found, local_pf,     &
+                                         fill_value, nzb_do, nzt_do )
+                resorted = .TRUE.
+             ENDIF
+
              IF ( .NOT. found  .AND.  radiation )  THEN
                 CALL radiation_data_output_3d( av, do3d(av,if), found,         &
@@ -750,30 +748,12 @@
              ENDIF
 
-!
-!--          Gust module quantities
-             IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
-                CALL gust_data_output_3d( av, do3d(av,if), found, local_pf,    &
-                                          nzb_do, nzt_do )
-                resorted = .TRUE.
-             ENDIF
-
-!
-!--          Chemistry model output
-             IF ( .NOT. found  .AND.  air_chemistry )  THEN
-                CALL chem_data_output_3d( av, do3d(av,if), found,              &
-                                          local_pf, fill_value, nzb_do, nzt_do )
-                resorted = .TRUE.
-             ENDIF
-
-!
-!--          Plant canopy model output
-             IF ( .NOT. found  .AND.  plant_canopy )  THEN
-                CALL pcm_data_output_3d( av, do3d(av,if), found, local_pf,     &
-                                         fill_value, nzb_do, nzt_do )
-                resorted = .TRUE.
-             ENDIF
-
-!
-!--          User defined quantity
+             IF ( .NOT. found )  THEN
+                CALL tcm_data_output_3d( av, do3d(av,if), found, local_pf,     &
+                                         nzb_do, nzt_do )
+                resorted = .TRUE.
+             ENDIF
+
+!
+!--          User defined quantities
              IF ( .NOT. found )  THEN
                 CALL user_data_output_3d( av, do3d(av,if), found, local_pf,    &
Index: palm/trunk/SOURCE/eqn_state_seawater.f90
===================================================================
--- palm/trunk/SOURCE/eqn_state_seawater.f90	(revision 3291)
+++ 	(revision )
@@ -1,397 +1,0 @@
-!> @file eqn_state_seawater.f90
-!------------------------------------------------------------------------------!
-! This file is part of the PALM model system.
-!
-! PALM is free software: you can redistribute it and/or modify it under the
-! terms of the GNU General Public License as published by the Free Software
-! Foundation, either version 3 of the License, or (at your option) any later
-! version.
-!
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
-! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along with
-! PALM. If not, see <http://www.gnu.org/licenses/>.
-!
-! Copyright 1997-2018 Leibniz Universitaet Hannover
-!------------------------------------------------------------------------------!
-!
-! Current revisions:
-! -----------------
-! 
-! 
-! Former revisions:
-! -----------------
-! $Id$
-! unused variables removed
-! 
-! 2718 2018-01-02 08:49:38Z maronga
-! Corrected "Former revisions" section
-! 
-! 2696 2017-12-14 17:12:51Z kanani
-! Change in file header (GPL part)
-!
-! 2369 2017-08-22 15:20:37Z suehring
-! Bugfix, do not mask topography here, since density becomes zero, leading to
-! division by zero in production_e
-! 
-! 2233 2017-05-30 18:08:54Z suehring
-!
-! 2232 2017-05-30 17:47:52Z suehring
-! Adjustments to new topography and surface concept
-! 
-! 2031 2016-10-21 15:11:58Z knoop
-! renamed variable rho to rho_ocean
-! 
-! 2000 2016-08-20 18:09:15Z knoop
-! Forced header and separation lines into 80 columns
-! 
-! 1873 2016-04-18 14:50:06Z maronga
-! Module renamed (removed _mod)
-! 
-! 
-! 1850 2016-04-08 13:29:27Z maronga
-! Module renamed
-! 
-! 
-! 1682 2015-10-07 23:56:08Z knoop
-! Code annotations made doxygen readable 
-! 
-! 1353 2014-04-08 15:21:23Z heinze
-! REAL constants provided with KIND-attribute 
-!
-! 1320 2014-03-20 08:40:49Z raasch
-! ONLY-attribute added to USE-statements,
-! kind-parameters added to all INTEGER and REAL declaration statements, 
-! kinds are defined in new module kinds, 
-! revision history before 2012 removed,
-! comment fields (!:) to be used for variable explanations added to
-! all variable declaration statements 
-!
-! 1036 2012-10-22 13:43:42Z raasch
-! code put under GPL (PALM 3.9)
-!
-! 97 2007-06-21 08:23:15Z raasch
-! 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.
-!> eqn_state_seawater calculates the potential density referred at hyp(0).
-!> eqn_state_seawater_func calculates density.
-!------------------------------------------------------------------------------!
- MODULE eqn_state_seawater_mod
- 
-    
-    USE kinds
-
-    IMPLICIT NONE
-
-    PRIVATE
-    PUBLIC eqn_state_seawater, eqn_state_seawater_func
-
-    REAL(wp), 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(wp), DIMENSION(13), PARAMETER ::  den =                               &
-                          (/ 1.0D0,                  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
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Call for all grid points
-!------------------------------------------------------------------------------!
-    SUBROUTINE eqn_state_seawater
-
-       USE arrays_3d,                                                          &
-           ONLY:  hyp, prho, pt_p, rho_ocean, sa_p
-       USE indices,                                                            &
-           ONLY:  nxl, nxr, nyn, nys, nzb, nzt
-
-       USE surface_mod,                                                        &
-          ONLY :  bc_h
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i       !< running index x direction
-       INTEGER(iwp) ::  j       !< running index y direction
-       INTEGER(iwp) ::  k       !< running index z direction
-       INTEGER(iwp) ::  m       !< running index surface elements
-
-       REAL(wp) ::  pden   !<
-       REAL(wp) ::  pnom   !<
-       REAL(wp) ::  p1     !<
-       REAL(wp) ::  p2     !<
-       REAL(wp) ::  p3     !<
-       REAL(wp) ::  pt1    !<
-       REAL(wp) ::  pt2    !<
-       REAL(wp) ::  pt3    !<
-       REAL(wp) ::  pt4    !<
-       REAL(wp) ::  sa1    !<
-       REAL(wp) ::  sa15   !<
-       REAL(wp) ::  sa2    !<
-       
-                       
-
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt
-!
-!--             Pressure is needed in dbar
-                p1 = hyp(k) * 1E-4_wp
-                p2 = p1 * p1
-                p3 = p2 * p1
-
-!
-!--             Temperature needed in degree Celsius
-                pt1 = pt_p(k,j,i) - 273.15_wp
-                pt2 = pt1 * pt1
-                pt3 = pt1 * pt2
-                pt4 = pt2 * pt2
-
-                sa1  = sa_p(k,j,i)
-                sa15 = sa1 * SQRT( sa1 )
-                sa2  = sa1 * sa1
-
-                pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     +    &
-                       nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 +    &
-                       nom(7)*sa2
-
-                pden = 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
-!
-!--             Potential density (without pressure terms)
-                prho(k,j,i) = pnom / pden 
-
-                pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  +    &
-                       nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2
-
-                pden = pden +             den(11)*p1     + den(12)*p2*pt3 +    &
-                       den(13)*p3*pt1
-
-!
-!--             In-situ density
-                rho_ocean(k,j,i) = pnom / pden 
-
-             ENDDO
-!
-!--          Neumann conditions are assumed at top boundary
-             prho(nzt+1,j,i)      = prho(nzt,j,i)
-             rho_ocean(nzt+1,j,i) = rho_ocean(nzt,j,i)
-
-          ENDDO
-       ENDDO
-!
-!--    Neumann conditions at up/downward-facing surfaces
-       !$OMP PARALLEL DO PRIVATE( i, j, k )
-       DO  m = 1, bc_h(0)%ns
-          i = bc_h(0)%i(m)            
-          j = bc_h(0)%j(m)
-          k = bc_h(0)%k(m)
-          prho(k-1,j,i)      = prho(k,j,i)
-          rho_ocean(k-1,j,i) = rho_ocean(k,j,i)
-       ENDDO
-!
-!--    Downward facing surfaces
-       !$OMP PARALLEL DO PRIVATE( i, j, k )
-       DO  m = 1, bc_h(1)%ns
-          i = bc_h(1)%i(m)            
-          j = bc_h(1)%j(m)
-          k = bc_h(1)%k(m)
-          prho(k+1,j,i)      = prho(k,j,i)
-          rho_ocean(k+1,j,i) = rho_ocean(k,j,i)
-       ENDDO
-
-    END SUBROUTINE eqn_state_seawater
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Call for grid point i,j
-!------------------------------------------------------------------------------!
-    SUBROUTINE eqn_state_seawater_ij( i, j )
-
-       USE arrays_3d,                                                          &
-           ONLY:  hyp, prho, pt_p, rho_ocean, sa_p
-           
-       USE indices,                                                            &
-           ONLY:  nzb, nzt
-
-       USE surface_mod,                                                        &
-          ONLY :  bc_h
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i       !< running index x direction
-       INTEGER(iwp) ::  j       !< running index y direction
-       INTEGER(iwp) ::  k       !< running index z direction
-       INTEGER(iwp) ::  m       !< running index surface elements
-       INTEGER(iwp) ::  surf_e  !< End index of surface elements at (j,i)-gridpoint
-       INTEGER(iwp) ::  surf_s  !< Start index of surface elements at (j,i)-gridpoint
-
-       REAL(wp) ::  pden   !<
-       REAL(wp) ::  pnom   !<
-       REAL(wp) ::  p1     !<
-       REAL(wp) ::  p2     !<
-       REAL(wp) ::  p3     !<
-       REAL(wp) ::  pt1    !<
-       REAL(wp) ::  pt2    !<
-       REAL(wp) ::  pt3    !<
-       REAL(wp) ::  pt4    !<
-       REAL(wp) ::  sa1    !<
-       REAL(wp) ::  sa15   !<
-       REAL(wp) ::  sa2    !<
-
-       DO  k = nzb+1, nzt
-!
-!--       Pressure is needed in dbar
-          p1 = hyp(k) * 1E-4_wp
-          p2 = p1 * p1
-          p3 = p2 * p1
-
-!
-!--       Temperature needed in degree Celsius
-          pt1 = pt_p(k,j,i) - 273.15_wp
-          pt2 = pt1 * pt1
-          pt3 = pt1 * pt2
-          pt4 = pt2 * pt2
-
-          sa1  = sa_p(k,j,i)
-          sa15 = sa1 * SQRT( sa1 )
-          sa2  = sa1 * sa1
-
-          pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     +          &
-                 nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 +          &
-                 nom(7)*sa2
-
-          pden = 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
-!
-!--       Potential density (without pressure terms)
-          prho(k,j,i) = pnom / pden 
-
-          pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  +          &
-                 nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2
-          pden = pden +             den(11)*p1     + den(12)*p2*pt3 +          &
-                 den(13)*p3*pt1
-
-!
-!--       In-situ density
-          rho_ocean(k,j,i) = pnom / pden 
-
-
-       ENDDO
-!
-!--    Neumann conditions at up/downward-facing walls
-       surf_s = bc_h(0)%start_index(j,i)   
-       surf_e = bc_h(0)%end_index(j,i)   
-       DO  m = surf_s, surf_e
-          k                  = bc_h(0)%k(m)
-          prho(k-1,j,i)      = prho(k,j,i)
-          rho_ocean(k-1,j,i) = rho_ocean(k,j,i)
-       ENDDO
-!
-!--    Downward facing surfaces
-       surf_s = bc_h(1)%start_index(j,i)   
-       surf_e = bc_h(1)%end_index(j,i)   
-       DO  m = surf_s, surf_e
-          k                  = bc_h(1)%k(m)
-          prho(k+1,j,i)      = prho(k,j,i)
-          rho_ocean(k+1,j,i) = rho_ocean(k,j,i)
-       ENDDO
-!
-!--    Neumann condition are assumed at top boundary
-       prho(nzt+1,j,i)      = prho(nzt,j,i)
-       rho_ocean(nzt+1,j,i) = rho_ocean(nzt,j,i)
-
-    END SUBROUTINE eqn_state_seawater_ij
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Equation of state as a function
-!------------------------------------------------------------------------------!
-    REAL(wp) FUNCTION eqn_state_seawater_func( p, pt, sa )
-
-       IMPLICIT NONE
-
-       REAL(wp) ::  p      !<
-       REAL(wp) ::  p1     !<
-       REAL(wp) ::  p2     !<
-       REAL(wp) ::  p3     !<
-       REAL(wp) ::  pt     !<
-       REAL(wp) ::  pt1    !<
-       REAL(wp) ::  pt2    !<
-       REAL(wp) ::  pt3    !<
-       REAL(wp) ::  pt4    !<
-       REAL(wp) ::  sa     !<
-       REAL(wp) ::  sa15   !<
-       REAL(wp) ::  sa2    !<
-
-!
-!--    Pressure is needed in dbar
-       p1 = p  * 1E-4_wp
-       p2 = p1 * p1
-       p3 = p2 * p1
-
-!
-!--    Temperature needed in degree Celsius
-       pt1 = pt - 273.15_wp
-       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 3291)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! ocean renamed ocean_mode
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -296,8 +299,8 @@
                 dt_3d, humidity, initializing_actions, land_surface,           &
                 large_scale_forcing, large_scale_subsidence, max_pr_user,      &
-                message_string, neutral, ocean, passive_scalar, simulated_time,&
-                simulated_time_at_begin, use_subsidence_tendencies,            &
-                use_surface_fluxes, use_top_fluxes, ws_scheme_mom,             &
-                ws_scheme_sca
+                message_string, neutral, ocean_mode, passive_scalar,           &
+                simulated_time, simulated_time_at_begin,                       &
+                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
+                ws_scheme_mom, ws_scheme_sca
         
     USE cpulog,                                                                &
@@ -432,5 +435,5 @@
 !--       According to the Neumann bc for the horizontal velocity components,
 !--       the corresponding fluxes has to satisfiy the same bc.
-          IF ( ocean )  THEN
+          IF ( ocean_mode )  THEN
              sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
              sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
@@ -459,5 +462,5 @@
              sums_l(:,17,i)                        = sums_wspts_ws_l(:,i)      &
                                            * heatflux_output_conversion  ! w*pt*
-             IF ( ocean          ) sums_l(:,66,i)  = sums_wssas_ws_l(:,i) ! w*sa*
+             IF ( ocean_mode     ) sums_l(:,66,i)  = sums_wssas_ws_l(:,i) ! w*sa*
              IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)       &
                                            * waterflux_output_conversion  ! w*q*
@@ -490,5 +493,5 @@
 !
 !--    Horizontally averaged profile of salinity
-       IF ( ocean )  THEN
+       IF ( ocean_mode )  THEN
           !$OMP DO
           DO  i = nxl, nxr
@@ -563,5 +566,5 @@
              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
+             IF ( ocean_mode )  THEN
                 sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
              ENDIF
@@ -592,5 +595,5 @@
        CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL,  &
                            MPI_SUM, comm2d, ierr )
-       IF ( ocean )  THEN
+       IF ( ocean_mode )  THEN
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
           CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb,       &
@@ -623,5 +626,5 @@
        sums(:,2) = sums_l(:,2,0)
        sums(:,4) = sums_l(:,4,0)
-       IF ( ocean )  sums(:,23) = sums_l(:,23,0)
+       IF ( ocean_mode )  sums(:,23) = sums_l(:,23,0)
        IF ( humidity ) THEN
           sums(:,44) = sums_l(:,44,0)
@@ -648,5 +651,5 @@
 !
 !--    Salinity
-       IF ( ocean )  THEN
+       IF ( ocean_mode )  THEN
           sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
           hom(:,1,23,sr) = sums(:,23)             ! sa
@@ -930,5 +933,5 @@
 !
 !--             Salinity flux w"sa"
-                IF ( ocean )  THEN
+                IF ( ocean_mode )  THEN
                    sums_l(k,65,tn) = sums_l(k,65,tn)                           &
                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
@@ -1000,5 +1003,5 @@
                          sums_l(k+ki,61,tn) = sums_l(k+ki,61,tn) + &
                                     0.0_wp * rmask(j,i,sr)        ! v"pt"
-                         IF ( ocean )  THEN
+                         IF ( ocean_mode )  THEN
                             sums_l(k+ki,65,tn) = sums_l(k+ki,65,tn) + &
                                        surf_def_h(l)%sasws(m) * rmask(j,i,sr)  ! w"sa"
@@ -1053,5 +1056,5 @@
                    sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
                                     0.0_wp * rmask(j,i,sr)        ! v"pt"
-                   IF ( ocean )  THEN
+                   IF ( ocean_mode )  THEN
                       sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
                                        surf_lsm_h%sasws(m) * rmask(j,i,sr)  ! w"sa"
@@ -1104,5 +1107,5 @@
                    sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
                                     0.0_wp * rmask(j,i,sr)        ! v"pt"
-                   IF ( ocean )  THEN
+                   IF ( ocean_mode )  THEN
                       sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
                                        surf_usm_h%sasws(m) * rmask(j,i,sr)  ! w"sa"
@@ -1276,5 +1279,5 @@
                                     0.0_wp * rmask(j,i,sr)        ! v"pt" 
 
-                IF ( ocean )  THEN
+                IF ( ocean_mode )  THEN
                    sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
                                        surf_def_h(2)%sasws(m) * rmask(j,i,sr)  ! w"sa"
@@ -1336,5 +1339,5 @@
 !--             Salinity flux and density (density does not belong to here,
 !--             but so far there is no other suitable place to calculate)
-                IF ( ocean )  THEN
+                IF ( ocean_mode )  THEN
                    IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
                       pts = 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +          &
@@ -1562,5 +1565,5 @@
 !
 !--    Density at top follows Neumann condition
-       IF ( ocean )  THEN
+       IF ( ocean_mode )  THEN
           sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
           sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
@@ -2032,5 +2035,5 @@
           first = .TRUE.
 
-          IF ( ocean )  THEN
+          IF ( ocean_mode )  THEN
              DO  k = nzt, nzb+1, -1
                 IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
@@ -2080,5 +2083,5 @@
           dptdz_threshold = 0.2_wp / 100.0_wp
 
-          IF ( ocean )  THEN
+          IF ( ocean_mode )  THEN
              DO  k = nzt+1, nzb+5, -1
                 IF ( dptdz(k) > dptdz_threshold  .AND.                         &
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 3291)
+++ palm/trunk/SOURCE/header.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! changes concerning modularization of ocean option
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -425,4 +428,8 @@
     USE netcdf_interface,                                                      &
         ONLY:  netcdf_data_format, netcdf_data_format_string, netcdf_deflate
+
+    USE ocean_mod,                                                             &
+        ONLY:  ibc_sa_t, prho_reference, sa_surface, sa_vertical_gradient,     &
+               sa_vertical_gradient_level, sa_vertical_gradient_level_ind
 
     USE particle_attributes,                                                   &
@@ -578,5 +585,5 @@
     ENDIF
     IF ( nested_run )  run_classification = 'nested ' // run_classification(1:63)
-    IF ( ocean )  THEN
+    IF ( ocean_mode )  THEN
        run_classification = 'ocean - ' // run_classification(1:61)
     ELSE
@@ -733,5 +740,5 @@
     IF ( use_upstream_for_tke )  WRITE ( io, 143 )
     IF ( rayleigh_damping_factor /= 0.0_wp )  THEN
-       IF ( .NOT. ocean )  THEN
+       IF ( .NOT. ocean_mode )  THEN
           WRITE ( io, 123 )  'above', rayleigh_damping_height, &
                rayleigh_damping_factor
@@ -815,5 +822,5 @@
 !
 !-- Computational grid
-    IF ( .NOT. ocean )  THEN
+    IF ( .NOT. ocean_mode )  THEN
        WRITE ( io, 250 )  dx, dy
        
@@ -1176,5 +1183,5 @@
     ENDIF
 
-    IF ( ocean )  THEN
+    IF ( ocean_mode )  THEN
        r_lower = 'sa(0)    = sa(1)'
        IF ( ibc_sa_t == 0 )  THEN
@@ -1253,5 +1260,5 @@
           WRITE ( io, 316 )
        ENDIF
-       IF ( ocean  .AND.  constant_top_salinityflux )                          &
+       IF ( ocean_mode  .AND.  constant_top_salinityflux )                          &
           WRITE ( io, 309 )  top_salinityflux
        IF ( humidity       )  WRITE ( io, 315 )
@@ -1412,5 +1419,5 @@
 !-- Initial salinity profile
 !-- Building output strings, starting with surface salinity
-    IF ( ocean )  THEN
+    IF ( ocean_mode )  THEN
        WRITE ( temperatures, '(F6.2)' )  sa_surface
        gradients = '------'
@@ -1895,5 +1902,5 @@
     WRITE ( io, 412 )  TRIM( reference_state )
     IF ( use_single_reference_value )  THEN
-       IF ( ocean )  THEN
+       IF ( ocean_mode )  THEN
           WRITE ( io, 413 )  prho_reference
        ELSE
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 3291)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! changes concerning modularization of ocean option
+! 
+! 3289 2018-09-28 10:23:58Z suehring
 ! Introduce module parameter for number of inflow profiles
 ! 
@@ -549,4 +552,7 @@
     USE netcdf_data_input_mod,                                                 &
         ONLY:  init_3d, netcdf_data_input_init_3d
+
+    USE ocean_mod,                                                             &
+        ONLY:  ocean_init, ocean_init_arrays
     
     USE particle_attributes,                                                   &
@@ -772,23 +778,4 @@
                  s_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
                  s_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-#endif
-    ENDIF
-
-    IF ( ocean )  THEN
-#if defined( __nopointer )
-       ALLOCATE( prho(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
-                 rho_ocean(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
-                 sa(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                            &
-                 sa_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
-                 tsa_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-#else
-       ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                        &
-                 rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                         &
-                 sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
-                 sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
-                 sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-       prho => prho_1
-       rho_ocean  => rho_1  ! routines calc_mean_profile and diffusion_e require
-                      ! density to be apointer
 #endif
     ENDIF
@@ -1010,9 +997,6 @@
        s => s_1;  s_p => s_2;  ts_m => s_3
     ENDIF    
-
-    IF ( ocean )  THEN
-       sa => sa_1;  sa_p => sa_2;  tsa_m => sa_3
-    ENDIF
 #endif
+
 !
 !-- Initialize arrays for turbulence closure
@@ -1022,24 +1006,11 @@
     CALL init_surface_arrays
 !
-!-- Allocate microphysics module arrays
-    IF ( bulk_cloud_model )  THEN
-       CALL bcm_init_arrays
-    ENDIF
-!
-!-- Allocate land surface model arrays
-    IF ( land_surface )  THEN
-       CALL lsm_init_arrays
-    ENDIF
-
-!
-!-- Allocate wind turbine model arrays
-    IF ( wind_turbine )  THEN
-       CALL wtm_init_arrays
-    ENDIF
-!
-!-- Allocate gust module arrays
-    IF ( gust_module_enabled )  THEN
-       CALL gust_init_arrays
-    ENDIF
+!-- Allocate arrays for other modules
+    IF ( bulk_cloud_model    )  CALL bcm_init_arrays
+    IF ( gust_module_enabled )  CALL gust_init_arrays
+    IF ( land_surface        )  CALL lsm_init_arrays
+    IF ( ocean_mode          )  CALL ocean_init_arrays
+    IF ( wind_turbine        )  CALL wtm_init_arrays
+    IF ( uv_exposure         )  CALL uvem_init_arrays
 
 !
@@ -1053,9 +1024,4 @@
     IF ( uv_exposure )  THEN
        CALL uvem_init
-    ENDIF
-!
-!-- Allocate uv exposure arrays
-    IF ( uv_exposure )  THEN
-       CALL uvem_init_arrays
     ENDIF
 
@@ -1230,11 +1196,4 @@
           ENDIF
 
-          IF ( ocean )  THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   sa(:,j,i) = sa_init
-                ENDDO
-             ENDDO
-          ENDIF
 !
 !--       Set velocity components at non-atmospheric / oceanic grid points to 
@@ -1380,11 +1339,4 @@
           ENDIF
 
-          IF ( ocean )  THEN
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   sa(:,j,i) = sa_init
-                ENDDO
-             ENDDO
-          ENDIF
 !
 !--       Compute initial temperature field and other constants used in case
@@ -1449,11 +1401,4 @@
        hom(:,1,7,:)  = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
 
-
-!
-!--    Store initial salinity profile
-       IF ( ocean )  THEN
-          hom(:,1,26,:)  = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 )
-       ENDIF
-
        IF ( humidity )  THEN
 !
@@ -1557,9 +1502,4 @@
        ENDIF       
 
-       IF ( ocean )  THEN
-          tsa_m = 0.0_wp
-          sa_p  = sa
-       ENDIF
-       
        CALL location_message( 'finished', .TRUE. )
 
@@ -1810,6 +1750,4 @@
        ENDIF
        IF ( passive_scalar )  s_p  = s
-       IF ( ocean          )  sa_p = sa
-
 !
 !--    Allthough tendency arrays are set in prognostic_equations, they have
@@ -1821,5 +1759,4 @@
        ENDIF
        IF ( passive_scalar )  ts_m  = 0.0_wp
-       IF ( ocean          )  tsa_m = 0.0_wp
 !
 !--    Initialize synthetic turbulence generator in case of restart.
@@ -2241,16 +2178,13 @@
     IF ( dt_dvrp /= 9999999.9_wp )  CALL init_dvrp
 
-    IF ( ocean )  THEN
-!
-!--    Initialize quantities needed for the ocean model
-       CALL init_ocean
-
-    ENDIF
-!
-!--    Initialize quantities for handling cloud physics
-!--    This routine must be called before lpm_init, because
-!--    otherwise, array d_exner, needed in data_output_dvrp (called by
-!--    lpm_init) is not defined.
-    IF ( .NOT. ocean )  THEN
+!
+!-- Initialize quantities needed for the ocean model
+    IF ( ocean_mode )  CALL ocean_init
+
+!
+!-- Initialize quantities for handling cloud physics.
+!-- This routine must be called before lpm_init, becaus otherwise,
+!-- array d_exner, needed in data_output_dvrp (called by lpm_init) is not defined.
+    IF ( .NOT. ocean_mode )  THEN
 
        ALLOCATE( hyp(nzb:nzt+1) )
@@ -2454,5 +2388,6 @@
     rdf_sc = 0.0_wp
     IF ( rayleigh_damping_factor /= 0.0_wp )  THEN
-       IF (  .NOT.  ocean )  THEN
+
+       IF (  .NOT.  ocean_mode )  THEN
           DO  k = nzb+1, nzt
              IF ( zu(k) >= rayleigh_damping_height )  THEN
@@ -2464,4 +2399,7 @@
           ENDDO
        ELSE
+!
+!--       In ocean mode, rayleigh damping is applied in the lower part of the
+!--       model domain
           DO  k = nzt, nzb+1, -1
              IF ( zu(k) <= rayleigh_damping_height )  THEN
@@ -2473,4 +2411,5 @@
           ENDDO
        ENDIF
+
     ENDIF
     IF ( scalar_rayleigh_damping )  rdf_sc = rdf
Index: palm/trunk/SOURCE/init_grid.f90
===================================================================
--- palm/trunk/SOURCE/init_grid.f90	(revision 3291)
+++ palm/trunk/SOURCE/init_grid.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! ocean renamed ocean_mode
+! 
+! 3241 2018-09-12 15:02:00Z raasch
 ! unused variables removed
 ! 
@@ -359,5 +362,5 @@
                dz_stretch_level_start, ibc_uv_b, message_string,               &
                momentum_advec, number_stretch_level_end,                       &
-               number_stretch_level_start,ocean, psolver, scalar_advec,        &
+               number_stretch_level_start, ocean_mode, psolver, scalar_advec,  &
                topography, use_surface_fluxes
          
@@ -509,5 +512,5 @@
 !
 !-- Define the vertical grid levels
-    IF ( .NOT. ocean )  THEN
+    IF ( .NOT. ocean_mode )  THEN
     
 !
@@ -1293,5 +1296,5 @@
 
     USE control_parameters,                                                    &
-        ONLY:  bc_lr_cyc, bc_ns_cyc, message_string, ocean
+        ONLY:  bc_lr_cyc, bc_ns_cyc, message_string, ocean_mode
 
     USE indices,                                                               &
@@ -1386,5 +1389,5 @@
 !-- In order to map topography on PALM grid also in case of ocean simulations,
 !-- pre-calculate an offset value.
-    ocean_offset = MERGE( zw(0), 0.0_wp, ocean )
+    ocean_offset = MERGE( zw(0), 0.0_wp, ocean_mode )
 !
 !-- Reference buildings on top of orography. This is not necessary
@@ -1638,5 +1641,5 @@
 !-- Topography input via ASCII format. 
     ELSE
-       ocean_offset     = MERGE( zw(0), 0.0_wp, ocean )
+       ocean_offset     = MERGE( zw(0), 0.0_wp, ocean_mode )
        topo_3d          = IBSET( topo_3d, 0 )
        topo_3d(nzb,:,:) = IBCLR( topo_3d(nzb,:,:), 0 )
Index: palm/trunk/SOURCE/init_masks.f90
===================================================================
--- palm/trunk/SOURCE/init_masks.f90	(revision 3291)
+++ palm/trunk/SOURCE/init_masks.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! ocean renamed ocean_mode
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -145,5 +148,5 @@
                mask_x_loop, mask_xyz_dimension, mask_y, mask_y_loop, mask_z,   &
                mask_z_loop, max_masks,  message_string, mid,                   &
-               passive_scalar, ocean, varnamelength
+               passive_scalar, ocean_mode, varnamelength
 
 
@@ -412,7 +415,7 @@
 
              CASE ( 'rho_ocean' )
-                IF ( .NOT. ocean )  THEN
-                   WRITE ( message_string, * ) 'output of "', TRIM( var ),     &
-                        '" requires ocean = .TRUE.'
+                IF ( .NOT. ocean_mode )  THEN
+                   WRITE ( message_string, * ) 'output of "', TRIM( var ),     &
+                        '" requires ocean mode'
                    CALL message( 'init_masks', 'PA0109', 1, 2, 0, 6, 0 )
                 ENDIF
@@ -428,7 +431,7 @@
 
              CASE ( 'sa' )
-                IF ( .NOT. ocean )  THEN
-                   WRITE ( message_string, * ) 'output of "', TRIM( var ),     &
-                        '" requires ocean = .TRUE.'
+                IF ( .NOT. ocean_mode )  THEN
+                   WRITE ( message_string, * ) 'output of "', TRIM( var ),     &
+                        '" requires ocean mode'
                    CALL message( 'init_masks', 'PA0109', 1, 2, 0, 6, 0 )
                 ENDIF
Index: palm/trunk/SOURCE/init_ocean.f90
===================================================================
--- palm/trunk/SOURCE/init_ocean.f90	(revision 3291)
+++ 	(revision )
@@ -1,213 +1,0 @@
-!> @file init_ocean.f90
-!------------------------------------------------------------------------------!
-! This file is part of the PALM model system.
-!
-! PALM is free software: you can redistribute it and/or modify it under the
-! terms of the GNU General Public License as published by the Free Software
-! Foundation, either version 3 of the License, or (at your option) any later
-! version.
-!
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
-! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along with
-! PALM. If not, see <http://www.gnu.org/licenses/>.
-!
-! Copyright 1997-2018 Leibniz Universitaet Hannover
-!------------------------------------------------------------------------------!
-!
-! Current revisions:
-! -----------------
-! 
-! 
-! Former revisions:
-! ------------------
-! $Id$
-! Modularization of all bulk cloud physics code components
-! 
-! 3270 2018-09-24 06:28:20Z raasch
-! bugfix for r2502: wrong calculation of hyp(nzt) in second iterative step
-! 
-! 2846 2018-03-01 08:48:47Z raasch
-! FORTRAN bugfix for r2845
-! 
-! 2845 2018-03-01 08:32:34Z raasch
-! bugfix: set kinematic viscosity for sea water
-! 
-! 2718 2018-01-02 08:49:38Z maronga
-! Corrected "Former revisions" section
-! 
-! 2696 2017-12-14 17:12:51Z kanani
-! Change in file header (GPL part)
-!
-! 2502 2017-09-26 15:51:59Z gronemeier
-! Bugfix: use equation of state for seawater to calculate rho_ocean_init at nzt
-!
-! 2195 2017-03-23 08:15:17Z raasch
-!
-! 2194 2017-03-23 08:03:04Z raasch
-! bugfix: density is now used in single-value reference state instead of
-! potential density
-!
-! 2031 2016-10-21 15:11:58Z knoop
-! renamed variable rho_init to rho_ocean_init
-! 
-! 2000 2016-08-20 18:09:15Z knoop
-! Forced header and separation lines into 80 columns
-! 
-! 1682 2015-10-07 23:56:08Z knoop
-! Code annotations made doxygen readable 
-! 
-! 1353 2014-04-08 15:21:23Z heinze
-! REAL constants provided with KIND-attribute
-!
-! 1322 2014-03-20 16:38:49Z raasch
-! REAL constants defined as wp_kind
-!
-! 1320 2014-03-20 08:40:49Z raasch
-! ONLY-attribute added to USE-statements,
-! kind-parameters added to all INTEGER and REAL declaration statements,
-! kinds are defined in new module kinds,
-! revision history before 2012 removed,
-! comment fields (!:) to be used for variable explanations added to
-! all variable declaration statements
-!
-! 1179 2013-06-14 05:57:58Z raasch
-! Initial density profile is stored in array hom
-!
-! 1036 2012-10-22 13:43:42Z raasch
-! code put under GPL (PALM 3.9)
-!
-! 97 2007-06-21 08:23:15Z raasch
-! Initial revision
-!
-! Description:
-! ------------
-!> Initialization of quantities needed for the ocean version
-!------------------------------------------------------------------------------!
- SUBROUTINE init_ocean
- 
-
-    USE arrays_3d,                                                             &
-        ONLY:  dzu, hyp, pt_init, ref_state, sa_init, zu, zw
-
-    USE basic_constants_and_equations_mod,                                     &
-        ONLY:  g
-
-    USE control_parameters,                                                    &
-        ONLY:  molecular_viscosity, prho_reference, rho_surface,               &
-               rho_reference, surface_pressure, use_single_reference_value
-
-    USE eqn_state_seawater_mod,                                                &
-        ONLY:  eqn_state_seawater, eqn_state_seawater_func
-
-    USE indices,                                                               &
-        ONLY:  nzb, nzt
-
-    USE kinds
-
-    USE pegrid
-
-    USE statistics,                                                            &
-        ONLY:  hom, statistic_regions
-
-    IMPLICIT NONE
-
-    INTEGER(iwp) ::  k !<
-    INTEGER(iwp) ::  n !<
-
-    REAL(wp)     ::  pt_l !<
-    REAL(wp)     ::  sa_l !<
-
-    REAL(wp), DIMENSION(nzb:nzt+1) ::  rho_ocean_init !<
-
-    ALLOCATE( hyp(nzb:nzt+1) )
-
-!
-!-- Set water density near the ocean surface
-    rho_surface = 1027.62_wp
-
-!
-!-- Set kinematic viscosity to sea water at 20C
-!-- WARNING: This value is especially used for calculating terminal fall
-!--          velocities in the LPM, so lpm_init should always be called after
-!--          init_ocean!
-    molecular_viscosity = 1.05E-6_wp
-
-!
-!-- Calculate initial vertical profile of hydrostatic pressure (in Pa)
-!-- and the reference density (used later in buoyancy term)
-!-- First step: Calculate pressure using reference density
-    hyp(nzt+1) = surface_pressure * 100.0_wp
-    hyp(nzt)   = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1)
-    rho_ocean_init(nzt)   = rho_surface
-    rho_ocean_init(nzt+1) = rho_surface  ! only required for output
-
-    DO  k = nzt-1, 1, -1
-       hyp(k) = hyp(k+1) + rho_surface * g * dzu(k)
-    ENDDO
-    hyp(0) = hyp(1) + rho_surface * g * dzu(1)
-
-!
-!-- Second step: Iteratively calculate in situ density (based on presssure)
-!-- and pressure (based on in situ density)
-    DO  n = 1, 5
-
-       rho_reference = rho_surface * 0.5_wp * dzu(nzt+1)
-
-       DO  k = nzt, 0, -1
-
-          sa_l = 0.5_wp * ( sa_init(k) + sa_init(k+1) )
-          pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) )
-
-          rho_ocean_init(k) = eqn_state_seawater_func( hyp(k), pt_l, sa_l )
-
-          rho_reference = rho_reference + rho_ocean_init(k) * dzu(k+1)
-
-       ENDDO
-
-       rho_reference = rho_reference / ( zw(nzt) - zu(nzb) )
-
-       hyp(nzt) = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1)
-       DO  k = nzt-1, 0, -1
-          hyp(k) = hyp(k+1) + g * 0.5_wp * ( rho_ocean_init(k)                 &
-                                           + rho_ocean_init(k+1) ) * dzu(k+1)
-       ENDDO
-
-    ENDDO
-
-!
-!-- Calculate the reference potential density
-    prho_reference = 0.0_wp
-    DO  k = 0, nzt
-
-       sa_l = 0.5_wp * ( sa_init(k) + sa_init(k+1) )
-       pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) )
-
-       prho_reference = prho_reference + dzu(k+1) * &
-                        eqn_state_seawater_func( 0.0_wp, pt_l, sa_l )
-
-    ENDDO
-
-    prho_reference = prho_reference / ( zu(nzt) - zu(nzb) )
-
-!
-!-- Calculate the 3d array of initial in situ and potential density,
-!-- based on the initial temperature and salinity profile
-    CALL eqn_state_seawater
-
-!
-!-- Store initial density profile
-    hom(:,1,77,:)  = SPREAD( rho_ocean_init(:), 2, statistic_regions+1 )
-
-!
-!-- Set the reference state to be used in the buoyancy terms
-    IF ( use_single_reference_value )  THEN
-       ref_state(:) = rho_reference
-    ELSE
-       ref_state(:) = rho_ocean_init(:)
-    ENDIF
-
-
- END SUBROUTINE init_ocean
Index: palm/trunk/SOURCE/init_vertical_profiles.f90
===================================================================
--- palm/trunk/SOURCE/init_vertical_profiles.f90	(revision 3294)
+++ palm/trunk/SOURCE/init_vertical_profiles.f90	(revision 3294)
@@ -0,0 +1,148 @@
+!> @file ocean_mod.f90
+!------------------------------------------------------------------------------!
+! This file is part of the PALM model system.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation, either version 3 of the License, or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 2017-2018 Leibniz Universitaet Hannover
+!--------------------------------------------------------------------------------!
+!
+! Current revisions:
+! -----------------
+! 
+! 
+! Former revisions:
+! -----------------
+! $Id$
+! split from check_parameters as separate file to avoid circular dependency
+! with ocean_mod
+!
+! 
+!
+!
+! Authors:
+! --------
+! @author Siegfried Raasch
+!
+! Description:
+! ------------
+!> Inititalizes the vertical profiles of scalar quantities.
+!------------------------------------------------------------------------------!
+ SUBROUTINE init_vertical_profiles( vertical_gradient_level_ind,               &
+                                    vertical_gradient_level,                   &
+                                    vertical_gradient, initial_profile,        &
+                                    surface_value, bc_top_gradient )
+
+    USE arrays_3d,                                                             &
+        ONLY:  dzu, zu
+
+    USE control_parameters,                                                    &
+        ONLY:  ocean_mode
+
+    USE indices,                                                               &
+        ONLY:  nz, nzt
+
+    USE kinds
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i  !< loop counter
+    INTEGER(iwp) ::  k  !< loop counter
+    INTEGER(iwp), DIMENSION(1:10) ::  vertical_gradient_level_ind  !< vertical grid indices for gradient levels
+
+    REAL(wp)     ::  bc_top_gradient  !< model top gradient
+    REAL(wp)     ::  gradient         !< vertica gradient of the respective quantity
+    REAL(wp)     ::  surface_value    !< surface value of the respecitve quantity
+
+    REAL(wp), DIMENSION(0:nz+1) ::  initial_profile          !< initialisation profile
+    REAL(wp), DIMENSION(1:10)   ::  vertical_gradient        !< given vertical gradient
+    REAL(wp), DIMENSION(1:10)   ::  vertical_gradient_level  !< given vertical gradient level
+
+    i = 1
+    gradient = 0.0_wp
+
+    IF ( .NOT. ocean_mode )  THEN
+
+       vertical_gradient_level_ind(1) = 0
+       DO  k = 1, nzt+1
+          IF ( i < 11 )  THEN
+             IF ( vertical_gradient_level(i) < zu(k)  .AND.            &
+                  vertical_gradient_level(i) >= 0.0_wp )  THEN
+                gradient = vertical_gradient(i) / 100.0_wp
+                vertical_gradient_level_ind(i) = k - 1
+                i = i + 1
+             ENDIF
+          ENDIF
+          IF ( gradient /= 0.0_wp )  THEN
+             IF ( k /= 1 )  THEN
+                initial_profile(k) = initial_profile(k-1) + dzu(k) * gradient
+             ELSE
+                initial_profile(k) = initial_profile(k-1) + dzu(k) * gradient
+             ENDIF
+          ELSE
+             initial_profile(k) = initial_profile(k-1)
+          ENDIF
+!
+!--       Avoid negative values of scalars
+          IF ( initial_profile(k) < 0.0_wp )  THEN
+             initial_profile(k) = 0.0_wp
+          ENDIF
+       ENDDO
+
+    ELSE
+
+!
+!--    In ocean mode, profiles are constructed starting from the ocean surface,
+!--    which is at the top of the model domain
+       vertical_gradient_level_ind(1) = nzt+1
+       DO  k = nzt, 0, -1
+          IF ( i < 11 )  THEN
+             IF ( vertical_gradient_level(i) > zu(k)  .AND.            &
+                  vertical_gradient_level(i) <= 0.0_wp )  THEN
+                gradient = vertical_gradient(i) / 100.0_wp
+                vertical_gradient_level_ind(i) = k + 1
+                i = i + 1
+             ENDIF
+          ENDIF
+          IF ( gradient /= 0.0_wp )  THEN
+             IF ( k /= nzt )  THEN
+                initial_profile(k) = initial_profile(k+1) - dzu(k+1) * gradient
+             ELSE
+                initial_profile(k)   = surface_value - 0.5_wp * dzu(k+1) *     &
+                                                       gradient
+                initial_profile(k+1) = surface_value + 0.5_wp * dzu(k+1) *     &
+                                                       gradient
+             ENDIF
+          ELSE
+             initial_profile(k) = initial_profile(k+1)
+          ENDIF
+!
+!--       Avoid negative values of scalars
+          IF ( initial_profile(k) < 0.0_wp )  THEN
+             initial_profile(k) = 0.0_wp
+          ENDIF
+       ENDDO
+
+    ENDIF
+
+!
+!-- In case of no given gradients, choose zero gradient conditions
+    IF ( vertical_gradient_level(1) == -999999.9_wp )  THEN
+       vertical_gradient_level(1) = 0.0_wp
+    ENDIF
+!
+!-- Store gradient at the top boundary for possible Neumann boundary condition
+    bc_top_gradient  = ( initial_profile(nzt+1) - initial_profile(nzt) ) /     &
+                       dzu(nzt+1)
+
+ END SUBROUTINE init_vertical_profiles
Index: palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90
===================================================================
--- palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90	(revision 3291)
+++ palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! ocean renamed ocean_mode
+! 
+! 3241 2018-09-12 15:02:00Z raasch
 ! unused variables removed
 ! 
@@ -92,5 +95,5 @@
                large_scale_forcing, large_scale_subsidence, lsf_surf, lsf_vert,&
                lsf_exception, message_string, nesting_offline, neutral,        &
-               nudging, passive_scalar, pt_surface, ocean, q_surface,          &
+               nudging, passive_scalar, pt_surface, ocean_mode, q_surface,     &
                surface_heatflux, surface_pressure, surface_waterflux,          &
                topography, use_subsidence_tendencies
@@ -696,7 +699,7 @@
        ENDIF
 
-       IF ( large_scale_forcing  .AND.  ocean )  THEN
+       IF ( large_scale_forcing  .AND.  ocean_mode )  THEN
           message_string = 'The usage of large scale forcing from external &'//&
-                        'file LSF_DATA is not implemented for ocean runs'
+                        'file LSF_DATA is not implemented for ocean mode'
           CALL message( 'check_parameters', 'PA0378', 1, 2, 0, 6, 0 )
        ENDIF
Index: palm/trunk/SOURCE/lpm_init.f90
===================================================================
--- palm/trunk/SOURCE/lpm_init.f90	(revision 3291)
+++ palm/trunk/SOURCE/lpm_init.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! ocean renamed ocean_mode
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -231,5 +234,5 @@
     USE control_parameters,                                                    &
         ONLY:  cloud_droplets, constant_flux_layer, current_timestep_number,   &
-               dt_3d, dz, initializing_actions, message_string, ocean,         &
+               dt_3d, dz, initializing_actions, message_string, ocean_mode,    &
                simulated_time
 
@@ -328,5 +331,5 @@
 !-- In case of oceans runs, the vertical index calculations need an offset,
 !-- because otherwise the k indices will become negative
-    IF ( ocean )  THEN
+    IF ( ocean_mode )  THEN
        offset_ocean_nzt    = nzt
        offset_ocean_nzt_m1 = nzt - 1
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 3291)
+++ palm/trunk/SOURCE/modules.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! ocean renamed ocean_mode
+! 
+! 3289 2018-09-28 10:23:58Z suehring
 ! +num_mean_inflow_profiles
 ! 
@@ -1109,5 +1112,4 @@
     CHARACTER (LEN=20)   ::  bc_s_b = 'dirichlet'                         !< namelist parameter
     CHARACTER (LEN=20)   ::  bc_s_t = 'initial_gradient'                  !< namelist parameter
-    CHARACTER (LEN=20)   ::  bc_sa_t = 'neumann'                          !< namelist parameter
     CHARACTER (LEN=20)   ::  bc_uv_b = 'dirichlet'                        !< namelist parameter
     CHARACTER (LEN=20)   ::  bc_uv_t = 'dirichlet'                        !< namelist parameter
@@ -1180,5 +1182,4 @@
     INTEGER(iwp) ::  ibc_s_b                           !< integer flag for bc_s_b
     INTEGER(iwp) ::  ibc_s_t                           !< integer flag for bc_s_t
-    INTEGER(iwp) ::  ibc_sa_t                          !< integer flag for bc_sa_t
     INTEGER(iwp) ::  ibc_uv_b                          !< integer flag for bc_uv_b
     INTEGER(iwp) ::  ibc_uv_t                          !< integer flag for bc_uv_t
@@ -1242,5 +1243,4 @@
     INTEGER(iwp) ::  q_vertical_gradient_level_ind(10) = -9999   !< grid index values of q_vertical_gradient_level(s)
     INTEGER(iwp) ::  s_vertical_gradient_level_ind(10) = -9999   !< grid index values of s_vertical_gradient_level(s)   
-    INTEGER(iwp) ::  sa_vertical_gradient_level_ind(10) = -9999  !< grid index values of sa_vertical_gradient_level(s)
     INTEGER(iwp) ::  section(100,3)                              !< collective array for section_xy/xz/yz
     INTEGER(iwp) ::  section_xy(100) = -9999                     !< namelist parameter
@@ -1293,5 +1293,5 @@
     LOGICAL ::  constant_top_heatflux = .TRUE.                   !< heat flux at domain top constant?
     LOGICAL ::  constant_top_momentumflux = .FALSE.              !< momentum flux at domain topconstant?
-    LOGICAL ::  constant_top_salinityflux = .TRUE.               !< salinity flux at ocean domain top?
+    LOGICAL ::  constant_top_salinityflux = .TRUE.               !< constant salinity flux at ocean surface
     LOGICAL ::  constant_top_scalarflux = .TRUE.                 !< passive-scalar flux at domain top constant?
     LOGICAL ::  constant_scalarflux = .TRUE.                     !< passive-scalar flux at surfaces constant?
@@ -1328,5 +1328,5 @@
     LOGICAL ::  neutral = .FALSE.                                !< namelist parameter
     LOGICAL ::  nudging = .FALSE.                                !< namelist parameter
-    LOGICAL ::  ocean = .FALSE.                                  !< namelist parameter
+    LOGICAL ::  ocean_mode = .FALSE.                             !< namelist parameter
     LOGICAL ::  passive_scalar = .FALSE.                         !< namelist parameter
     LOGICAL ::  plant_canopy = .FALSE.                           !< switch for use of plant canopy model
@@ -1380,5 +1380,5 @@
                                                                !< (galilei transformation)
     REAL(wp) ::  alpha_surface = 0.0_wp                        !< namelist parameter
-    REAL(wp) ::  atmos_ocean_sign = 1.0_wp                     !< vertical-grid conversion factor 
+    REAL(wp) ::  atmos_ocean_sign = 1.0_wp                     !< vertical-grid conversion factor
                                                                !< (=1.0 in atmosphere, =-1.0 in ocean)
     REAL(wp) ::  averaging_interval = 0.0_wp                   !< namelist parameter
@@ -1454,5 +1454,4 @@
     REAL(wp) ::  particle_maximum_age = 9999999.9_wp           !< namelist parameter
     REAL(wp) ::  prandtl_number = 1.0_wp                       !< namelist parameter
-    REAL(wp) ::  prho_reference                                !< reference state of potential density
     REAL(wp) ::  pt_damping_factor = 0.0_wp                    !< namelist parameter
     REAL(wp) ::  pt_damping_width = 0.0_wp                     !< namelist parameter
@@ -1472,5 +1471,4 @@
     REAL(wp) ::  rho_surface                                   !< surface value of density
     REAL(wp) ::  roughness_length = 0.1_wp                     !< namelist parameter
-    REAL(wp) ::  sa_surface = 35.0_wp                          !< namelist parameter
     REAL(wp) ::  simulated_time = 0.0_wp                       !< elapsed simulated time
     REAL(wp) ::  simulated_time_at_begin                       !< elapsed simulated time of previous run (job chain)
@@ -1550,6 +1548,4 @@
     REAL(wp) ::  s_vertical_gradient(10) = 0.0_wp                  !< namelist parameter
     REAL(wp) ::  s_vertical_gradient_level(10) = -999999.9_wp    !< namelist parameter
-    REAL(wp) ::  sa_vertical_gradient(10) = 0.0_wp                 !< namelist parameter
-    REAL(wp) ::  sa_vertical_gradient_level(10) = -999999.9_wp   !< namelist parameter
     REAL(wp) ::  skip_time_domask(max_masks) = 9999999.9_wp        !< namelist parameter
     REAL(wp) ::  threshold(20) = 0.0_wp                            !< namelist parameter
Index: palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3291)
+++ palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! changes concerning modularization of ocean option
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -568,6 +571,6 @@
                do3d_time_count, domask_time_count, end_time, land_surface,     &
                mask_size_l, mask_i, mask_i_global, mask_j, mask_j_global,      &
-               mask_k_global, message_string, mid, ntdim_2d_xy,                &
-               ntdim_2d_xz, ntdim_2d_yz, ntdim_3d, nz_do3d, plant_canopy,      &
+               mask_k_global, message_string, mid, ntdim_2d_xy, ntdim_2d_xz,   &
+               ntdim_2d_yz, ntdim_3d, nz_do3d, ocean_mode, plant_canopy,       &
                run_description_header, section, simulated_time,                &
                simulated_time_at_begin, skip_time_data_output_av,              &
@@ -589,4 +592,7 @@
     USE land_surface_model_mod,                                                &
         ONLY: lsm_define_netcdf_grid, nzb_soil, nzt_soil, nzs, zs
+
+    USE ocean_mod,                                                             &
+        ONLY:  ocean_define_netcdf_grid
 
     USE pegrid
@@ -943,5 +949,5 @@
                 CASE ( 'e', 'lpt', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'pt',   &
                        'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
-                       'rho_ocean', 's', 'sa', 'vpt'  )
+                       's', 'vpt'  )
 
                    grid_x = 'x'
@@ -971,23 +977,33 @@
 
 !             
-!--       Block of urban surface model outputs   
+!--             Block of urban surface model outputs
                 CASE ( 'usm_output' )
 
-                   CALL usm_define_netcdf_grid( domask(mid,av,i), found, &
+                   CALL usm_define_netcdf_grid( domask( mid,av,i), found,      &
                                                         grid_x, grid_y, grid_z )
 
                 CASE DEFAULT
-
-                   CALL tcm_define_netcdf_grid( domask(mid,av,i), found, &
-                                                        grid_x, grid_y, grid_z )
-
-!
-!--                Check for land surface quantities
+!
+!--                Check for quantities defined in other modules
+                   IF ( .NOT. found  .AND.  air_chemistry )  THEN
+                      CALL chem_define_netcdf_grid( domask(mid,av,i), found,   &
+                                                    grid_x, grid_y, grid_z )
+                   ENDIF
+
+                   IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
+                      CALL gust_define_netcdf_grid( domask(mid,av,i), found,   &
+                                                    grid_x, grid_y, grid_z )
+                   ENDIF
+
                    IF ( land_surface )  THEN
                       CALL lsm_define_netcdf_grid( domask(mid,av,i), found,    &
                                                    grid_x, grid_y, grid_z )
                    ENDIF
-!
-!--                Check for plant canopy quantities
+
+                   IF ( .NOT. found  .AND.  ocean_mode )  THEN
+                      CALL ocean_define_netcdf_grid( domask(mid,av,i), found,  &
+                                                     grid_x, grid_y, grid_z )
+                   ENDIF
+
                    IF ( .NOT. found  .AND.  plant_canopy )  THEN
                       CALL pcm_define_netcdf_grid( domask(mid,av,i), found,    &
@@ -995,6 +1011,4 @@
                    ENDIF
 
-!
-!--                Check for radiation quantities
                    IF ( .NOT. found  .AND.  radiation )  THEN
                       CALL radiation_define_netcdf_grid( domask(mid,av,i),     &
@@ -1003,21 +1017,9 @@
                    ENDIF
 
-!
-!--                Check for gust module quantities
-                   IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
-                      CALL gust_define_netcdf_grid( domask(mid,av,i), found,   &
-                                                    grid_x, grid_y, grid_z )
-                   ENDIF
-
-!
-!--                Check for chemistry quantities                   
-                   IF ( .NOT. found  .AND.  air_chemistry )  THEN
-                      CALL chem_define_netcdf_grid( domask(mid,av,i),          &
-                                                    found, grid_x, grid_y,     &
-                                                    grid_z )
-                   ENDIF
-
-!
-!--                Check for user-defined quantities
+                   CALL tcm_define_netcdf_grid( domask( mid,av,i), found,      &
+                                                        grid_x, grid_y, grid_z )
+
+!
+!--                Now check for user-defined quantities
                    IF ( .NOT. found )  THEN
                       CALL user_define_netcdf_grid( domask(mid,av,i), found,   &
@@ -4806,5 +4808,5 @@
 !
 !--       Define y-axis
-          CALL netcdf_create_dim(id_set_sp, 'k_y', ny/2, id_dim_y_sp, 282 )
+          CALL netcdf_create_dim( id_set_sp, 'k_y', ny/2, id_dim_y_sp, 282 )
           CALL netcdf_create_var( id_set_sp, (/ id_dim_y_sp /), 'k_y',         &
                                   NF90_DOUBLE, id_var_y_sp, 'm-1', '', 283,    &
Index: palm/trunk/SOURCE/ocean_mod.f90
===================================================================
--- palm/trunk/SOURCE/ocean_mod.f90	(revision 3294)
+++ palm/trunk/SOURCE/ocean_mod.f90	(revision 3294)
@@ -0,0 +1,1794 @@
+!> @file ocean_mod.f90
+!------------------------------------------------------------------------------!
+! This file is part of the PALM model system.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation, either version 3 of the License, or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 2017-2018 Leibniz Universitaet Hannover
+!--------------------------------------------------------------------------------!
+!
+! Current revisions:
+! -----------------
+! 
+! 
+! Former revisions:
+! -----------------
+! $Id$
+! initial revision
+! 
+! 
+! 
+!
+! Authors:
+! --------
+! @author Siegfried Raasch
+!
+! Description:
+! ------------
+!> This module contains relevant code for PALM's ocean option, e.g. the
+!> prognostic equation for salinity, the equation of state for seawater,
+!> and the Craik Leibovich force (Stokes force)
+!>
+!------------------------------------------------------------------------------!
+ MODULE ocean_mod
+ 
+
+#if defined( __nopointer )
+    USE arrays_3d,                                                             &
+        ONLY:  prho, rho_ocean, sa, sa_init, sa_p, tsa_m
+#else
+    USE arrays_3d,                                                             &
+        ONLY:  prho, prho_1, rho_ocean, rho_1, sa, sa_init, sa_1, sa_2, sa_3,  &
+               sa_p, tsa_m
+#endif
+
+    USE control_parameters,                                                    &
+        ONLY:  atmos_ocean_sign, bottom_salinityflux,                          &
+               constant_top_salinityflux, ocean_mode, top_salinityflux,        &
+               wall_salinityflux
+
+    USE kinds
+
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=20) ::  bc_sa_t = 'neumann'  !< namelist parameter
+
+    INTEGER(iwp) ::  ibc_sa_t   !< integer flag for bc_sa_t
+
+    INTEGER(iwp) ::  sa_vertical_gradient_level_ind(10) = -9999  !< grid index values of sa_vertical_gradient_level(s)
+
+    LOGICAL ::  stokes_force = .FALSE.              !< switch to switch on the Stokes force
+
+    REAL(wp) ::  langmuir_number = 0.34_wp    !< turbulent Langmuir number, default from Li et al., 2005
+    REAL(wp) ::  prho_reference               !< reference state of potential density at ocean surface
+    REAL(wp) ::  sa_surface = 35.0_wp         !< surface salinity, namelist parameter
+    REAL(wp) ::  sa_vertical_gradient(10) = 0.0_wp               !< namelist parameter
+    REAL(wp) ::  sa_vertical_gradient_level(10) = -999999.9_wp   !< namelist parameter
+    REAL(wp) ::  stokes_waveheight =  1.0_wp  !< typical wave height assumed for Stokes velocity
+    REAL(wp) ::  stokes_wavelength = 40.0_wp  !< typical wavelength assumed for Stokes velocity
+
+    REAL(wp), 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 /)
+                          !< constants used in equation of state for seawater
+
+    REAL(wp), DIMENSION(13), PARAMETER ::  den =                               &
+                          (/ 1.0D0,                  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 /)
+                          !< constants used in equation of state for seawater
+
+    SAVE
+
+    PUBLIC ::  bc_sa_t, ibc_sa_t, langmuir_number, prho_reference, sa_surface, &
+               sa_vertical_gradient, sa_vertical_gradient_level,               &
+               sa_vertical_gradient_level_ind, stokes_force,                   &
+               stokes_waveheight, stokes_wavelength
+
+
+    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
+
+    INTERFACE ocean_check_parameters
+       MODULE PROCEDURE ocean_check_parameters
+    END INTERFACE ocean_check_parameters
+
+    INTERFACE ocean_check_data_output
+       MODULE PROCEDURE ocean_check_data_output
+    END INTERFACE ocean_check_data_output
+
+    INTERFACE ocean_check_data_output_pr
+       MODULE PROCEDURE ocean_check_data_output_pr
+    END INTERFACE ocean_check_data_output_pr
+
+    INTERFACE ocean_define_netcdf_grid
+       MODULE PROCEDURE ocean_define_netcdf_grid
+    END INTERFACE ocean_define_netcdf_grid
+
+    INTERFACE ocean_data_output_2d
+       MODULE PROCEDURE ocean_data_output_2d
+    END INTERFACE ocean_data_output_2d
+
+    INTERFACE ocean_data_output_3d
+       MODULE PROCEDURE ocean_data_output_3d
+    END INTERFACE ocean_data_output_3d
+
+    INTERFACE ocean_init
+       MODULE PROCEDURE ocean_init
+    END INTERFACE ocean_init
+
+    INTERFACE ocean_init_arrays
+       MODULE PROCEDURE ocean_init_arrays
+    END INTERFACE ocean_init_arrays
+
+    INTERFACE ocean_parin
+       MODULE PROCEDURE ocean_parin
+    END INTERFACE ocean_parin
+
+    INTERFACE ocean_prognostic_equations
+       MODULE PROCEDURE ocean_prognostic_equations
+       MODULE PROCEDURE ocean_prognostic_equations_ij
+    END INTERFACE ocean_prognostic_equations
+
+    INTERFACE ocean_swap_timelevel
+       MODULE PROCEDURE ocean_swap_timelevel
+    END INTERFACE ocean_swap_timelevel
+
+    INTERFACE ocean_rrd_global
+       MODULE PROCEDURE ocean_rrd_global
+    END INTERFACE ocean_rrd_global
+
+    INTERFACE ocean_rrd_local
+       MODULE PROCEDURE ocean_rrd_local
+    END INTERFACE ocean_rrd_local
+
+    INTERFACE ocean_wrd_global
+       MODULE PROCEDURE ocean_wrd_global
+    END INTERFACE ocean_wrd_global
+
+    INTERFACE ocean_wrd_local
+       MODULE PROCEDURE ocean_wrd_local
+    END INTERFACE ocean_wrd_local
+
+    INTERFACE ocean_3d_data_averaging
+       MODULE PROCEDURE ocean_3d_data_averaging
+    END INTERFACE ocean_3d_data_averaging
+
+    SAVE
+
+    PRIVATE
+!
+!-- Add INTERFACES that must be available to other modules (alphabetical order)
+    PUBLIC eqn_state_seawater, ocean_check_data_output,                        &
+           ocean_check_data_output_pr, ocean_check_parameters,                 &
+           ocean_data_output_2d, ocean_data_output_3d,                         &
+           ocean_define_netcdf_grid, ocean_init, ocean_init_arrays,            &
+           ocean_parin, ocean_prognostic_equations, ocean_swap_timelevel,      &
+           ocean_rrd_global, ocean_rrd_local, ocean_wrd_global,                &
+           ocean_wrd_local, ocean_3d_data_averaging
+
+
+ CONTAINS
+
+!------------------------------------------------------------------------------!
+! 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.
+!> eqn_state_seawater calculates the potential density referred at hyp(0).
+!> eqn_state_seawater_func calculates density.
+!> TODO: so far, routine is not adjusted to use topography
+!------------------------------------------------------------------------------!
+ SUBROUTINE eqn_state_seawater
+
+    USE arrays_3d,                                                             &
+        ONLY:  hyp, prho, pt_p, rho_ocean, sa_p
+    USE indices,                                                               &
+        ONLY:  nxl, nxr, nyn, nys, nzb, nzt
+
+    USE surface_mod,                                                           &
+        ONLY :  bc_h
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i       !< running index x direction
+    INTEGER(iwp) ::  j       !< running index y direction
+    INTEGER(iwp) ::  k       !< running index z direction
+    INTEGER(iwp) ::  m       !< running index surface elements
+
+    REAL(wp) ::  pden   !< temporary scalar user for calculating density
+    REAL(wp) ::  pnom   !< temporary scalar user for calculating density
+    REAL(wp) ::  p1     !< temporary scalar user for calculating density
+    REAL(wp) ::  p2     !< temporary scalar user for calculating density
+    REAL(wp) ::  p3     !< temporary scalar user for calculating density
+    REAL(wp) ::  pt1    !< temporary scalar user for calculating density
+    REAL(wp) ::  pt2    !< temporary scalar user for calculating density
+    REAL(wp) ::  pt3    !< temporary scalar user for calculating density
+    REAL(wp) ::  pt4    !< temporary scalar user for calculating density
+    REAL(wp) ::  sa1    !< temporary scalar user for calculating density
+    REAL(wp) ::  sa15   !< temporary scalar user for calculating density
+    REAL(wp) ::  sa2    !< temporary scalar user for calculating density
+
+
+    DO  i = nxl, nxr
+       DO  j = nys, nyn
+          DO  k = nzb+1, nzt
+!
+!--          Pressure is needed in dbar
+             p1 = hyp(k) * 1E-4_wp
+             p2 = p1 * p1
+             p3 = p2 * p1
+
+!
+!--          Temperature needed in degree Celsius
+             pt1 = pt_p(k,j,i) - 273.15_wp
+             pt2 = pt1 * pt1
+             pt3 = pt1 * pt2
+             pt4 = pt2 * pt2
+
+             sa1  = sa_p(k,j,i)
+             sa15 = sa1 * SQRT( sa1 )
+             sa2  = sa1 * sa1
+
+             pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     +       &
+                    nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 +       &
+                    nom(7)*sa2
+
+             pden = 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
+!
+!--          Potential density (without pressure terms)
+             prho(k,j,i) = pnom / pden
+
+             pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  +       &
+                    nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2
+
+             pden = pden +             den(11)*p1     + den(12)*p2*pt3 +       &
+                    den(13)*p3*pt1
+
+!
+!--          In-situ density
+             rho_ocean(k,j,i) = pnom / pden
+
+          ENDDO
+
+!
+!--       Neumann conditions are assumed at top boundary and currently also at
+!--       bottom boundary (see comment lines below)
+          prho(nzt+1,j,i)      = prho(nzt,j,i)
+          rho_ocean(nzt+1,j,i) = rho_ocean(nzt,j,i)
+
+       ENDDO
+    ENDDO
+!
+!-- Neumann conditions at up/downward-facing surfaces
+    !$OMP PARALLEL DO PRIVATE( i, j, k )
+    DO  m = 1, bc_h(0)%ns
+       i = bc_h(0)%i(m)
+       j = bc_h(0)%j(m)
+       k = bc_h(0)%k(m)
+       prho(k-1,j,i)      = prho(k,j,i)
+       rho_ocean(k-1,j,i) = rho_ocean(k,j,i)
+    ENDDO
+!
+!-- Downward facing surfaces
+    !$OMP PARALLEL DO PRIVATE( i, j, k )
+    DO  m = 1, bc_h(1)%ns
+       i = bc_h(1)%i(m)
+       j = bc_h(1)%j(m)
+       k = bc_h(1)%k(m)
+       prho(k+1,j,i)      = prho(k,j,i)
+       rho_ocean(k+1,j,i) = rho_ocean(k,j,i)
+    ENDDO
+
+ END SUBROUTINE eqn_state_seawater
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Same as above, but for grid point i,j
+!------------------------------------------------------------------------------!
+ SUBROUTINE eqn_state_seawater_ij( i, j )
+
+    USE arrays_3d,                                                             &
+        ONLY:  hyp, prho, pt_p, rho_ocean, sa_p
+
+    USE indices,                                                               &
+        ONLY:  nzb, nzt
+
+    USE surface_mod,                                                           &
+       ONLY :  bc_h
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i       !< running index x direction
+    INTEGER(iwp) ::  j       !< running index y direction
+    INTEGER(iwp) ::  k       !< running index z direction
+    INTEGER(iwp) ::  m       !< running index surface elements
+    INTEGER(iwp) ::  surf_e  !< end index of surface elements at (j,i)-gridpoint
+    INTEGER(iwp) ::  surf_s  !< start index of surface elements at (j,i)-gridpoint
+
+    REAL(wp) ::  pden   !< temporary scalar user for calculating density
+    REAL(wp) ::  pnom   !< temporary scalar user for calculating density
+    REAL(wp) ::  p1     !< temporary scalar user for calculating density
+    REAL(wp) ::  p2     !< temporary scalar user for calculating density
+    REAL(wp) ::  p3     !< temporary scalar user for calculating density
+    REAL(wp) ::  pt1    !< temporary scalar user for calculating density
+    REAL(wp) ::  pt2    !< temporary scalar user for calculating density
+    REAL(wp) ::  pt3    !< temporary scalar user for calculating density
+    REAL(wp) ::  pt4    !< temporary scalar user for calculating density
+    REAL(wp) ::  sa1    !< temporary scalar user for calculating density
+    REAL(wp) ::  sa15   !< temporary scalar user for calculating density
+    REAL(wp) ::  sa2    !< temporary scalar user for calculating density
+
+    DO  k = nzb+1, nzt
+!
+!--    Pressure is needed in dbar
+       p1 = hyp(k) * 1E-4_wp
+       p2 = p1 * p1
+       p3 = p2 * p1
+!
+!--    Temperature needed in degree Celsius
+       pt1 = pt_p(k,j,i) - 273.15_wp
+       pt2 = pt1 * pt1
+       pt3 = pt1 * pt2
+       pt4 = pt2 * pt2
+
+       sa1  = sa_p(k,j,i)
+       sa15 = sa1 * SQRT( sa1 )
+       sa2  = sa1 * sa1
+
+       pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     +             &
+              nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + nom(7)*sa2
+
+       pden = 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
+!
+!--    Potential density (without pressure terms)
+       prho(k,j,i) = pnom / pden
+
+       pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  +             &
+              nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2
+       pden = pden +             den(11)*p1     + den(12)*p2*pt3 +             &
+              den(13)*p3*pt1
+
+!
+!--    In-situ density
+       rho_ocean(k,j,i) = pnom / pden
+
+    ENDDO
+!
+!-- Neumann conditions at up/downward-facing walls
+    surf_s = bc_h(0)%start_index(j,i)
+    surf_e = bc_h(0)%end_index(j,i)
+    DO  m = surf_s, surf_e
+       k = bc_h(0)%k(m)
+       prho(k-1,j,i)      = prho(k,j,i)
+       rho_ocean(k-1,j,i) = rho_ocean(k,j,i)
+    ENDDO
+!
+!-- Downward facing surfaces
+    surf_s = bc_h(1)%start_index(j,i)
+    surf_e = bc_h(1)%end_index(j,i)
+    DO  m = surf_s, surf_e
+       k = bc_h(1)%k(m)
+       prho(k+1,j,i)      = prho(k,j,i)
+       rho_ocean(k+1,j,i) = rho_ocean(k,j,i)
+    ENDDO
+!
+!-- Neumann condition are assumed at top boundary
+    prho(nzt+1,j,i)      = prho(nzt,j,i)
+    rho_ocean(nzt+1,j,i) = rho_ocean(nzt,j,i)
+
+ END SUBROUTINE eqn_state_seawater_ij
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Equation of state for seawater as a function
+!------------------------------------------------------------------------------!
+ REAL(wp) FUNCTION eqn_state_seawater_func( p, pt, sa )
+
+    IMPLICIT NONE
+
+    REAL(wp) ::  p      !< temporary scalar user for calculating density
+    REAL(wp) ::  p1     !< temporary scalar user for calculating density
+    REAL(wp) ::  p2     !< temporary scalar user for calculating density
+    REAL(wp) ::  p3     !< temporary scalar user for calculating density
+    REAL(wp) ::  pt     !< temporary scalar user for calculating density
+    REAL(wp) ::  pt1    !< temporary scalar user for calculating density
+    REAL(wp) ::  pt2    !< temporary scalar user for calculating density
+    REAL(wp) ::  pt3    !< temporary scalar user for calculating density
+    REAL(wp) ::  pt4    !< temporary scalar user for calculating density
+    REAL(wp) ::  sa     !< temporary scalar user for calculating density
+    REAL(wp) ::  sa15   !< temporary scalar user for calculating density
+    REAL(wp) ::  sa2    !< temporary scalar user for calculating density
+
+!
+!-- Pressure is needed in dbar
+    p1 = p  * 1.0E-4_wp
+    p2 = p1 * p1
+    p3 = p2 * p1
+
+!
+!-- Temperature needed in degree Celsius
+    pt1 = pt - 273.15_wp
+    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
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Reads the ocean parameters namelist
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_parin
+
+    USE control_parameters,                                                    &
+        ONLY:  message_string
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
+
+
+    NAMELIST /ocean_parameters/  bc_sa_t, bottom_salinityflux,                 &
+             langmuir_number, sa_surface, sa_vertical_gradient,                &
+             sa_vertical_gradient_level, stokes_force, stokes_waveheight,      &
+             stokes_wavelength, top_salinityflux, wall_salinityflux
+
+!
+!-- Try to find the namelist
+    REWIND ( 11 )
+    line = ' '
+    DO WHILE ( INDEX( line, '&ocean_parameters' ) == 0 )
+       READ ( 11, '(A)', END=12 )  line
+    ENDDO
+    BACKSPACE ( 11 )
+
+!
+!-- Read namelist
+    READ ( 11, ocean_parameters, ERR = 10 )
+!
+!-- Set switch that enables PALM's ocean mode
+    ocean_mode = .TRUE.
+
+    GOTO 12
+
+ 10 BACKSPACE( 11 )
+    READ( 11 , '(A)') line
+    CALL parin_fail_message( 'ocean_parameters', line )
+
+ 12 CONTINUE
+
+ END SUBROUTINE ocean_parin
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check parameters routine for the ocean mode
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_check_parameters
+
+    USE control_parameters,                                                    &
+        ONLY:  coupling_char, coupling_mode, message_string, use_top_fluxes
+
+    IMPLICIT NONE
+
+
+!
+!-- Check ocean setting
+    IF ( TRIM( coupling_mode ) == 'uncoupled'  .AND.                           &
+         TRIM( coupling_char ) == '_O' .AND.                                   &
+         .NOT. ocean_mode )  THEN
+
+!
+!--    Check whether an (uncoupled) atmospheric run has been declared as an
+!--    ocean run (this setting is done via palmrun-option -y)
+       message_string = 'ocean = .F. does not allow coupling_char = "' //      &
+                        TRIM( coupling_char ) // '" set by palmrun-option "-y"'
+       CALL message( 'check_parameters', 'PA0317', 1, 2, 0, 6, 0 )
+
+    ENDIF
+
+!
+!-- Ocean version must use flux boundary conditions at the top
+    IF ( .NOT. use_top_fluxes )  THEN
+       message_string = 'use_top_fluxes must be .TRUE. in ocean mode'
+       CALL message( 'ocean_check_parameters', 'PA0042', 1, 2, 0, 6, 0 )
+    ENDIF
+
+!
+!-- Boundary conditions for salinity
+    IF ( bc_sa_t == 'dirichlet' )  THEN
+       ibc_sa_t = 0
+    ELSEIF ( bc_sa_t == 'neumann' )  THEN
+       ibc_sa_t = 1
+    ELSE
+       message_string = 'unknown boundary condition: bc_sa_t = "' //           &
+                        TRIM( bc_sa_t ) // '"'
+       CALL message( 'ocean_check_parameters', 'PA0068', 1, 2, 0, 6, 0 )
+    ENDIF
+
+    IF ( top_salinityflux == 9999999.9_wp )  constant_top_salinityflux = .FALSE.
+
+    IF ( ibc_sa_t == 1  .AND.  top_salinityflux == 9999999.9_wp )  THEN
+       message_string = 'boundary condition: bc_sa_t = "' //                   &
+                        TRIM( bc_sa_t ) // '" requires to set top_salinityflux'
+       CALL message( 'ocean_check_parameters', 'PA0069', 1, 2, 0, 6, 0 )
+    ENDIF
+
+!
+!-- 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_wp )  THEN
+       message_string = 'boundary condition: bc_sa_t = "' //                   &
+                        TRIM( bc_sa_t ) // '" is not allowed with ' //         &
+                        'top_salinityflux /= 0.0'
+       CALL message( 'ocean_check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
+    ENDIF
+
+ END SUBROUTINE ocean_check_parameters
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check data output.
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_check_data_output( var, unit )
+ 
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  unit     !< unit of output variable
+    CHARACTER (LEN=*) ::  var      !< name of output variable
+
+
+    SELECT CASE ( TRIM( var ) )
+
+       CASE ( 'rho_ocean' )
+          unit = 'kg/m3'
+
+       CASE ( 'sa' )
+          unit = 'psu'
+
+       CASE DEFAULT
+          unit = 'illegal'
+
+    END SELECT
+
+ END SUBROUTINE ocean_check_data_output
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check data output of profiles
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_check_data_output_pr( variable, var_count, unit, dopr_unit )
+
+    USE arrays_3d,                                                             &
+        ONLY:  zu, zw
+
+    USE control_parameters,                                                    &
+        ONLY:  data_output_pr, message_string
+
+    USE indices
+
+    USE profil_parameter
+
+    USE statistics
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  unit      !<
+    CHARACTER (LEN=*) ::  variable  !<
+    CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
+
+    INTEGER(iwp) ::  var_count     !<
+
+    SELECT CASE ( TRIM( variable ) )
+
+       CASE ( 'prho' )
+          dopr_index(var_count) = 71
+          dopr_unit             = 'kg/m3'
+          hom(:,2,71,:) = SPREAD( zu, 2, statistic_regions+1 )
+          unit = dopr_unit
+
+       CASE ( 'rho_ocean' )
+          dopr_index(var_count) = 64
+          dopr_unit             = 'kg/m3'
+          hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
+          IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
+             dopr_initial_index(var_count) = 77
+             hom(:,2,77,:)             = SPREAD( zu, 2, statistic_regions+1 )
+             hom(nzb,2,77,:)           = 0.0_wp    ! because zu(nzb) is negative
+             data_output_pr(var_count) = data_output_pr(var_count)(2:)
+          ENDIF
+          unit = dopr_unit
+
+       CASE ( 'sa', '#sa' )
+          dopr_index(var_count) = 23
+          dopr_unit             = 'psu'
+          hom(:,2,23,:) = SPREAD( zu, 2, statistic_regions+1 )
+          IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
+             dopr_initial_index(var_count) = 26
+             hom(:,2,26,:)             = SPREAD( zu, 2, statistic_regions+1 )
+             hom(nzb,2,26,:)           = 0.0_wp    ! because zu(nzb) is negative
+             data_output_pr(var_count) = data_output_pr(var_count)(2:)
+          ENDIF
+          unit = dopr_unit
+
+       CASE ( 'w"sa"' )
+          dopr_index(var_count) = 65
+          dopr_unit             = 'psu m/s'
+          hom(:,2,65,:) = SPREAD( zw, 2, statistic_regions+1 )
+          unit = dopr_unit
+
+       CASE ( 'w*sa*' )
+          dopr_index(var_count) = 66
+          dopr_unit             = 'psu m/s'
+          hom(:,2,66,:) = SPREAD( zw, 2, statistic_regions+1 )
+          unit = dopr_unit
+
+       CASE ( 'wsa' )
+          dopr_index(var_count) = 67
+          dopr_unit             = 'psu m/s'
+          hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 )
+          unit = dopr_unit
+
+       CASE DEFAULT
+          unit = 'illegal'
+
+    END SELECT
+
+
+ END SUBROUTINE ocean_check_data_output_pr
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Define appropriate grid for netcdf variables.
+!> It is called out from subroutine netcdf.
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
+    
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !< x grid of output variable
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y   !< y grid of output variable
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z   !< z grid of output variable
+    CHARACTER (LEN=*), INTENT(IN)  ::  var      !< name of output variable
+
+    LOGICAL, INTENT(OUT) ::  found   !< flag if output variable is found
+
+    found  = .TRUE.
+
+!
+!-- Check for the grid
+    SELECT CASE ( TRIM( var ) )
+
+       CASE ( 'rho_ocean', 'sa' )
+          grid_x = 'x'
+          grid_y = 'y'
+          grid_z = 'zu'
+
+       CASE DEFAULT
+          found  = .FALSE.
+          grid_x = 'none'
+          grid_y = 'none'
+          grid_z = 'none'
+
+    END SELECT
+
+ END SUBROUTINE ocean_define_netcdf_grid
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Average 3D data.
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_3d_data_averaging( mode, variable )
+ 
+
+    USE arrays_3d,                                                             &
+        ONLY:  rho_ocean, sa
+
+    USE averaging,                                                             &
+        ONLY:  rho_ocean_av, sa_av
+
+    USE control_parameters,                                                    &
+        ONLY:  average_count_3d
+
+    USE indices,                                                               &
+        ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  mode       !< flag defining mode 'allocate', 'sum' or 'average'
+    CHARACTER (LEN=*) ::  variable   !< name of variable
+
+    INTEGER(iwp) ::  i   !< loop index
+    INTEGER(iwp) ::  j   !< loop index
+    INTEGER(iwp) ::  k   !< loop index
+
+    IF ( mode == 'allocate' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+
+          CASE ( 'rho_ocean' )
+             IF ( .NOT. ALLOCATED( rho_ocean_av ) )  THEN
+                ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             ENDIF
+             rho_ocean_av = 0.0_wp
+
+          CASE ( 'sa' )
+             IF ( .NOT. ALLOCATED( sa_av ) )  THEN
+                ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             ENDIF
+             sa_av = 0.0_wp
+
+          CASE DEFAULT
+             CONTINUE
+
+       END SELECT
+
+    ELSEIF ( mode == 'sum' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+
+          CASE ( 'rho_ocean' )
+             IF ( ALLOCATED( rho_ocean_av ) )  THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      DO  k = nzb, nzt+1
+                         rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) +           &
+                                               rho_ocean(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          CASE ( 'sa' )
+             IF ( ALLOCATED( sa_av ) )  THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      DO  k = nzb, nzt+1
+                         sa_av(k,j,i) = sa_av(k,j,i) + sa(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          CASE DEFAULT
+             CONTINUE
+
+       END SELECT
+
+    ELSEIF ( mode == 'average' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+
+          CASE ( 'rho_ocean' )
+             IF ( ALLOCATED( rho_ocean_av ) ) THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      DO  k = nzb, nzt+1
+                         rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) /           &
+                                               REAL( average_count_3d, KIND=wp )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          CASE ( 'sa' )
+             IF ( ALLOCATED( sa_av ) ) THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      DO  k = nzb, nzt+1
+                         sa_av(k,j,i) = sa_av(k,j,i) /                         &
+                                        REAL( average_count_3d, KIND=wp )
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+       END SELECT
+
+    ENDIF
+
+ END SUBROUTINE ocean_3d_data_averaging
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Define 2D output variables.
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_data_output_2d( av, variable, found, grid, mode, local_pf,   &
+                                  nzb_do, nzt_do )
+ 
+    USE arrays_3d,                                                             &
+        ONLY:  rho_ocean, sa
+
+    USE averaging,                                                             &
+        ONLY:  rho_ocean_av, sa_av
+
+    USE indices,                                                               &
+        ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  grid       !< name of vertical grid
+    CHARACTER (LEN=*) ::  mode       !< either 'xy', 'xz' or 'yz'
+    CHARACTER (LEN=*) ::  variable   !< name of variable
+
+    INTEGER(iwp) ::  av        !< flag for (non-)average output
+    INTEGER(iwp) ::  flag_nr   !< number of masking flag
+    INTEGER(iwp) ::  i         !< loop index
+    INTEGER(iwp) ::  j         !< loop index
+    INTEGER(iwp) ::  k         !< loop index
+    INTEGER(iwp) ::  nzb_do    !< vertical output index (bottom)
+    INTEGER(iwp) ::  nzt_do    !< vertical output index (top)
+
+    LOGICAL ::  found   !< flag if output variable is found
+    LOGICAL ::  resorted  !< flag if output is already resorted
+
+    REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
+
+    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< local
+       !< array to which output data is resorted to
+
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
+    
+    found = .TRUE.
+    resorted = .FALSE.
+!
+!-- Set masking flag for topography for not resorted arrays
+    flag_nr = 0
+
+    SELECT CASE ( TRIM( variable ) )
+
+       CASE ( 'rho_ocean_xy', 'rho_ocean_xz', 'rho_ocean_yz' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => rho_ocean
+          ELSE
+             IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN
+                ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                rho_ocean_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => rho_ocean_av
+          ENDIF
+
+       CASE ( 'sa_xy', 'sa_xz', 'sa_yz' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => sa
+          ELSE
+             IF ( .NOT. ALLOCATED( sa_av ) ) THEN
+                ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                sa_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => sa_av
+          ENDIF
+          IF ( mode == 'xy' ) grid = 'zu'
+
+       CASE DEFAULT
+          found = .FALSE.
+          grid  = 'none'
+
+    END SELECT
+
+    IF ( found .AND. .NOT. resorted )  THEN
+!       DO  i = nxl, nxr
+!          DO  j = nys, nyn
+!             DO  k = nzb_do, nzt_do
+!                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
+!                                         REAL( fill_value, KIND = wp ),        &
+!                                         BTEST( wall_flags_0(k,j,i), flag_nr ) )
+!             ENDDO
+!          ENDDO
+!       ENDDO
+    ENDIF
+ 
+ END SUBROUTINE ocean_data_output_2d
+
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Define 3D output variables.
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
+ 
+
+    USE arrays_3d,                                                             &
+        ONLY:  rho_ocean, sa
+
+    USE averaging,                                                             &
+        ONLY:  rho_ocean_av, sa_av
+
+    USE indices,                                                               &
+        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,           &
+               wall_flags_0
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  variable   !< name of variable
+
+    INTEGER(iwp) ::  av        !< flag for (non-)average output
+    INTEGER(iwp) ::  flag_nr   !< number of masking flag
+    INTEGER(iwp) ::  i         !< loop index
+    INTEGER(iwp) ::  j         !< loop index
+    INTEGER(iwp) ::  k         !< loop index
+    INTEGER(iwp) ::  nzb_do    !< lower limit of the data output (usually 0)
+    INTEGER(iwp) ::  nzt_do    !< vertical upper limit of the data output (usually nz_do3d)
+
+    LOGICAL ::  found     !< flag if output variable is found
+    LOGICAL ::  resorted  !< flag if output is already resorted
+
+    REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
+
+    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< local
+                                  !< array to which output data is resorted to
+
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
+
+    found = .TRUE.
+    resorted = .FALSE.
+!
+!-- Set masking flag for topography for not resorted arrays
+    flag_nr = 0
+
+    SELECT CASE ( TRIM( variable ) )
+
+       CASE ( 'rho_ocean' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => rho_ocean
+          ELSE
+             IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN
+                ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                rho_ocean_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => rho_ocean_av
+          ENDIF
+
+       CASE ( 'sa' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => sa
+          ELSE
+             IF ( .NOT. ALLOCATED( sa_av ) ) THEN
+                ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                sa_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => sa_av
+          ENDIF
+
+       CASE DEFAULT
+          found = .FALSE.
+
+    END SELECT
+
+
+    IF ( found  .AND.  .NOT. resorted )  THEN
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_do, nzt_do
+                local_pf(i,j,k) = MERGE(                                       &
+                                   to_be_resorted(k,j,i),                      &
+                                   REAL( fill_value, KIND = wp ),              &
+                                   BTEST( wall_flags_0(k,j,i), flag_nr ) )
+             ENDDO
+          ENDDO
+       ENDDO
+       resorted = .TRUE.
+    ENDIF
+
+ END SUBROUTINE ocean_data_output_3d
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Allocate arrays and assign pointers.
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_init_arrays
+
+    USE indices,                                                               &
+        ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt
+
+    USE pmc_interface,                                                         &
+        ONLY:  nested_run
+
+    IMPLICIT NONE
+
+#if defined( __nopointer )
+    ALLOCATE( prho(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
+              rho_ocean(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
+              sa(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                            &
+              sa_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
+              tsa_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+#else
+    ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                        &
+              rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                         &
+              sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
+              sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
+              sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
+    prho => prho_1
+    rho_ocean  => rho_1  ! routines calc_mean_profile and diffusion_e require
+                         ! density to be a pointer
+#endif
+
+#if ! defined( __nopointer )
+!
+!-- Initial assignment of pointers
+    sa  => sa_1;   sa_p  => sa_2;   tsa_m  => sa_3
+#endif
+
+ END SUBROUTINE ocean_init_arrays
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialization of quantities needed for the ocean mode
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_init
+
+
+    USE arrays_3d,                                                             &
+        ONLY:  dzu, hyp, pt_init, ref_state, zu, zw
+
+    USE basic_constants_and_equations_mod,                                     &
+        ONLY:  g
+
+    USE control_parameters,                                                    &
+        ONLY:  initializing_actions, molecular_viscosity, rho_surface,         &
+               rho_reference, surface_pressure, use_single_reference_value
+
+    USE indices,                                                               &
+        ONLY:  nxl, nxlg, nxrg, nyng, nys, nysg, nzb, nzt
+
+    USE kinds
+
+    USE pegrid
+
+    USE statistics,                                                            &
+        ONLY:  hom, statistic_regions
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i  !< loop index
+    INTEGER(iwp) ::  j  !< loop index
+    INTEGER(iwp) ::  k  !< loop index
+    INTEGER(iwp) ::  n  !< loop index
+
+    REAL(wp)     ::  dum   !< dummy argument
+    REAL(wp)     ::  pt_l  !< local scalar for pt used in equation of state function
+    REAL(wp)     ::  sa_l  !< local scalar for sa used in equation of state function
+
+    REAL(wp), DIMENSION(nzb:nzt+1) ::  rho_ocean_init  !< local array for initial density
+
+    ALLOCATE( hyp(nzb:nzt+1) )
+
+
+!
+!-- In case of no restart run, calculate the inital salinity profilevcusing the
+!-- given salinity gradients
+    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
+
+       sa_init = sa_surface
+!
+!--    Last arguments gives back the gradient at top level to be used as
+!--    possible Neumann boundary condition. This is not realized for the ocean
+!--    mode, therefore a dummy argument is used.
+       CALL init_vertical_profiles( sa_vertical_gradient_level_ind,          &
+                                    sa_vertical_gradient_level,              &
+                                    sa_vertical_gradient, sa_init,           &
+                                    sa_surface, dum )
+    ENDIF
+
+!
+!-- Initialize required 3d-arrays
+    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
+         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
+!
+!--    Initialization via computed 1D-model profiles
+       IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 )  THEN
+
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+                sa(:,j,i) = sa_init
+             ENDDO
+          ENDDO
+
+       ENDIF
+!
+!--    Store initial profiles for output purposes etc.
+!--    Store initial salinity profile
+       hom(:,1,26,:)  = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 )
+!
+!--    Initialize old and new time levels.
+       tsa_m = 0.0_wp
+       sa_p  = sa
+
+    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
+
+!
+!--    Initialize new time levels (only done in order to set boundary values
+!--    including ghost points)
+       sa_p = sa
+!
+!--    Allthough tendency arrays are set in prognostic_equations, they have
+!--    have to be predefined here because they are used (but multiplied with 0)
+!--    there before they are set.
+       tsa_m = 0.0_wp
+
+    ENDIF
+
+!
+!-- Set water density near the ocean surface
+    rho_surface = 1027.62_wp
+
+!
+!-- Set kinematic viscosity to sea water at 20C.
+!-- This changes the default value that is given for air!
+    molecular_viscosity = 1.05E-6_wp
+
+!
+!-- Change sign of buoyancy/stability terms because density gradient is used
+!-- instead of the potential temperature gradient to calculate the buoyancy
+    atmos_ocean_sign = -1.0_wp
+
+!
+!-- Calculate initial vertical profile of hydrostatic pressure (in Pa)
+!-- and the reference density (used later in buoyancy term)
+!-- First step: Calculate pressure using reference density
+    hyp(nzt+1) = surface_pressure * 100.0_wp
+    hyp(nzt)   = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1)
+    rho_ocean_init(nzt)   = rho_surface
+    rho_ocean_init(nzt+1) = rho_surface  ! only required for output
+
+    DO  k = nzt-1, 1, -1
+       hyp(k) = hyp(k+1) + rho_surface * g * dzu(k)
+    ENDDO
+    hyp(0) = hyp(1) + rho_surface * g * dzu(1)
+
+!
+!-- Second step: Iteratively calculate in situ density (based on presssure)
+!-- and pressure (based on in situ density)
+    DO  n = 1, 5
+
+       rho_reference = rho_surface * 0.5_wp * dzu(nzt+1)
+
+       DO  k = nzt, 0, -1
+
+          sa_l = 0.5_wp * ( sa_init(k) + sa_init(k+1) )
+          pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) )
+
+          rho_ocean_init(k) = eqn_state_seawater_func( hyp(k), pt_l, sa_l )
+
+          rho_reference = rho_reference + rho_ocean_init(k) * dzu(k+1)
+
+       ENDDO
+
+       rho_reference = rho_reference / ( zw(nzt) - zu(nzb) )
+
+       hyp(nzt) = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1)
+       DO  k = nzt-1, 0, -1
+          hyp(k) = hyp(k+1) + g * 0.5_wp * ( rho_ocean_init(k)                 &
+                                           + rho_ocean_init(k+1) ) * dzu(k+1)
+       ENDDO
+
+    ENDDO
+
+!
+!-- Calculate the reference potential density
+    prho_reference = 0.0_wp
+    DO  k = 0, nzt
+
+       sa_l = 0.5_wp * ( sa_init(k) + sa_init(k+1) )
+       pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) )
+
+       prho_reference = prho_reference + dzu(k+1) * &
+                        eqn_state_seawater_func( 0.0_wp, pt_l, sa_l )
+
+    ENDDO
+
+    prho_reference = prho_reference / ( zu(nzt) - zu(nzb) )
+
+!
+!-- Calculate the 3d array of initial in situ and potential density,
+!-- based on the initial temperature and salinity profile
+    CALL eqn_state_seawater
+
+!
+!-- Store initial density profile
+    hom(:,1,77,:)  = SPREAD( rho_ocean_init(:), 2, statistic_regions+1 )
+
+!
+!-- Set the reference state to be used in the buoyancy terms
+    IF ( use_single_reference_value )  THEN
+       ref_state(:) = rho_reference
+    ELSE
+       ref_state(:) = rho_ocean_init(:)
+    ENDIF
+
+
+ END SUBROUTINE ocean_init
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Prognostic equation for salinity.
+!> Vector-optimized version
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_prognostic_equations
+
+    USE advec_s_bc_mod,                                                        &
+        ONLY:  advec_s_bc
+
+    USE advec_s_pw_mod,                                                        &
+        ONLY:  advec_s_pw
+
+    USE advec_s_up_mod,                                                        &
+        ONLY:  advec_s_up
+
+    USE advec_ws,                                                              &
+        ONLY:  advec_s_ws
+
+    USE arrays_3d,                                                             &
+        ONLY:  rdf_sc, tend, tsa_m
+
+    USE control_parameters,                                                    &
+        ONLY:  dt_3d, intermediate_timestep_count, intermediate_timestep_count_max,   &
+               scalar_advec, timestep_scheme, tsc, ws_scheme_sca
+
+    USE cpulog,                                                                &
+        ONLY:  cpu_log, log_point
+
+    USE diffusion_s_mod,                                                       &
+        ONLY:  diffusion_s
+
+    USE indices,                                                               &
+        ONLY:  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
+
+    USE surface_mod,                                                           &
+        ONLY:  surf_def_v, surf_def_h, surf_lsm_h, surf_lsm_v, surf_usm_h,     &
+               surf_usm_v
+
+    USE user_actions_mod,                                                      &
+        ONLY:  user_actions
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i       !< loop index
+    INTEGER(iwp) ::  j       !< loop index
+    INTEGER(iwp) ::  k       !< loop index
+
+    REAL(wp)     ::  sbt     !< weighting factor for sub-time step
+
+!
+!-- Compute prognostic equations for the ocean mode
+!-- First, start with salinity
+    CALL cpu_log( log_point(37), 'sa-equation', 'start' )
+
+!
+!-- sa-tendency terms with communication
+    sbt = tsc(2)
+    IF ( scalar_advec == 'bc-scheme' )  THEN
+
+       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--       Bott-Chlond scheme always uses Euler time step. Thus:
+          sbt = 1.0_wp
+       ENDIF
+       tend = 0.0_wp
+       CALL advec_s_bc( sa, 'sa' )
+
+    ENDIF
+
+!
+!-- sa-tendency terms with no communication
+    IF ( scalar_advec /= 'bc-scheme' )  THEN
+       tend = 0.0_wp
+       IF ( timestep_scheme(1:5) == 'runge' )  THEN
+          IF ( ws_scheme_sca )  THEN
+             CALL advec_s_ws( sa, 'sa' )
+          ELSE
+             CALL advec_s_pw( sa )
+          ENDIF
+       ELSE
+          CALL advec_s_up( sa )
+       ENDIF
+    ENDIF
+
+    CALL diffusion_s( sa,                                                      &
+                      surf_def_h(0)%sasws, surf_def_h(1)%sasws,                &
+                      surf_def_h(2)%sasws,                                     &
+                      surf_lsm_h%sasws,    surf_usm_h%sasws,                   &
+                      surf_def_v(0)%sasws, surf_def_v(1)%sasws,                &
+                      surf_def_v(2)%sasws, surf_def_v(3)%sasws,                &
+                      surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws,                &
+                      surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws,                &
+                      surf_usm_v(0)%sasws, surf_usm_v(1)%sasws,                &
+                      surf_usm_v(2)%sasws, surf_usm_v(3)%sasws )
+
+    CALL user_actions( 'sa-tendency' )
+
+!
+!-- Prognostic equation for salinity
+    DO  i = nxl, nxr
+       DO  j = nys, nyn
+          DO  k = nzb+1, nzt
+             sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +         &
+                                                   tsc(3) * tsa_m(k,j,i) )     &
+                                               - tsc(5) * rdf_sc(k) *          &
+                                                 ( sa(k,j,i) - sa_init(k) )    &
+                                       )                                       &
+                                         * MERGE( 1.0_wp, 0.0_wp,              &
+                                               BTEST( wall_flags_0(k,j,i), 0 ) &
+                                                )
+             IF ( sa_p(k,j,i) < 0.0_wp )  sa_p(k,j,i) = 0.1_wp * 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+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+1, nzt
+                   tsa_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
+                                     5.3125_wp * tsa_m(k,j,i)
+                ENDDO
+             ENDDO
+         ENDDO
+       ENDIF
+    ENDIF
+
+    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' )
+
+ END SUBROUTINE ocean_prognostic_equations
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Prognostic equations for ocean mode (so far, salinity only)
+!> Cache-optimized version
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_prognostic_equations_ij( i, j, i_omp_start, tn )
+
+    USE advec_s_pw_mod,                                                        &
+        ONLY:  advec_s_pw
+
+    USE advec_s_up_mod,                                                        &
+        ONLY:  advec_s_up
+
+    USE advec_ws,                                                              &
+        ONLY:  advec_s_ws
+
+    USE arrays_3d,                                                             &
+        ONLY:  diss_l_sa, diss_s_sa, flux_l_sa, flux_s_sa, rdf_sc, tend, tsa_m
+
+    USE control_parameters,                                                    &
+        ONLY:  dt_3d, intermediate_timestep_count,                             &
+               intermediate_timestep_count_max, timestep_scheme, tsc,          &
+               ws_scheme_sca
+
+    USE diffusion_s_mod,                                                       &
+        ONLY:  diffusion_s
+
+    USE indices,                                                               &
+        ONLY:  nzb, nzt, wall_flags_0
+
+    USE surface_mod,                                                           &
+        ONLY:  surf_def_v, surf_def_h, surf_lsm_h, surf_lsm_v, surf_usm_h,     &
+               surf_usm_v
+
+    USE user_actions_mod,                                                      &
+        ONLY:  user_actions
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i             !< loop index x direction
+    INTEGER(iwp) ::  i_omp_start   !< first loop index of i-loop in calling    &
+                                   !< routine prognostic_equations
+    INTEGER(iwp) ::  j             !< loop index y direction
+    INTEGER(iwp) ::  k             !< loop index z direction
+    INTEGER(iwp) ::  tn            !< task number of openmp task
+
+!
+!-- Compute prognostic equations for the ocean mode
+!-- First, start with tendency-terms for salinity
+    tend(:,j,i) = 0.0_wp
+    IF ( timestep_scheme(1:5) == 'runge' ) &
+    THEN
+       IF ( ws_scheme_sca )  THEN
+          CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa,  diss_s_sa, flux_l_sa,   &
+                           diss_l_sa, i_omp_start, tn  )
+       ELSE
+          CALL advec_s_pw( i, j, sa )
+       ENDIF
+    ELSE
+       CALL advec_s_up( i, j, sa )
+    ENDIF
+    CALL diffusion_s( i, j, sa,                                       &
+                      surf_def_h(0)%sasws, surf_def_h(1)%sasws,       &
+                      surf_def_h(2)%sasws,                            &
+                      surf_lsm_h%sasws,    surf_usm_h%sasws,          &
+                      surf_def_v(0)%sasws, surf_def_v(1)%sasws,       &
+                      surf_def_v(2)%sasws, surf_def_v(3)%sasws,       &
+                      surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws,       &
+                      surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws,       &
+                      surf_usm_v(0)%sasws, surf_usm_v(1)%sasws,       &
+                      surf_usm_v(2)%sasws, surf_usm_v(3)%sasws )
+
+    CALL user_actions( i, j, 'sa-tendency' )
+
+!
+!-- Prognostic equation for salinity
+    DO  k = nzb+1, nzt
+
+       sa_p(k,j,i) = sa(k,j,i) + ( dt_3d *                                     &
+                                           ( tsc(2) * tend(k,j,i) +            &
+                                             tsc(3) * tsa_m(k,j,i) )           &
+                                 - tsc(5) * rdf_sc(k)                          &
+                                          * ( sa(k,j,i) - sa_init(k) )         &
+                                 ) * MERGE( 1.0_wp, 0.0_wp,                    &
+                                            BTEST( wall_flags_0(k,j,i), 0 ) )
+
+       IF ( sa_p(k,j,i) < 0.0_wp )  sa_p(k,j,i) = 0.1_wp * 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+1, nzt
+             tsa_m(k,j,i) = tend(k,j,i)
+          ENDDO
+       ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )&
+       THEN
+          DO  k = nzb+1, nzt
+             tsa_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
+                               5.3125_wp * tsa_m(k,j,i)
+          ENDDO
+       ENDIF
+    ENDIF
+
+!
+!-- Calculate density by the equation of state for seawater
+    CALL eqn_state_seawater( i, j )
+
+ END SUBROUTINE ocean_prognostic_equations_ij
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Swapping of timelevels.
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_swap_timelevel( mod_count )
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(IN) ::  mod_count  !< flag defining where pointers point to
+
+#if defined( __nopointer )
+
+    sa = sa_p
+
+#else
+    
+    SELECT CASE ( mod_count )
+
+       CASE ( 0 )
+          sa => sa_1;    sa_p => sa_2
+
+       CASE ( 1 )
+          sa => sa_2;    sa_p => sa_1
+
+    END SELECT
+
+#endif
+
+ END SUBROUTINE ocean_swap_timelevel
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine reads the respective restart data for the ocean module.
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_rrd_global( found )
+
+
+    USE control_parameters,                                                    &
+        ONLY: length, restart_string
+
+
+    IMPLICIT NONE
+
+    LOGICAL, INTENT(OUT)  ::  found
+
+
+    found = .TRUE.
+
+    SELECT CASE ( restart_string(1:length) )
+
+       CASE ( 'bc_sa_t' )
+          READ ( 13 )  bc_sa_t
+
+       CASE ( 'bottom_salinityflux' )
+          READ ( 13 )  bottom_salinityflux
+
+       CASE ( 'langmuir_number' )
+          READ ( 13 )  langmuir_number
+
+       CASE ( 'sa_init' )
+          READ ( 13 )  sa_init
+
+       CASE ( 'sa_surface' )
+          READ ( 13 )  sa_surface
+
+       CASE ( 'sa_vertical_gradient' )
+          READ ( 13 )  sa_vertical_gradient
+
+       CASE ( 'sa_vertical_gradient_level' )
+          READ ( 13 )  sa_vertical_gradient_level
+
+       CASE ( 'stokes_force' )
+          READ ( 13 )  stokes_force
+
+       CASE ( 'stokes_waveheight' )
+          READ ( 13 )  stokes_waveheight
+
+       CASE ( 'stokes_wavelength' )
+          READ ( 13 )  stokes_wavelength
+
+       CASE ( 'top_salinityflux' )
+          READ ( 13 )  top_salinityflux
+
+       CASE ( 'wall_salinityflux' )
+          READ ( 13 )  wall_salinityflux
+
+       CASE DEFAULT
+
+          found = .FALSE.
+
+    END SELECT
+
+ END SUBROUTINE ocean_rrd_global
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine reads the respective restart data for the ocean module.
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,        &
+                             nxr_on_file, nynf, nync, nyn_on_file, nysf,       &
+                             nysc, nys_on_file, tmp_2d, tmp_3d, found )
+
+    USE averaging,                                                             &
+        ONLY:  rho_ocean_av, sa_av
+
+    USE control_parameters,                                                    &
+        ONLY:  length, restart_string
+
+    USE indices,                                                               &
+        ONLY:  nbgp, nxlg, nxrg, nyng, nysg, nzb, nzt
+
+    USE pegrid
+
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i               !<
+    INTEGER(iwp) ::  k               !<
+    INTEGER(iwp) ::  nxlc            !<
+    INTEGER(iwp) ::  nxlf            !<
+    INTEGER(iwp) ::  nxl_on_file     !<
+    INTEGER(iwp) ::  nxrc            !<
+    INTEGER(iwp) ::  nxrf            !<
+    INTEGER(iwp) ::  nxr_on_file     !<
+    INTEGER(iwp) ::  nync            !<
+    INTEGER(iwp) ::  nynf            !<
+    INTEGER(iwp) ::  nyn_on_file     !<
+    INTEGER(iwp) ::  nysc            !<
+    INTEGER(iwp) ::  nysf            !<
+    INTEGER(iwp) ::  nys_on_file     !<
+
+    LOGICAL, INTENT(OUT)  ::  found
+
+    REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d   !<
+    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
+
+
+    found = .TRUE.
+
+    SELECT CASE ( restart_string(1:length) )
+
+       CASE ( 'rho_ocean_av' )
+          IF ( .NOT. ALLOCATED( rho_ocean_av ) )  THEN
+             ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( k == 1 )  READ ( 13 )  tmp_3d
+          rho_ocean_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =            &
+                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+       CASE ( 'sa' )
+          IF ( k == 1 )  READ ( 13 )  tmp_3d
+          sa(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                      &
+                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+       CASE ( 'sa_av' )
+          IF ( .NOT. ALLOCATED( sa_av ) )  THEN
+              ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( k == 1 )  READ ( 13 )  tmp_3d
+          sa_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                   &
+                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+       CASE DEFAULT
+          found = .FALSE.
+
+    END SELECT
+
+ END SUBROUTINE ocean_rrd_local
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine writes the respective restart data for the ocean module.
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_wrd_global
+
+
+    IMPLICIT NONE
+
+    CALL wrd_write_string( 'bc_sa_t' )
+    WRITE ( 14 )  bc_sa_t
+
+    CALL wrd_write_string( 'bottom_salinityflux' )
+    WRITE ( 14 )  bottom_salinityflux
+
+    CALL wrd_write_string( 'langmuir_number' )
+    WRITE ( 14 )  langmuir_number
+
+    CALL wrd_write_string( 'sa_init' )
+    WRITE ( 14 )  sa_init
+
+    CALL wrd_write_string( 'sa_surface' )
+    WRITE ( 14 )  sa_surface
+
+    CALL wrd_write_string( 'sa_vertical_gradient' )
+    WRITE ( 14 )  sa_vertical_gradient
+
+    CALL wrd_write_string( 'sa_vertical_gradient_level' )
+    WRITE ( 14 )  sa_vertical_gradient_level
+
+    CALL wrd_write_string( 'stokes_force' )
+    WRITE ( 14 )  stokes_force
+
+    CALL wrd_write_string( 'stokes_waveheight' )
+    WRITE ( 14 )  stokes_waveheight
+
+    CALL wrd_write_string( 'stokes_wavelength' )
+    WRITE ( 14 )  stokes_wavelength
+
+    CALL wrd_write_string( 'top_salinityflux' )
+    WRITE ( 14 )  top_salinityflux
+
+    CALL wrd_write_string( 'wall_salinityflux' )
+    WRITE ( 14 )  wall_salinityflux
+
+ END SUBROUTINE ocean_wrd_global
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine writes the respective restart data for the ocean module.
+!------------------------------------------------------------------------------!
+ SUBROUTINE ocean_wrd_local
+
+    USE averaging,                                                             &
+        ONLY:  rho_ocean_av, sa_av
+
+    IMPLICIT NONE
+
+    IF ( ALLOCATED( rho_ocean_av ) )  THEN
+       CALL wrd_write_string( 'rho_ocean_av' )
+       WRITE ( 14 )  rho_ocean_av
+    ENDIF
+
+    CALL wrd_write_string( 'sa' )
+    WRITE ( 14 )  sa
+
+    IF ( ALLOCATED( sa_av ) )  THEN
+       CALL wrd_write_string( 'sa_av' )
+       WRITE ( 14 )  sa_av
+    ENDIF
+
+ END SUBROUTINE ocean_wrd_local
+
+
+ END MODULE ocean_mod
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 3291)
+++ palm/trunk/SOURCE/parin.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! changes concerning modularization of ocean option
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -471,4 +474,7 @@
         ONLY:  netcdf_data_format, netcdf_deflate, netcdf_precision
 
+    USE ocean_mod,                                                             &
+        ONLY:  ocean_parin
+
     USE pegrid
                
@@ -527,6 +533,6 @@
     NAMELIST /inipar/  alpha_surface, approximation, 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_sa_t, bc_uv_b, bc_uv_t,                 &
-             bottom_salinityflux, building_height, building_length_x,          &
+             bc_q_t,bc_s_b, bc_s_t, bc_uv_b, bc_uv_t,                 &
+             building_height, building_length_x,          &
              building_length_y, building_wall_left, building_wall_south,       &
              calc_soil_moisture_during_spinup,                                 &
@@ -559,6 +565,6 @@
              most_method, nesting_offline,                  &
              netcdf_precision, neutral, ngsrb,                                 &
-             nsor, nsor_ini, nudging, nx, ny, nz, ocean, omega, omega_sor,     &
-             outflow_source_plane, passive_scalar,                             &
+             nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega,           &
+             omega_sor, outflow_source_plane, passive_scalar,                  &
              prandtl_number, psolver, pt_damping_factor,        &
              pt_damping_width, pt_reference, pt_surface,                       &
@@ -571,6 +577,6 @@
              recycling_width, recycling_yshift,                                &
              reference_state, residual_limit,                                  &
-             roughness_length, sa_surface,                                     &
-             sa_vertical_gradient, sa_vertical_gradient_level, scalar_advec,   &
+             roughness_length,                                                 &
+             scalar_advec,   &
              scalar_rayleigh_damping,                              &
              spinup_time, spinup_pt_amplitude, spinup_pt_mean,                 &
@@ -581,5 +587,5 @@
              s_vertical_gradient_level, time_utc_init, timestep_scheme,        &
              topography, topography_grid_convention, top_heatflux,             &
-             top_momentumflux_u, top_momentumflux_v, top_salinityflux,         &
+             top_momentumflux_u, top_momentumflux_v,                           &
              top_scalarflux, transpose_compute_overlap,                        &
              tunnel_height, tunnel_length, tunnel_width_x, tunnel_width_y,     &
@@ -593,28 +599,28 @@
              vg_vertical_gradient_level, v_bulk, v_profile,&
              wall_adjustment, wall_heatflux, wall_humidityflux,                &
-             wall_salinityflux, wall_scalarflux, y_shift, zeta_max, zeta_min,  &
+             wall_scalarflux, y_shift, zeta_max, zeta_min,  &
              z0h_factor
 
-    NAMELIST /initialization_parameters/  alpha_surface,         &
+    NAMELIST /initialization_parameters/  alpha_surface,                       &
              approximation, 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_sa_t, bc_uv_b, bc_uv_t,                 &
-             bottom_salinityflux, building_height, building_length_x,          &
+             bc_q_t,bc_s_b, bc_s_t, bc_uv_b, bc_uv_t,                          &
+             building_height, building_length_x,                               &
              building_length_y, building_wall_left, building_wall_south,       &
              calc_soil_moisture_during_spinup,                                 &
-             call_psolver_at_all_substeps,  &
+             call_psolver_at_all_substeps,                                     &
              canyon_height,                                                    &
              canyon_width_x, canyon_width_y, canyon_wall_left,                 &
-             canyon_wall_south, cfl_factor, cloud_droplets,   &
-             cloud_top_radiation,                 &
-             collective_wait, complex_terrain,           &
+             canyon_wall_south, cfl_factor, cloud_droplets,                    &
+             cloud_top_radiation,                                              &
+             collective_wait, complex_terrain,                                 &
              conserve_volume_flow,                                             &
              conserve_volume_flow_mode, constant_flux_layer,                   &
-             coupling_start_time,             &
+             coupling_start_time,                                              &
              cycle_mg, damp_level_1d,                                          &
              data_output_during_spinup,                                        &
              day_of_year_init,                                                 &
              dissipation_1d,                                                   &
-             dp_external, dp_level_b, dp_smooth, dpdxy,    &
+             dp_external, dp_level_b, dp_smooth, dpdxy,                        &
              dt, dt_pr_1d, dt_run_control_1d, dt_spinup, dx, dy, dz, dz_max,   &
              dz_stretch_factor, dz_stretch_level, dz_stretch_level_start,      &
@@ -626,12 +632,12 @@
              initializing_actions, km_constant,                                &
              large_scale_forcing, large_scale_subsidence, latitude,            &
-             longitude,                                 &
+             longitude,                                                        &
              loop_optimization, lsf_exception, masking_method, mg_cycles,      &
              mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
-             most_method, nesting_offline,                  &
+             most_method, nesting_offline,                                     &
              netcdf_precision, neutral, ngsrb,                                 &
-             nsor, nsor_ini, nudging, nx, ny, nz, ocean, omega, omega_sor,     &
-             outflow_source_plane, passive_scalar,                             &
-             prandtl_number, psolver, pt_damping_factor,        &
+             nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega,           &
+             omega_sor, outflow_source_plane, passive_scalar,                  &
+             prandtl_number, psolver, pt_damping_factor,                       &
              pt_damping_width, pt_reference, pt_surface,                       &
              pt_surface_initial_change, pt_vertical_gradient,                  &
@@ -643,7 +649,6 @@
              recycling_width, recycling_yshift,                                &
              reference_state, residual_limit,                                  &
-             roughness_length, sa_surface,                                     &
-             sa_vertical_gradient, sa_vertical_gradient_level, scalar_advec,   &
-             scalar_rayleigh_damping,                              &
+             roughness_length, scalar_advec,                                   &
+             scalar_rayleigh_damping,                                          &
              spinup_time, spinup_pt_amplitude, spinup_pt_mean,                 &
              statistic_regions, subs_vertical_gradient,                        &
@@ -653,5 +658,5 @@
              s_vertical_gradient_level, time_utc_init, timestep_scheme,        &
              topography, topography_grid_convention, top_heatflux,             &
-             top_momentumflux_u, top_momentumflux_v, top_salinityflux,         &
+             top_momentumflux_u, top_momentumflux_v,                           &
              top_scalarflux, transpose_compute_overlap,                        &
              tunnel_height, tunnel_length, tunnel_width_x, tunnel_width_y,     &
@@ -663,8 +668,7 @@
              use_free_convection_scaling,                                      &
              uv_heights, u_bulk, u_profile, vg_surface, vg_vertical_gradient,  &
-             vg_vertical_gradient_level, v_bulk, v_profile,&
+             vg_vertical_gradient_level, v_bulk, v_profile,                    &
              wall_adjustment, wall_heatflux, wall_humidityflux,                &
-             wall_salinityflux, wall_scalarflux, y_shift, zeta_max, zeta_min,  &
-             z0h_factor
+             wall_scalarflux, y_shift, zeta_max, zeta_min, z0h_factor
              
     NAMELIST /d3par/  averaging_interval, averaging_interval_pr,               &
@@ -858,56 +862,19 @@
 
 !
-!--       Check if land surface model is used and read &lsm_par if required
+!--       Check for module namelists and read them
           CALL lsm_parin
-
-!
-!--       Check if gust module is used and read &gust_par if required
           CALL bcm_parin
-
-!
-!--       Check if urban surface model is used and read &urban_surface_par if required
           CALL usm_parin
-
-!
-!--       Check if spectra shall be calculated and read spectra_par if required
           CALL spectra_parin
-
-!
-!--       Check if radiation model is used and read &radiation_par if required
           CALL radiation_parin
-
-!
-!--       Check if gust module is used and read &gust_par if required
           CALL gust_parin
-
-!
-!--       Check if multi agent system is used and read &agents_par if required
           CALL mas_parin
-
-!
-!--       Check if plant canopy model is used and read &canopy_par if required
+          CALL ocean_parin
           CALL pcm_parin
-
-!
-!--       Read control parameters for optionally used model software packages
           CALL package_parin
-
-!
-!--       Check if wind turbine model is used and read &wind_turbine_par if
-!--       required
           CALL wtm_parin
-!
-!--       Check if virtual flights should be carried out and read &flight_par
-!--       if required
           CALL flight_parin
-!
-!--       Check if synthetic turbulence generator is used and read stg_par if
-!--       required
           CALL stg_parin
-!
-!--       Read chemistry variables
           CALL chem_parin
-!
-!--       Check if uv exposure model is used and read &uvexposure_par
           CALL uvem_parin
 !
Index: palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 3291)
+++ palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! ocean renamed ocean_mode
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -685,5 +688,5 @@
 
        USE control_parameters,                                                 &
-           ONLY: humidity, message_string, ocean, urban_surface
+           ONLY: humidity, message_string, ocean_mode, urban_surface
 
        USE netcdf_data_input_mod,                                              &
@@ -732,5 +735,5 @@
           gradient = 0.0_wp
 
-          IF (  .NOT.  ocean )  THEN
+          IF (  .NOT.  ocean_mode )  THEN
 
              lad(0) = lad_surface
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 3291)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! changes concerning modularization of ocean option
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -277,4 +280,33 @@
 
 
+    USE advec_s_bc_mod,                                                        &
+        ONLY:  advec_s_bc
+
+    USE advec_s_pw_mod,                                                        &
+        ONLY:  advec_s_pw
+
+    USE advec_s_up_mod,                                                        &
+        ONLY:  advec_s_up
+
+    USE advec_u_pw_mod,                                                        &
+        ONLY:  advec_u_pw
+
+    USE advec_u_up_mod,                                                        &
+        ONLY:  advec_u_up
+
+    USE advec_v_pw_mod,                                                        &
+        ONLY:  advec_v_pw
+
+    USE advec_v_up_mod,                                                        &
+        ONLY:  advec_v_up
+
+    USE advec_w_pw_mod,                                                        &
+        ONLY:  advec_w_pw
+
+    USE advec_w_up_mod,                                                        &
+        ONLY:  advec_w_up
+
+    USE advec_ws,                                                              &
+        ONLY:  advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws
 
     USE arrays_3d,                                                             &
@@ -287,8 +319,7 @@
                flux_l_sa, nc, nc_p, nr, nr_p, pt, ptdf_x, ptdf_y, pt_init,     &
                pt_p, prho, q, q_init, q_p, qc, qc_p, qr, qr_p, rdf, rdf_sc,    &
-               ref_state, rho_ocean, s,  s_init, s_p, sa, sa_init, sa_p, tend, &
-               te_m, tnc_m,  tnr_m, tpt_m, tq_m, tqc_m, tqr_m, ts_m, tsa_m,    &
-               tu_m, tv_m, tw_m, u, ug, u_init, u_p, v, vg, vpt, v_init, v_p,  &
-               w, w_p
+               ref_state, rho_ocean, s, s_init, s_p, tend, te_m, tnc_m,        &
+               tnr_m, tpt_m, tq_m, tqc_m, tqr_m, ts_m, tu_m, tv_m, tw_m, u,    &
+               ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
 
     USE bulk_cloud_model_mod,                                                  &
@@ -297,13 +328,19 @@
                microphysics_morrison, microphysics_seifert
 
+    USE buoyancy_mod,                                                          &
+        ONLY:  buoyancy
+
+    USE calc_radiation_mod,                                                    &
+        ONLY:  calc_radiation
+
     USE chemistry_model_mod,                                                   &
         ONLY:  chem_integrate, chem_prognostic_equations,                      &
                chem_species, nspec, nvar, spc_names
            
+    USE chem_modules,                                                          &
+        ONLY:  call_chem_at_all_substeps, chem_gasphase_on
+
     USE chem_photolysis_mod,                                                   &
         ONLY:  photolysis_control
-
-    USE chem_modules,                                                          &
-        ONLY:  call_chem_at_all_substeps, chem_gasphase_on
 
     USE control_parameters,                                                    &
@@ -315,7 +352,5 @@
                large_scale_subsidence,                  &
                neutral, nudging, &
-               ocean, passive_scalar, plant_canopy,                            &
-               prho_reference, prho_reference,                                 &
-               prho_reference, pt_reference, pt_reference, pt_reference,       &
+               ocean_mode, passive_scalar, plant_canopy, pt_reference,         &
                scalar_advec, scalar_advec, simulated_time, sloping_surface,    &
                timestep_scheme, tsc, use_subsidence_tendencies,                &
@@ -323,10 +358,21 @@
                ws_scheme_sca, urban_surface, land_surface
 
+    USE coriolis_mod,                                                          &
+        ONLY:  coriolis
+
     USE cpulog,                                                                &
         ONLY:  cpu_log, log_point, log_point_s
 
-
-    USE eqn_state_seawater_mod,                                                &
-        ONLY:  eqn_state_seawater
+    USE diffusion_s_mod,                                                       &
+        ONLY:  diffusion_s
+
+    USE diffusion_u_mod,                                                       &
+        ONLY:  diffusion_u
+
+    USE diffusion_v_mod,                                                       &
+        ONLY:  diffusion_v
+
+    USE diffusion_w_mod,                                                       &
+        ONLY:  diffusion_w
 
     USE indices,                                                               &
@@ -334,59 +380,11 @@
                nzb, nzt, wall_flags_0
 
-    USE advec_ws,                                                              &
-        ONLY:  advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws
-
-    USE advec_s_bc_mod,                                                        &
-        ONLY:  advec_s_bc
-
-    USE advec_s_pw_mod,                                                        &
-        ONLY:  advec_s_pw
-
-    USE advec_s_up_mod,                                                        &
-        ONLY:  advec_s_up
-
-    USE advec_u_pw_mod,                                                        &
-        ONLY:  advec_u_pw
-
-    USE advec_u_up_mod,                                                        &
-        ONLY:  advec_u_up
-
-    USE advec_v_pw_mod,                                                        &
-        ONLY:  advec_v_pw
-
-    USE advec_v_up_mod,                                                        &
-        ONLY:  advec_v_up
-
-    USE advec_w_pw_mod,                                                        &
-        ONLY:  advec_w_pw
-
-    USE advec_w_up_mod,                                                        &
-        ONLY:  advec_w_up
-
-    USE buoyancy_mod,                                                          &
-        ONLY:  buoyancy
-
-    USE calc_radiation_mod,                                                    &
-        ONLY:  calc_radiation
-
-    USE coriolis_mod,                                                          &
-        ONLY:  coriolis
-
-    USE diffusion_s_mod,                                                       &
-        ONLY:  diffusion_s
-
-    USE diffusion_u_mod,                                                       &
-        ONLY:  diffusion_u
-
-    USE diffusion_v_mod,                                                       &
-        ONLY:  diffusion_v
-
-    USE diffusion_w_mod,                                                       &
-        ONLY:  diffusion_w
-
     USE kinds
 
     USE lsf_nudging_mod,                                                       &
         ONLY:  ls_advec, nudge
+
+    USE ocean_mod,                                                             &
+        ONLY:  eqn_state_seawater, ocean_prognostic_equations
 
     USE plant_canopy_model_mod,                                                &
@@ -403,4 +401,8 @@
         ONLY:  subsidence
 
+    USE surface_mod,                                                           &
+        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
+                surf_usm_v
+
     USE turbulence_closure_mod,                                                &
         ONLY:  tcm_prognostic
@@ -408,8 +410,4 @@
     USE user_actions_mod,                                                      &
         ONLY:  user_actions
-
-    USE surface_mod,                                                           &
-        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
-                surf_usm_v 
 
     USE wind_turbine_model_mod,                                                &
@@ -699,5 +697,5 @@
 
           IF ( .NOT. neutral )  THEN
-             IF ( ocean )  THEN
+             IF ( ocean_mode )  THEN
                 CALL buoyancy( i, j, rho_ocean, 3 )
              ELSE
@@ -844,71 +842,4 @@
                 ENDIF
              ENDIF
-
-          ENDIF
-
-!
-!--       If required, compute prognostic equation for salinity
-          IF ( ocean )  THEN
-
-!
-!--          Tendency-terms for salinity
-             tend(:,j,i) = 0.0_wp
-             IF ( timestep_scheme(1:5) == 'runge' ) &
-             THEN
-                IF ( ws_scheme_sca )  THEN
-                    CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa,  &
-                                diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn  )
-                ELSE
-                    CALL advec_s_pw( i, j, sa )
-                ENDIF
-             ELSE
-                CALL advec_s_up( i, j, sa )
-             ENDIF
-             CALL diffusion_s( i, j, sa,                                       &
-                               surf_def_h(0)%sasws, surf_def_h(1)%sasws,       &
-                               surf_def_h(2)%sasws,                            &
-                               surf_lsm_h%sasws,    surf_usm_h%sasws,          & 
-                               surf_def_v(0)%sasws, surf_def_v(1)%sasws,       &
-                               surf_def_v(2)%sasws, surf_def_v(3)%sasws,       &
-                               surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws,       &
-                               surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws,       &
-                               surf_usm_v(0)%sasws, surf_usm_v(1)%sasws,       &
-                               surf_usm_v(2)%sasws, surf_usm_v(3)%sasws )
-
-             CALL user_actions( i, j, 'sa-tendency' )
-
-!
-!--          Prognostic equation for salinity
-             DO  k = nzb+1, nzt
-                sa_p(k,j,i) = sa(k,j,i) + ( dt_3d *                            &
-                                                  ( tsc(2) * tend(k,j,i) +     &
-                                                    tsc(3) * tsa_m(k,j,i) )    &
-                                                  - tsc(5) * rdf_sc(k)         &
-                                                   * ( sa(k,j,i) - sa_init(k) )&
-                                          ) * MERGE( 1.0_wp, 0.0_wp,           &
-                                                BTEST( wall_flags_0(k,j,i), 0 )&
-                                                   )
-                IF ( sa_p(k,j,i) < 0.0_wp )  sa_p(k,j,i) = 0.1_wp * 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+1, nzt
-                      tsa_m(k,j,i) = tend(k,j,i)
-                   ENDDO
-                ELSEIF ( intermediate_timestep_count < &
-                         intermediate_timestep_count_max )  THEN
-                   DO  k = nzb+1, nzt
-                      tsa_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
-                                        5.3125_wp * tsa_m(k,j,i)
-                   ENDDO  
-                ENDIF
-             ENDIF
-
-!
-!--          Calculate density by the equation of state for seawater
-             CALL eqn_state_seawater( i, j )
 
           ENDIF
@@ -1324,11 +1255,12 @@
 
 !
-!--       If required, compute prognostic equation for chemical quantites
+!--       Calculate prognostic equation for chemical quantites
           IF ( air_chemistry )  THEN
+             !> TODO: remove time measurement since it slows down performance because it will be called extremely often
              CALL cpu_log( log_point(83), '(chem advec+diff+prog)', 'start' )
 !
 !--          Loop over chemical species
              DO  lsp = 1, nvar                         
-                CALL chem_prognostic_equations ( chem_species(lsp)%conc_p,     &
+                CALL chem_prognostic_equations( chem_species(lsp)%conc_p,      &
                                      chem_species(lsp)%conc,                   &
                                      chem_species(lsp)%tconc_m,                &
@@ -1342,8 +1274,13 @@
 
              CALL cpu_log( log_point(83), '(chem advec+diff+prog)', 'stop' )             
-          ENDIF   ! Chemicals equations                     
-
-       ENDDO
-    ENDDO
+          ENDIF   ! Chemical equations
+!
+!--       Calculate prognostic equations for the ocean
+          IF ( ocean_mode )  THEN
+             CALL ocean_prognostic_equations( i, j, i_omp_start, tn )
+          ENDIF
+
+       ENDDO  ! loop over j
+    ENDDO  ! loop over i
     !$OMP END PARALLEL
 
@@ -1576,5 +1513,5 @@
 
     IF ( .NOT. neutral )  THEN
-       IF ( ocean )  THEN
+       IF ( ocean_mode )  THEN
           CALL buoyancy( rho_ocean, 3 )
        ELSE
@@ -1770,105 +1707,4 @@
 
 !
-!-- If required, compute prognostic equation for salinity
-    IF ( ocean )  THEN
-
-       CALL cpu_log( log_point(37), 'sa-equation', 'start' )
-
-!
-!--    sa-tendency terms with communication
-       sbt = tsc(2)
-       IF ( scalar_advec == 'bc-scheme' )  THEN
-
-          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
-!
-!--          Bott-Chlond scheme always uses Euler time step. Thus:
-             sbt = 1.0_wp
-          ENDIF
-          tend = 0.0_wp
-          CALL advec_s_bc( sa, 'sa' )
-
-       ENDIF
-
-!
-!--    sa-tendency terms with no communication
-       IF ( scalar_advec /= 'bc-scheme' )  THEN
-          tend = 0.0_wp
-          IF ( timestep_scheme(1:5) == 'runge' )  THEN
-             IF ( ws_scheme_sca )  THEN
-                 CALL advec_s_ws( sa, 'sa' )
-             ELSE
-                 CALL advec_s_pw( sa )
-             ENDIF
-          ELSE
-             CALL advec_s_up( sa )
-          ENDIF
-       ENDIF
-
-       CALL diffusion_s( sa,                                                   &
-                         surf_def_h(0)%sasws, surf_def_h(1)%sasws,             &
-                         surf_def_h(2)%sasws,                                  &
-                         surf_lsm_h%sasws,    surf_usm_h%sasws,                &
-                         surf_def_v(0)%sasws, surf_def_v(1)%sasws,             &
-                         surf_def_v(2)%sasws, surf_def_v(3)%sasws,             &
-                         surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws,             &
-                         surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws,             &
-                         surf_usm_v(0)%sasws, surf_usm_v(1)%sasws,             &
-                         surf_usm_v(2)%sasws, surf_usm_v(3)%sasws )
-
-       CALL user_actions( 'sa-tendency' )
-
-!
-!--    Prognostic equation for salinity
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt
-                sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +      &
-                                                   tsc(3) * tsa_m(k,j,i) )     &
-                                                 - tsc(5) * rdf_sc(k) *        &
-                                                 ( sa(k,j,i) - sa_init(k) )    &
-                                          )                                    &
-                                   * MERGE( 1.0_wp, 0.0_wp,                    &
-                                             BTEST( wall_flags_0(k,j,i), 0 )   &
-                                          )
-                IF ( sa_p(k,j,i) < 0.0_wp )  sa_p(k,j,i) = 0.1_wp * 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+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+1, nzt
-                      tsa_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
-                                        5.3125_wp * tsa_m(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-       ENDIF
-
-       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
-
-!
 !-- If required, compute prognostic equation for total water content
     IF ( humidity )  THEN
@@ -2074,6 +1910,6 @@
 
           CALL cpu_log( log_point(67), 'qc-equation', 'stop' )
+
           CALL cpu_log( log_point(68), 'nc-equation', 'start' )
-
 !
 !--       Calculate prognostic equation for cloud drop concentration
@@ -2461,8 +2297,10 @@
     ENDIF
 
+!
+!-- Calculate prognostic equations for turbulence closure
     CALL tcm_prognostic()
 
 !
-!-- If required, compute prognostic equation for chemical quantites
+!-- Calculate prognostic equation for chemical quantites
     IF ( air_chemistry )  THEN
        CALL cpu_log( log_point(83), '(chem advec+diff+prog)', 'start' )
@@ -2470,9 +2308,9 @@
 !--    Loop over chemical species
        DO  lsp = 1, nvar                         
-          CALL chem_prognostic_equations ( chem_species(lsp)%conc_p,           &
-                                           chem_species(lsp)%conc,             &
-                                           chem_species(lsp)%tconc_m,          &
-                                           chem_species(lsp)%conc_pr_init,     &
-                                           lsp )       
+          CALL chem_prognostic_equations( chem_species(lsp)%conc_p,            &
+                                          chem_species(lsp)%conc,              &
+                                          chem_species(lsp)%tconc_m,           &
+                                          chem_species(lsp)%conc_pr_init,      &
+                                          lsp )
        ENDDO
 
@@ -2480,4 +2318,7 @@
     ENDIF   ! Chemicals equations
 
+!
+!-- Calculate prognostic equations for the ocean
+    IF ( ocean_mode )  CALL ocean_prognostic_equations()
 
  END SUBROUTINE prognostic_equations_vector
Index: palm/trunk/SOURCE/read_restart_data_mod.f90
===================================================================
--- palm/trunk/SOURCE/read_restart_data_mod.f90	(revision 3291)
+++ palm/trunk/SOURCE/read_restart_data_mod.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! changes concerning modularization of ocean option
+! 
+! 3289 2018-09-28 10:23:58Z suehring
 ! Introduce module parameter for number of inflow profiles
 ! 
@@ -77,5 +80,4 @@
 
     USE control_parameters
-     
 
     IMPLICIT NONE
@@ -99,10 +101,10 @@
 
 
-    PUBLIC rrd_global, rrd_read_parts_of_global, rrd_local,      & 
-           rrd_skip_global
+    PUBLIC rrd_global, rrd_read_parts_of_global, rrd_local, rrd_skip_global
 
 
  CONTAINS
 
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
@@ -115,5 +117,5 @@
        USE arrays_3d,                                                          &
            ONLY:  inflow_damping_factor, mean_inflow_profiles, pt_init,        &
-                  q_init, ref_state, s_init, sa_init, u_init, ug, v_init, vg
+                  q_init, ref_state, sa_init, s_init, u_init, ug, v_init, vg
 
        USE bulk_cloud_model_mod,                                               &
@@ -140,4 +142,7 @@
        USE netcdf_interface,                                                   &
            ONLY:  netcdf_precision, output_for_t0
+
+       USE ocean_mod,                                                          &
+           ONLY:  ocean_rrd_global
 
        USE particle_attributes,                                                &
@@ -177,7 +182,6 @@
 
        CALL check_open( 13 )
-
-!
-!-- Make version number check first
+!
+!--    Make version number check first
        READ ( 13 )  length
        READ ( 13 )  restart_string(1:length)
@@ -196,6 +200,6 @@
 
 !
-!-- Read number of PEs and horizontal index bounds of all PEs used in previous
-!-- run
+!--    Read number of PEs and horizontal index bounds of all PEs used in the
+!--    previous run
        READ ( 13 )  length
        READ ( 13 )  restart_string(1:length)
@@ -223,7 +227,7 @@
 
 !
-!-- Read vertical number of gridpoints and number of different areas used
-!-- for computing statistics. Allocate arrays depending on these values,
-!-- which are needed for the following read instructions.
+!--    Read vertical number of gridpoints and number of different areas used
+!--    for computing statistics. Allocate arrays depending on these values,
+!--    which are needed for the following read instructions.
        READ ( 13 )  length
        READ ( 13 )  restart_string(1:length)
@@ -257,4 +261,5 @@
        ENDIF
        READ ( 13 )  statistic_regions
+
        IF ( .NOT. ALLOCATED( ug ) )  THEN
           ALLOCATE( ug(0:nz+1), u_init(0:nz+1), vg(0:nz+1),                    &
@@ -266,9 +271,9 @@
 
 !
-!-- Now read all control parameters:
-!-- Caution: When the following read instructions have been changed, the
-!-- -------  version number stored in the variable binary_version_global has to
-!--          be increased. The same changes must also be done in 
-!--          wrd_write_global.
+!--    Now read all control parameters:
+!--    Caution: When the following read instructions have been changed, the
+!--    -------  version number stored in the variable binary_version_global has
+!--             to be increased. The same changes must also be done in
+!--             wrd_write_global.
        READ ( 13 )  length
        READ ( 13 )  restart_string(1:length)
@@ -318,12 +323,8 @@
              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
@@ -496,6 +497,6 @@
                 READ ( 13 )  ny
                 ny_on_file = ny
-             CASE ( 'ocean' )
-                READ ( 13 )  ocean
+             CASE ( 'ocean_mode' )
+                READ ( 13 )  ocean_mode
              CASE ( 'old_dt' )
                 READ ( 13 )  old_dt
@@ -580,12 +581,4 @@
              CASE ( 's_vertical_gradient_level_ind' )
                 READ ( 13 )  s_vertical_gradient_level_ind
-             CASE ( 'sa_init' )
-                READ ( 13 )  sa_init
-             CASE ( 'sa_surface' )
-                READ ( 13 )  sa_surface
-             CASE ( 'sa_vertical_gradient' )
-                READ ( 13 )  sa_vertical_gradient
-             CASE ( 'sa_vertical_gradient_level' )
-                READ ( 13 )  sa_vertical_gradient_level
              CASE ( 'scalar_advec' )
                 READ ( 13 )  scalar_advec
@@ -662,6 +655,4 @@
              CASE ( 'top_momentumflux_v' )
                 READ ( 13 )  top_momentumflux_v
-             CASE ( 'top_salinityflux' )
-                READ ( 13 )  top_salinityflux
              CASE ( 'top_scalarflux' )
                 READ ( 13 )  top_scalarflux
@@ -754,6 +745,4 @@
              CASE ( 'wall_humidityflux' )
                 READ ( 13 )  wall_humidityflux
-             CASE ( 'wall_salinityflux' )
-                READ ( 13 )  wall_salinityflux
              CASE ( 'wall_scalarflux' )
                 READ ( 13 )  wall_scalarflux
@@ -769,18 +758,16 @@
                 READ ( 13 )  z_i
 
-
              CASE DEFAULT
-
-                IF ( .NOT. found ) CALL bcm_rrd_global( found )
-
-                IF ( .NOT. found ) CALL wtm_rrd_global( found )
-
-                IF ( .NOT. found ) CALL flight_rrd_global( found )
-
-                IF ( .NOT. found ) CALL stg_rrd_global ( found )
-
-                IF ( .NOT. found ) CALL gust_rrd_global( found )
-
-                IF ( .NOT. found ) CALL user_rrd_global( found )
+!
+!--             Read global variables from of other modules
+                IF ( .NOT. found )  CALL bcm_rrd_global( found )
+                IF ( .NOT. found )  CALL flight_rrd_global( found )
+                IF ( .NOT. found )  CALL gust_rrd_global( found )
+                IF ( .NOT. found )  CALL ocean_rrd_global( found )
+                IF ( .NOT. found )  CALL stg_rrd_global ( found )
+                IF ( .NOT. found )  CALL wtm_rrd_global( found )
+!
+!--             Read user-defined global variables
+                IF ( .NOT. found )  CALL user_rrd_global( found )
 
                 IF ( .NOT. found )  THEN
@@ -1049,7 +1036,9 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  e, kh, km, p, pt, q, ql, &
-               s, sa, u, u_m_l, u_m_n, u_m_r, u_m_s, v, v_m_l, v_m_n,      &
-               v_m_r, v_m_s, vpt, w, w_m_l, w_m_n, w_m_r, w_m_s
+        ONLY:  e, kh, km, p, pt, q, ql, s, u, u_m_l, u_m_n, u_m_r, u_m_s,      &
+               v, v_m_l, v_m_n, v_m_r, v_m_s, vpt, w, w_m_l, w_m_n, w_m_r, w_m_s
+
+    USE bulk_cloud_model_mod,                                                  &
+        ONLY :  bcm_rrd_local
 
     USE averaging
@@ -1073,6 +1062,6 @@
         ONLY:  lsm_rrd_local
 
-    USE bulk_cloud_model_mod,                                                  &
-        ONLY :  bcm_rrd_local
+    USE ocean_mod,                                                             &
+        ONLY:  ocean_rrd_local
 
     USE particle_attributes,                                                   &
@@ -1599,12 +1588,4 @@
                    DEALLOCATE( tmp_2d_id_random, tmp_2d_seq_random )
 
-                CASE ( 'rho_ocean_av' )
-                   IF ( .NOT. ALLOCATED( rho_ocean_av ) )  THEN
-                      ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   rho_ocean_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
-                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-                           
                 CASE ( 's' )
                    IF ( k == 1 )  READ ( 13 )  tmp_3d
@@ -1618,17 +1599,4 @@
                    IF ( k == 1 )  READ ( 13 )  tmp_3d
                    s_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =           &
-                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'sa' )
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   sa(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
-                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'sa_av' )
-                   IF ( .NOT. ALLOCATED( sa_av ) )  THEN
-                      ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   sa_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
                       tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -1884,5 +1852,5 @@
 
 !
-!--                Read microphysics module restart data
+!--                Read restart data of other modules
                    IF ( .NOT. found ) CALL bcm_rrd_local( i, k, nxlf,          &
                                            nxlc, nxl_on_file, nxrf, nxrc,      &
@@ -1891,22 +1859,16 @@
                                            nys_on_file, tmp_2d, tmp_3d, found )
 
-!
-!--                Read surface related variables
-                   IF ( .NOT. found ) CALL surface_rrd_local( i, k, nxlf,      &
-                                           nxlc, nxl_on_file, nxrf, nxrc,      & 
+                   IF ( .NOT. found ) CALL chem_rrd_local( i, k, nxlf,         &
+                                           nxlc, nxl_on_file, nxrf, nxrc,      &
                                            nxr_on_file, nynf, nync,            & 
                                            nyn_on_file, nysf, nysc,            &
-                                           nys_on_file, found )
-
-!
-!--                Read urban surface restart data
-                   IF ( .NOT. found ) CALL usm_rrd_local( i, k, nxlf,          & 
+                                           nys_on_file, tmp_3d, found )
+
+                   IF ( .NOT. found ) CALL gust_rrd_local( i, k, nxlf,         &
                                            nxlc, nxl_on_file, nxrf, nxrc,      &
                                            nxr_on_file, nynf, nync,            &
                                            nyn_on_file, nysf, nysc,            &
-                                           nys_on_file, found )
-
-!
-!--                Read land surface restart data
+                                           nys_on_file, tmp_2d, tmp_3d, found )
+
                    IF ( .NOT. found ) CALL lsm_rrd_local( i, k, nxlf,          & 
                                            nxlc, nxl_on_file, nxrf, nxrc,      &
@@ -1915,6 +1877,10 @@
                                            nys_on_file, tmp_2d, found )
 
-!
-!--                Read radiation restart data
+                   IF ( .NOT. found ) CALL ocean_rrd_local( i, k, nxlf,        &
+                                           nxlc, nxl_on_file, nxrf, nxrc,      &
+                                           nxr_on_file, nynf, nync,            &
+                                           nyn_on_file, nysf, nysc,            &
+                                           nys_on_file, tmp_2d, tmp_3d, found )
+
                    IF ( .NOT. found ) CALL radiation_rrd_local( i, k, nxlf,    &
                                            nxlc, nxl_on_file, nxrf, nxrc,      &
@@ -1923,19 +1889,15 @@
                                            nys_on_file, tmp_2d, tmp_3d, found )
 
-!
-!--                Read chemistry restart data
-                   IF ( .NOT. found ) CALL chem_rrd_local( i, k, nxlf,         &
+                   IF ( .NOT. found ) CALL surface_rrd_local( i, k, nxlf,      &
                                            nxlc, nxl_on_file, nxrf, nxrc,      &
                                            nxr_on_file, nynf, nync,            & 
                                            nyn_on_file, nysf, nysc,            &
-                                           nys_on_file, tmp_3d, found )
-
-!
-!--                Read gust module restart data
-                   IF ( .NOT. found ) CALL gust_rrd_local( i, k, nxlf,         &
+                                           nys_on_file, found )
+
+                   IF ( .NOT. found ) CALL usm_rrd_local( i, k, nxlf,          &
                                            nxlc, nxl_on_file, nxrf, nxrc,      &
                                            nxr_on_file, nynf, nync,            &
                                            nyn_on_file, nysf, nysc,            &
-                                           nys_on_file, tmp_2d, tmp_3d, found )
+                                           nys_on_file, found )
 
 !
Index: palm/trunk/SOURCE/subsidence_mod.f90
===================================================================
--- palm/trunk/SOURCE/subsidence_mod.f90	(revision 3291)
+++ palm/trunk/SOURCE/subsidence_mod.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! ocean renamed ocean_mode
+! 
+! 3045 2018-05-28 07:55:41Z Giersch
 ! Error message revised
 ! 
@@ -126,5 +129,5 @@
 
        USE control_parameters,                                                 &
-           ONLY:  message_string, ocean, subs_vertical_gradient,               &
+           ONLY:  message_string, ocean_mode, subs_vertical_gradient,          &
                   subs_vertical_gradient_level, subs_vertical_gradient_level_i
 
@@ -147,7 +150,7 @@
        ENDIF 
 
-       IF ( ocean )  THEN
+       IF ( ocean_mode )  THEN
           message_string = 'Applying large scale vertical motion is not ' //   &
-                           'allowed for ocean runs'
+                           'allowed for ocean mode'
           CALL message( 'init_w_subsidence', 'PA0324', 2, 2, 0, 6, 0 )
        ENDIF
Index: palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 3291)
+++ palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! changes concerning modularization of ocean option
+! 
+! 3291 2018-09-28 11:33:03Z scharf
 ! corrected previous commit for 3D topography
 ! 
@@ -211,14 +214,13 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  dzw, e, heatflux_output_conversion, nc, nr, p, prr, pt,         &
-               q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, sa, u, v, vpt, w,      &
-               waterflux_output_conversion, d_exner
+        ONLY:  dzw, d_exner, e, heatflux_output_conversion, nc, nr, p, prr,    &
+               pt, q, qc, ql, ql_c, ql_v, qr, s, u, v, vpt, w,                 &
+               waterflux_output_conversion
 
     USE averaging,                                                             &
-        ONLY:  e_av, ghf_av, lpt_av, lwp_av,                                   &
-               ol_av, p_av, pc_av, pr_av, pt_av, q_av, ql_av,                  &
-               ql_c_av, ql_v_av, ql_vp_av, qsws_av, qv_av, r_a_av,             &
-               rho_ocean_av, s_av, sa_av, shf_av, ssws_av, ts_av, tsurf_av,    &
-               u_av, us_av, v_av, vpt_av, w_av, z0_av, z0h_av, z0q_av
+        ONLY:  e_av, ghf_av, lpt_av, lwp_av, ol_av, p_av, pc_av, pr_av, pt_av, &
+               q_av, ql_av, ql_c_av, ql_v_av, ql_vp_av, qsws_av, qv_av,        &
+               r_a_av, s_av, shf_av, ssws_av, ts_av, tsurf_av, u_av, us_av,    &
+               v_av, vpt_av, w_av, z0_av, z0h_av, z0q_av
 
     USE basic_constants_and_equations_mod,                                     &
@@ -232,6 +234,6 @@
 
     USE control_parameters,                                                    &
-        ONLY:  air_chemistry, average_count_3d, doav, doav_n,                  &
-               land_surface, rho_surface, urban_surface, uv_exposure,          &
+        ONLY:  air_chemistry, average_count_3d, doav, doav_n, land_surface,    &
+               ocean_mode, rho_surface, urban_surface, uv_exposure,            &
                varnamelength
 
@@ -249,4 +251,7 @@
     USE land_surface_model_mod,                                                &
         ONLY:  lsm_3d_data_averaging
+
+    USE ocean_mod,                                                             &
+        ONLY:  ocean_3d_data_averaging
 
     USE particle_attributes,                                                   &
@@ -412,10 +417,4 @@
                 r_a_av = 0.0_wp
 
-             CASE ( 'rho_ocean' )
-                IF ( .NOT. ALLOCATED( rho_ocean_av ) )  THEN
-                   ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                rho_ocean_av = 0.0_wp
-
              CASE ( 's' )
                 IF ( .NOT. ALLOCATED( s_av ) )  THEN
@@ -423,10 +422,4 @@
                 ENDIF
                 s_av = 0.0_wp
-
-             CASE ( 'sa' )
-                IF ( .NOT. ALLOCATED( sa_av ) )  THEN
-                   ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                sa_av = 0.0_wp
 
              CASE ( 'shf*' )
@@ -501,8 +494,8 @@
                 ENDIF
                 z0q_av = 0.0_wp
+
+             CASE ( 'usm_output' )
 !             
-!--          Block of urban surface model outputs 
-             CASE ( 'usm_output' )
-
+!--             Block of urban surface model outputs
                 CALL usm_average_3d_data( 'allocate', doav(ii) )
               
@@ -511,45 +504,37 @@
 
 !
-!--             Turbulence closure module
-                CALL tcm_3d_data_averaging( 'allocate', doav(ii) )
-
-!
-!--             Microphysics module quantities
+!--             Allocating and initializing data arrays for other modules
                 IF ( bulk_cloud_model )  THEN
                    CALL bcm_3d_data_averaging( 'allocate', doav(ii) )
                 ENDIF
 
-!
-!--             Land surface quantity
+                IF ( air_chemistry  .AND.  trimvar(1:3) == 'kc_')  THEN
+                   CALL chem_3d_data_averaging( 'allocate', doav(ii) )
+                ENDIF
+
+                IF ( gust_module_enabled )  THEN
+                   CALL gust_3d_data_averaging( 'allocate', doav(ii) )
+                ENDIF
+
                 IF ( land_surface )  THEN
                    CALL lsm_3d_data_averaging( 'allocate', doav(ii) )
                 ENDIF
 
-!
-!--             Radiation quantity
+                IF ( ocean_mode )  THEN
+                   CALL ocean_3d_data_averaging( 'allocate', doav(ii) )
+                ENDIF
+
                 IF ( radiation )  THEN
                    CALL radiation_3d_data_averaging( 'allocate', doav(ii) )
                 ENDIF
 
-!
-!--             Gust module quantities
-                IF ( gust_module_enabled )  THEN
-                   CALL gust_3d_data_averaging( 'allocate', doav(ii) )
-                ENDIF
-
-!
-!--             Chemical quantity                                           
-                IF ( air_chemistry  .AND.  trimvar(1:3) == 'kc_')  THEN
-                   CALL chem_3d_data_averaging( 'allocate', doav(ii) )
-                ENDIF
-
-!
-!--             UV exposure quantity
                 IF ( uv_exposure  .AND.  trimvar(1:5) == 'uvem_')  THEN
                    CALL uvem_3d_data_averaging( 'allocate', doav(ii) )
                 ENDIF
 
-!
-!--             User-defined quantity
+                CALL tcm_3d_data_averaging( 'allocate', doav(ii) )
+
+!
+!--             User-defined quantities
                 CALL user_3d_data_averaging( 'allocate', doav(ii) )
 
@@ -877,15 +862,4 @@
              ENDIF
 
-          CASE ( 'rho_ocean' )
-             IF ( ALLOCATED( rho_ocean_av ) ) THEN 
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) + rho_ocean(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF  
-
           CASE ( 's' )
              IF ( ALLOCATED( s_av ) ) THEN 
@@ -894,15 +868,4 @@
                       DO  k = nzb, nzt+1
                          s_av(k,j,i) = s_av(k,j,i) + s(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
-          CASE ( 'sa' )
-             IF ( ALLOCATED( sa_av ) ) THEN 
-                DO  i = nxlg, nxrg
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         sa_av(k,j,i) = sa_av(k,j,i) + sa(k,j,i)
                       ENDDO
                    ENDDO
@@ -1182,10 +1145,11 @@
                 ENDDO
              ENDIF
+
+          CASE ( 'usm_output' )
 !             
-!--       Block of urban surface model outputs. 
-!--       In case of urban surface variables it should be always checked
-!--       if respective arrays are allocated, at least in case of a restart 
-!--       run, as averaged usm arrays are not read from file at the moment. 
-          CASE ( 'usm_output' )
+!--          Block of urban surface model outputs.
+!--          In case of urban surface variables it should be always checked
+!--          if respective arrays are allocated, at least in case of a restart
+!--          run, as averaged usm arrays are not read from file at the moment.
              CALL usm_average_3d_data( 'allocate', doav(ii) )
              CALL usm_average_3d_data( 'sum', doav(ii) )
@@ -1193,39 +1157,31 @@
           CASE DEFAULT
 !
-!--          Turbulence closure module
-             CALL tcm_3d_data_averaging( 'sum', doav(ii) )
-
-!
-!--          Microphysics module quantities
+!--          Summing up data from other modules
              IF ( bulk_cloud_model )  THEN
                 CALL bcm_3d_data_averaging( 'sum', doav(ii) )
              ENDIF
 
-!
-!--          Land surface quantity
+             IF ( air_chemistry  .AND.  trimvar(1:3) == 'kc_')  THEN
+                CALL chem_3d_data_averaging( 'sum',doav(ii) )
+             ENDIF
+
+             IF ( gust_module_enabled )  THEN
+                CALL gust_3d_data_averaging( 'sum', doav(ii) )
+             ENDIF
+
              IF ( land_surface )  THEN
                 CALL lsm_3d_data_averaging( 'sum', doav(ii) )
              ENDIF
 
-!
-!--          Radiation quantity
+             IF ( ocean_mode )  THEN
+                CALL ocean_3d_data_averaging( 'sum', doav(ii) )
+             ENDIF
+
              IF ( radiation )  THEN
                 CALL radiation_3d_data_averaging( 'sum', doav(ii) )
              ENDIF
 
-!
-!--          Gust module quantities
-             IF ( gust_module_enabled )  THEN
-                CALL gust_3d_data_averaging( 'sum', doav(ii) )
-             ENDIF
-
-!
-!--          Chemical quantity
-             IF ( air_chemistry  .AND.  trimvar(1:3) == 'kc_')  THEN
-                CALL chem_3d_data_averaging( 'sum',doav(ii) )
-             ENDIF
-
-!
-!--          UV exposure quantity
+             CALL tcm_3d_data_averaging( 'sum', doav(ii) )
+
              IF ( uv_exposure )  THEN
                 CALL uvem_3d_data_averaging( 'sum', doav(ii) )
@@ -1233,5 +1189,5 @@
 
 !
-!--          User-defined quantity
+!--          User-defined quantities
              CALL user_3d_data_averaging( 'sum', doav(ii) )
 
Index: palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_mod.f90	(revision 3291)
+++ palm/trunk/SOURCE/surface_mod.f90	(revision 3294)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! changes concerning modularization of ocean option
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -1114,5 +1117,5 @@
 !
 !--    Salinity surface flux
-       IF ( ocean )  DEALLOCATE ( surfaces%sasws )
+       IF ( ocean_mode )  DEALLOCATE ( surfaces%sasws )
 
     END SUBROUTINE deallocate_surface_attributes_h
@@ -1241,5 +1244,5 @@
 !
 !--    Salinity surface flux
-       IF ( ocean )  ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
+       IF ( ocean_mode )  ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
 
     END SUBROUTINE allocate_surface_attributes_h
@@ -1306,5 +1309,5 @@
 !
 !--    Salinity flux
-       IF ( ocean )  DEALLOCATE ( surfaces%sasws )
+       IF ( ocean_mode )  DEALLOCATE ( surfaces%sasws )
 
     END SUBROUTINE deallocate_surface_attributes_h_top
@@ -1378,5 +1381,5 @@
 !
 !--    Salinity flux
-       IF ( ocean )  ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
+       IF ( ocean_mode )  ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
 
     END SUBROUTINE allocate_surface_attributes_h_top
@@ -1496,5 +1499,5 @@
 !
 !--    Salinity surface flux
-       IF ( ocean )  DEALLOCATE ( surfaces%sasws )
+       IF ( ocean_mode )  DEALLOCATE ( surfaces%sasws )
 
     END SUBROUTINE deallocate_surface_attributes_v
@@ -1621,5 +1624,5 @@
 !
 !--    Salinity surface flux
-       IF ( ocean )  ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
+       IF ( ocean_mode )  ALLOCATE ( surfaces%sasws(1:surfaces%ns) )
 
     END SUBROUTINE allocate_surface_attributes_v
@@ -2268,5 +2271,5 @@
                 ENDIF
 
-                IF ( ocean )  THEN
+                IF ( ocean_mode )  THEN
                    IF ( upward_facing )  THEN 
                       surf%sasws(num_h) = bottom_salinityflux * rho_air_zw(k-1)
@@ -2344,5 +2347,5 @@
 !
 !--          Prescribe top salinity flux
-             IF ( ocean .AND. constant_top_salinityflux)                       &
+             IF ( ocean_mode .AND. constant_top_salinityflux)                  &
                 surf%sasws(num_h) = top_salinityflux * rho_air_zw(nzt+1)
 !
@@ -2488,5 +2491,5 @@
 !--          So far, salinityflux at vertical surfaces is simply zero 
 !--          at the moment  
-             IF ( ocean )  surf%sasws(num_v) = wall_salinityflux(component)
+             IF ( ocean_mode )  surf%sasws(num_v) = wall_salinityflux(component)
 !
 !--          Increment wall indices
Index: palm/trunk/SOURCE/swap_timelevel.f90
===================================================================
--- palm/trunk/SOURCE/swap_timelevel.f90	(revision 3291)
+++ palm/trunk/SOURCE/swap_timelevel.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! changes concerning modularization of ocean option
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -147,13 +150,10 @@
         ONLY:  chem_swap_timelevel 
 
-    USE land_surface_model_mod,                                                &
-        ONLY: lsm_swap_timelevel
-
     USE cpulog,                                                                &
         ONLY: cpu_log, log_point
 
     USE control_parameters,                                                    &
-        ONLY:  air_chemistry, humidity, land_surface,                          &
-               neutral, ocean, passive_scalar, timestep_count, urban_surface
+        ONLY:  air_chemistry, humidity, land_surface, neutral, ocean_mode,     &
+               passive_scalar, timestep_count, urban_surface
 
     USE gust_mod,                                                              &
@@ -164,4 +164,10 @@
         ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt
 #endif
+
+    USE land_surface_model_mod,                                                &
+        ONLY: lsm_swap_timelevel
+
+    USE ocean_mod,                                                             &
+        ONLY:  ocean_swap_timelevel
 
     USE pmc_interface,                                                         &
@@ -202,10 +208,4 @@
     ENDDO
 
-    CALL tcm_swap_timelevel ( 0 )
-
-    IF ( ocean )  THEN
-       sa = sa_p
-    ENDIF
-
     IF ( humidity )  THEN
        q = q_p
@@ -214,20 +214,12 @@
     IF ( passive_scalar )  s = s_p             
 
-    IF ( humidity  .AND.  bulk_cloud_model )  THEN
-       CALL bcm_swap_timelevel ( 0 )
-    ENDIF
-
-    IF ( land_surface )  THEN
-       CALL lsm_swap_timelevel ( 0 )
-    ENDIF
-
-    IF ( urban_surface )  THEN
-       CALL usm_swap_timelevel ( 0 )
-    ENDIF
-
-    IF ( gust_module_enabled )  THEN
-       CALL gust_swap_timelevel ( 0 )
-    ENDIF
-
+!
+!-- Swapping the timelevel of other modules
+    IF ( humidity  .AND.  bulk_cloud_model )  CALL bcm_swap_timelevel( 0 )
+    IF ( gust_module_enabled )                CALL gust_swap_timelevel( 0 )
+    IF ( land_surface )                       CALL lsm_swap_timelevel( 0 )
+    IF ( ocean_mode )                         CALL ocean_swap_timelevl( 0 )
+    CALL tcm_swap_timelevel( 0 )
+    IF ( urban_surface )                      CALL usm_swap_timelevel( 0 )
 
     CALL cpu_log( log_point(28), 'swap_timelevel (nop)', 'stop' )
@@ -245,5 +237,5 @@
              pt => pt_1;  pt_p => pt_2
           ENDIF
-          IF ( ocean )  THEN
+          IF ( ocean_mode )  THEN
              sa => sa_1;  sa_p => sa_2
           ENDIF
@@ -254,6 +246,4 @@
              s => s_1;    s_p => s_2
           ENDIF
-
-          IF ( air_chemistry )  CALL chem_swap_timelevel(0)
 
           swap_level = 1
@@ -267,5 +257,5 @@
              pt => pt_2;  pt_p => pt_1
           ENDIF
-          IF ( ocean )  THEN
+          IF ( ocean_mode )  THEN
              sa => sa_2;  sa_p => sa_1
           ENDIF
@@ -277,9 +267,9 @@
           ENDIF
 
-          IF ( air_chemistry )  CALL chem_swap_timelevel(1)
-
           swap_level = 2
 
     END SELECT
+
+    IF ( air_chemistry )  CALL chem_swap_timelevel( MOD( timestep_count, 2) )
 
     CALL tcm_swap_timelevel ( MOD( timestep_count, 2) )
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 3291)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! changes concerning modularization of ocean option
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -369,4 +372,8 @@
                ref_state, rho_ocean, s, s_p, sa_p, tend, u, u_p, v, vpt,       &
                v_p, w, w_p
+
+    USE bulk_cloud_model_mod,                                                  &
+        ONLY: bulk_cloud_model, calc_liquid_water_content,                     &
+              collision_turbulence, microphysics_morrison, microphysics_seifert
 
     USE calc_mean_profile_mod,                                                 &
@@ -399,5 +406,5 @@
                multi_agent_system_end, multi_agent_system_start,               &
                nesting_offline, neutral, nr_timesteps_this_run, nudging,       &
-               ocean, passive_scalar, prho_reference, pt_reference,            &
+               ocean_mode, passive_scalar, pt_reference,                       &
                pt_slope_offset, random_heatflux, rans_mode,                    &
                rans_tke_e, run_coupled, simulated_time, simulated_time_chr,    &
@@ -443,8 +450,4 @@
                lsf_nesting_offline, lsf_nesting_offline_mass_conservation
 
-    USE bulk_cloud_model_mod,                                                  &
-        ONLY: bulk_cloud_model, calc_liquid_water_content,                     &
-              collision_turbulence, microphysics_morrison, microphysics_seifert
-
     USE netcdf_data_input_mod,                                                 &
         ONLY:  nest_offl, netcdf_data_input_lsf
@@ -452,4 +455,7 @@
     USE multi_agent_system_mod,                                                &
         ONLY:  agents_active, multi_agent_system
+
+    USE ocean_mod,                                                             &
+        ONLY:  prho_reference
 
     USE particle_attributes,                                                   &
@@ -662,5 +668,5 @@
                 ref_state(:)  = hom(:,1,4,0) ! this is used in the buoyancy term
              ENDIF
-             IF ( ocean )  THEN
+             IF ( ocean_mode )  THEN
                 CALL calc_mean_profile( rho_ocean, 64 )
                 ref_state(:)  = hom(:,1,64,0)
@@ -753,5 +759,5 @@
              ENDIF
           ENDIF
-          IF ( ocean )  THEN
+          IF ( ocean_mode )  THEN
              CALL exchange_horiz( sa_p, nbgp )
              CALL exchange_horiz( rho_ocean, nbgp )
@@ -1053,5 +1059,5 @@
              CALL cpu_log( log_point(17), 'diffusivities', 'start' )
              IF ( .NOT. humidity ) THEN
-                IF ( ocean )  THEN
+                IF ( ocean_mode )  THEN
                    CALL tcm_diffusivities( prho, prho_reference )
                 ELSE
Index: palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3291)
+++ palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! changes concerning modularization of ocean option
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -164,6 +167,6 @@
                initializing_actions, intermediate_timestep_count,              &
                intermediate_timestep_count_max, km_constant,                   &
-               les_dynamic, les_mw, ocean, plant_canopy, prandtl_number,       &
-               prho_reference, pt_reference, rans_mode, rans_tke_e, rans_tke_l,&
+               les_dynamic, les_mw, ocean_mode, plant_canopy, prandtl_number,  &
+               pt_reference, rans_mode, rans_tke_e, rans_tke_l,                &
                simulated_time,timestep_scheme, turbulence_closure,             &
                turbulent_inflow, use_upstream_for_tke, vpt_reference,          &
@@ -190,4 +193,7 @@
 
     USE kinds
+
+    USE ocean_mod,                                                             &
+        ONLY:  prho_reference
 
     USE pegrid
@@ -1116,5 +1122,5 @@
                 e  = e_init
              ELSE
-                IF ( .NOT. ocean )  THEN
+                IF ( .NOT. ocean_mode )  THEN
                    kh   = 0.01_wp   ! there must exist an initial diffusion, because 
                    km   = 0.01_wp   ! otherwise no TKE would be produced by the
@@ -1160,5 +1166,5 @@
              e  = e_init
           ELSE
-             IF ( .NOT. ocean )  THEN
+             IF ( .NOT. ocean_mode )  THEN
                 kh   = 0.01_wp   ! there must exist an initial diffusion, because 
                 km   = 0.01_wp   ! otherwise no TKE would be produced by the
@@ -2155,5 +2161,5 @@
 
        IF ( .NOT. humidity )  THEN
-          IF ( ocean )  THEN
+          IF ( ocean_mode )  THEN
              CALL diffusion_e( prho, prho_reference )
           ELSE
@@ -2423,5 +2429,5 @@
 
        IF ( .NOT. humidity )  THEN
-          IF ( ocean )  THEN
+          IF ( ocean_mode )  THEN
              CALL diffusion_e( i, j, prho, prho_reference )
           ELSE
@@ -2934,5 +2940,5 @@
           IF ( .NOT. humidity )  THEN
 
-             IF ( ocean )  THEN
+             IF ( ocean_mode )  THEN
 !
 !--             So far in the ocean no special treatment of density flux
@@ -3738,5 +3744,5 @@
        IF ( .NOT. humidity )  THEN
 
-          IF ( ocean )  THEN
+          IF ( ocean_mode )  THEN
 !
 !--          So far in the ocean no special treatment of density flux in
Index: palm/trunk/SOURCE/write_restart_data_mod.f90
===================================================================
--- palm/trunk/SOURCE/write_restart_data_mod.f90	(revision 3291)
+++ palm/trunk/SOURCE/write_restart_data_mod.f90	(revision 3294)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! changes concerning modularization of ocean option
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -95,5 +98,8 @@
        USE arrays_3d,                                                          &
            ONLY:  inflow_damping_factor, mean_inflow_profiles, pt_init,        &
-                  q_init, ref_state, s_init, sa_init, u_init, ug, v_init, vg
+                  q_init, ref_state, s_init, u_init, ug, v_init, vg
+
+       USE bulk_cloud_model_mod,                                               &
+           ONLY:  bulk_cloud_model, bcm_wrd_global
 
        USE date_and_time_mod,                                                  &
@@ -112,7 +118,4 @@
            ONLY:  nx, ny, nz
 
-       USE bulk_cloud_model_mod,                                               &
-           ONLY:  bulk_cloud_model, bcm_wrd_global
-
        USE model_1d_mod,                                                       &
            ONLY:  damp_level_1d, dt_pr_1d, dt_run_control_1d, end_time_1d
@@ -120,4 +123,7 @@
        USE netcdf_interface,                                                   &
            ONLY:  netcdf_precision, output_for_t0
+
+       USE ocean_mod,                                                          &
+           ONLY:  ocean_wrd_global
 
        USE particle_attributes,                                                &
@@ -238,7 +244,4 @@
        WRITE ( 14 )  bc_s_t
 
-       CALL wrd_write_string( 'bc_sa_t' ) 
-       WRITE ( 14 )  bc_sa_t
-
        CALL wrd_write_string( 'bc_uv_b' ) 
        WRITE ( 14 )  bc_uv_b
@@ -246,7 +249,4 @@
        CALL wrd_write_string( 'bc_uv_t' ) 
        WRITE ( 14 )  bc_uv_t
-
-       CALL wrd_write_string( 'bottom_salinityflux' ) 
-       WRITE ( 14 )  bottom_salinityflux
 
        CALL wrd_write_string( 'building_height' ) 
@@ -497,6 +497,6 @@
        WRITE ( 14 )  ny
 
-       CALL wrd_write_string( 'ocean' ) 
-       WRITE ( 14 )  ocean
+       CALL wrd_write_string( 'ocean_mode' ) 
+       WRITE ( 14 )  ocean_mode
 
        CALL wrd_write_string( 'old_dt' ) 
@@ -622,16 +622,4 @@
        CALL wrd_write_string( 's_vertical_gradient_level_ind' ) 
        WRITE ( 14 )  s_vertical_gradient_level_ind
-
-       CALL wrd_write_string( 'sa_init' ) 
-       WRITE ( 14 )  sa_init
-
-       CALL wrd_write_string( 'sa_surface' ) 
-       WRITE ( 14 )  sa_surface
-
-       CALL wrd_write_string( 'sa_vertical_gradient' ) 
-       WRITE ( 14 )  sa_vertical_gradient
-
-       CALL wrd_write_string( 'sa_vertical_gradient_level' ) 
-       WRITE ( 14 )  sa_vertical_gradient_level
 
        CALL wrd_write_string( 'scalar_advec' ) 
@@ -738,7 +726,4 @@
        WRITE ( 14 )  top_momentumflux_v
 
-       CALL wrd_write_string( 'top_salinityflux' ) 
-       WRITE ( 14 )  top_salinityflux
-
        CALL wrd_write_string( 'top_scalarflux' ) 
        WRITE ( 14 )  top_scalarflux
@@ -876,7 +861,4 @@
        WRITE ( 14 )  wall_humidityflux
 
-       CALL wrd_write_string( 'wall_salinityflux' ) 
-       WRITE ( 14 )  wall_salinityflux
-
        CALL wrd_write_string( 'wall_scalarflux' ) 
        WRITE ( 14 )  wall_scalarflux
@@ -897,30 +879,16 @@
        WRITE ( 14 )  z_i
 
-
-!    
-!-- The order of the calls can be different as in read_restart_data_mod where
-!-- the restart data will be read 
-
-!
-!-- If required, write restart data for gust module.
-       IF ( bulk_cloud_model )  CALL bcm_wrd_global
-!
-!-- If required, write restart data for wind turbine model.
-       IF ( wind_turbine )  CALL wtm_wrd_global
-
-!
-!-- If required, write restart data for virtual measurements.
-       IF ( virtual_flight )  CALL flight_wrd_global
-
-!
-!-- If required, write restart data for turbulence generator.
-       IF ( syn_turb_gen )  CALL stg_wrd_global
-
-!
-!-- If required, write restart data for gust module.
+!
+!--    Write restart data of the other modules
+!--    The order of the calls can be different as in read_restart_data_mod where
+!--    the restart data will be read
+       IF ( bulk_cloud_model )     CALL bcm_wrd_global
+       IF ( virtual_flight )       CALL flight_wrd_global
        IF ( gust_module_enabled )  CALL gust_wrd_global
-
-!
-!-- Write restart data for user data.
+       IF ( ocean_mode )           CALL ocean_wrd_global
+       IF ( syn_turb_gen )         CALL stg_wrd_global
+       IF ( wind_turbine )         CALL wtm_wrd_global
+!
+!--    Write restart data of the user-interface
        CALL user_wrd_global
 
@@ -940,5 +908,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  e, kh, km, p, pt, q, ql, s, sa, u, u_m_l, u_m_n, u_m_r,      &
+           ONLY:  e, kh, km, p, pt, q, ql, s, u, u_m_l, u_m_n, u_m_r,          &
                   u_m_s, v, v_m_l, v_m_n, v_m_r, v_m_s, vpt, w, w_m_l, w_m_n,  &
                   w_m_r, w_m_s
@@ -960,4 +928,7 @@
        USE land_surface_model_mod,                                             &
            ONLY:  lsm_wrd_local
+
+       USE ocean_mod,                                                          &
+           ONLY:  ocean_wrd_local
 
        USE particle_attributes,                                                &
@@ -1112,21 +1083,4 @@
        ENDIF
        
-       IF ( ocean )  THEN
-
-          IF ( ALLOCATED( rho_ocean_av ) )  THEN
-             CALL wrd_write_string( 'rho_ocean_av' ) 
-             WRITE ( 14 )  rho_ocean_av
-          ENDIF
-
-          CALL wrd_write_string( 'sa' ) 
-          WRITE ( 14 )  sa
-
-          IF ( ALLOCATED( sa_av ) )  THEN
-             CALL wrd_write_string( 'sa_av' ) 
-             WRITE ( 14 )  sa_av
-          ENDIF
-
-       ENDIF
-
        IF ( ALLOCATED( ql_c_av ) )  THEN
           CALL wrd_write_string( 'ql_c_av' ) 
@@ -1310,48 +1264,20 @@
 
 !
-!-- If required, write restart data for microphysics module.
-       IF ( bulk_cloud_model )  THEN
-          CALL bcm_wrd_local
-       ENDIF
-
-!
-!-- Write surface-related restart data.
+!--    Write restart data of other modules
+       IF ( bulk_cloud_model )     CALL bcm_wrd_local
+       IF ( air_chemistry )        CALL chem_wrd_local
+       IF ( gust_module_enabled )  CALL gust_wrd_local
+       IF ( land_surface )         CALL lsm_wrd_local
+       IF ( ocean_mode )           CALL ocean_wrd_local
        CALL surface_wrd_local
-
-!
-!-- If required, write restart data for urban surface model.
-       IF ( urban_surface )  THEN
-          CALL usm_wrd_local
-       ENDIF
-
-!
-!-- If required, write restart data for land surface model.
-       IF ( land_surface )  THEN
-          CALL lsm_wrd_local
-       ENDIF
-
-!
-!-- If required, write restart data for radiation model.
-       IF ( radiation )  THEN
-          CALL radiation_wrd_local
-       ENDIF
-
-!
-!-- If required, write restart data for chemistry model.
-       IF ( air_chemistry )  THEN
-          CALL chem_wrd_local
-       ENDIF
-
-!
-!-- If required, write restart data for gust module.
-       IF ( gust_module_enabled )  THEN
-          CALL gust_wrd_local
-       ENDIF
-
-!
-!-- Write user-related restart data.
+       IF ( radiation )            CALL radiation_wrd_local
+       IF ( urban_surface )        CALL usm_wrd_local
+
+!
+!--    Write user-defined restart data
        CALL user_wrd_local
 
-!--    Write end label.
+!
+!--    Write end label
        CALL wrd_write_string( '*** end ***' )
 
