Index: palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 377)
+++ palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 388)
@@ -125,4 +125,10 @@
 Errors:
 ------
+Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
+calculated in 5 iteration steps. (init_ocean)
+
+Bugfix: wrong sign in buoyancy production of ocean part in case of not using
+the reference density (only in 3D routine production_e) (production_e)
+
 Bugfix: output of averaged 2d/3d quantities requires that an avaraging
 interval has been set, respective error message is included (check_parameters)
@@ -183,3 +189,3 @@
 Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)
 
-advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list
+advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 377)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 388)
@@ -4,4 +4,5 @@
 ! Actual revisions:
 ! -----------------
+! Check profiles fpr prho and hyp.
 ! Bugfix: output of averaged 2d/3d quantities requires that an avaraging
 ! interval has been set, respective error message is included
@@ -933,5 +934,6 @@
 
 !
-!--    If required compute the profile of leaf area density used in the plant canopy model 
+!--    If required compute the profile of leaf area density used in the plant
+!--    canopy model 
        IF ( plant_canopy ) THEN
        
@@ -967,5 +969,6 @@
 
 !
-!--       In case of no given leaf area density gradients, choose a vanishing gradient
+!--       In case of no given leaf area density gradients, choose a vanishing
+!--       gradient
           IF ( lad_vertical_gradient_level(1) == -9999999.9 ) THEN
              lad_vertical_gradient_level(1) = 0.0
@@ -2089,7 +2092,14 @@
 
           CASE ( 'rho' )
-             dopr_index(i) = 64
-             dopr_unit(i)  = 'kg/m3'
-             hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
+             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 )
+             ENDIF
 
           CASE ( 'w"sa"' )
@@ -2150,4 +2160,21 @@
                 hom(:,2,70,:) = SPREAD( zu, 2, statistic_regions+1 )
              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
+             dopr_unit(i)  = 'kPa'
+             hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
 
           CASE DEFAULT
@@ -2929,6 +2956,6 @@
 !-- Check pressure gradient conditions
     IF ( dp_external .AND. conserve_volume_flow )  THEN
-       WRITE( message_string, * )  'Both dp_external and conserve_volume_flow', &
-            ' are .TRUE. but one of them must be .FALSE.'
+       WRITE( message_string, * )  'Both dp_external and conserve_volume_flo', &
+            'w are .TRUE. but one of them must be .FALSE.'
        CALL message( 'check_parameters', 'PA0150', 1, 2, 0, 6, 0 )
     ENDIF
@@ -2940,6 +2967,6 @@
        ENDIF
        IF ( .NOT. ANY( dpdxy /= 0.0 ) )  THEN
-          WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is zero',&
-               ', i.e. the external pressure gradient & will not be applied'
+          WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is ze', &
+               'ro, i.e. the external pressure gradient & will not be applied'
           CALL message( 'check_parameters', 'PA0152', 0, 1, 0, 6, 0 )
        ENDIF
Index: palm/trunk/SOURCE/eqn_state_seawater.f90
===================================================================
--- palm/trunk/SOURCE/eqn_state_seawater.f90	(revision 377)
+++ palm/trunk/SOURCE/eqn_state_seawater.f90	(revision 388)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! First constant in array den also defined as type double
+! Potential density is additionally calculated in eqn_state_seawater,
+! first constant in array den also defined as type double
 !
 ! Former revisions:
@@ -69,5 +70,5 @@
        INTEGER ::  i, j, k
 
-       REAL ::  p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
+       REAL ::  pden, pnom, p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
 
        DO  i = nxl, nxr
@@ -76,6 +77,4 @@
 !
 !--             Pressure is needed in dbar
-!                p1 = hyp(0) * 1E-4
-!                p1 = 0.0
                 p1 = hyp(k) * 1E-4
                 p2 = p1 * p1
@@ -93,22 +92,35 @@
                 sa2  = sa1 * sa1
 
-                rho(k,j,i) =                                                   &
-                        ( nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
-                          nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
-                          nom(7)*sa2       + nom(8)*p1      + nom(9)*p1*pt2  + &
-                          nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2   &
-                        ) /                                                    &
-                        ( den(1)           + den(2)*pt1     + den(3)*pt2     + &
-                          den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
-                          den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
-                          den(10)*sa15*pt2 + den(11)*p1     + den(12)*p2*pt3 + &
-                          den(13)*p3*pt1                                       &
-                        )
+                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(k,j,i) = pnom / pden
 
              ENDDO
 !
 !--          Neumann conditions are assumed at bottom and top boundary
-             rho(nzt+1,j,i)            = rho(nzt,j,i)
-             rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
+             prho(nzt+1,j,i)            = prho(nzt,j,i)
+             prho(nzb_s_inner(j,i),j,i) = prho(nzb_s_inner(j,i)+1,j,i)
+             rho(nzt+1,j,i)             = rho(nzt,j,i)
+             rho(nzb_s_inner(j,i),j,i)  = rho(nzb_s_inner(j,i)+1,j,i)
+
           ENDDO
        ENDDO
@@ -129,11 +141,9 @@
        INTEGER ::  i, j, k
 
-       REAL ::  p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
+       REAL ::  pden, pnom, p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
 
        DO  k = nzb_s_inner(j,i)+1, nzt
 !
 !--       Pressure is needed in dbar
-!          p1 = hyp(0) * 1E-4
-!          p1 = 0.0
           p1 = hyp(k) * 1E-4
           p2 = p1 * p1
@@ -151,20 +161,35 @@
           sa2  = sa1 * sa1
 
-          rho(k,j,i) = ( nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
-                         nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
-                         nom(7)*sa2       + nom(8)*p1      + nom(9)*p1*pt2  + &
-                         nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2   &
-                       ) /                                                    &
-                       ( den(1)           + den(2)*pt1     + den(3)*pt2     + &
-                         den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
-                         den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
-                         den(10)*sa15*pt2 + den(11)*p1     + den(12)*p2*pt3 + &
-                         den(13)*p3*pt1                                       &
-                       )
+          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(k,j,i) = pnom / pden
+
+
        ENDDO
+
 !
 !--    Neumann conditions are assumed at bottom and top boundary
-       rho(nzt+1,j,i)            = rho(nzt,j,i)
-       rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
+       prho(nzt+1,j,i)            = prho(nzt,j,i)
+       prho(nzb_s_inner(j,i),j,i) = prho(nzb_s_inner(j,i)+1,j,i)
+       rho(nzt+1,j,i)             = rho(nzt,j,i)
+       rho(nzb_s_inner(j,i),j,i)  = rho(nzb_s_inner(j,i)+1,j,i)
 
     END SUBROUTINE eqn_state_seawater_ij
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 377)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 388)
@@ -4,4 +4,6 @@
 ! Current revisions:
 ! -----------------
+! Vertical profiles of potential density and hydrostatic pressure are
+! calculated.
 ! Added missing timeseries calculation of w"q"(0), moved timeseries q* to the
 ! end.
@@ -531,5 +533,5 @@
                 IF ( humidity )  THEN
                    sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
-                                       qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
+                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
                    IF ( cloud_physics )  THEN
                       sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
@@ -546,5 +548,5 @@
                 IF ( passive_scalar )  THEN
                    sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
-                                       qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
+                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
                 ENDIF
              ENDIF
@@ -597,4 +599,6 @@
                                                        rmask(j,i,sr)
                    sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
+                                                       rmask(j,i,sr)
+                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) * &
                                                        rmask(j,i,sr)
                 ENDIF
@@ -644,5 +648,8 @@
 !
 !--    Density at top follows Neumann condition
-       IF ( ocean )  sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
+       IF ( ocean )  THEN
+          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
+          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
+       ENDIF
 
 !
@@ -897,4 +904,6 @@
        hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
        hom(:,1,70,sr) = sums(:,70)     ! q*2
+       hom(:,1,71,sr) = sums(:,71)     ! prho
+       hom(:,1,72,sr) = hyp * 1E-4     ! hyp in kPa
 
        hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
@@ -916,5 +925,4 @@
 !--    is less than 1.5 times the height where the heat flux becomes negative
 !--    (positive) for the first time.
-!--    NOTE: This criterion is still capable of improving!
        z_i(1) = 0.0
        first = .TRUE.
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 377)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 388)
@@ -7,4 +7,5 @@
 ! Current revisions:
 ! -----------------
+! Initialization of prho added.
 ! bugfix: correction of initial volume flow for non-flat topography
 ! bugfix: zero initialization of arrays within buildings for 'cyclic_fill'
@@ -234,10 +235,12 @@
        ALLOCATE( saswsb_1(nys-1:nyn+1,nxl-1:nxr+1), &
                  saswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
-       ALLOCATE( rho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-                 sa_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
-                 sa_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
+       ALLOCATE( prho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
+                 rho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
+                 sa_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
+                 sa_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
                  sa_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
-       rho => rho_1  ! routine calc_mean_profile requires density to be a
-                     ! pointer
+       prho => prho_1
+       rho  => rho_1  ! routines calc_mean_profile and diffusion_e require
+                      ! density to be apointer
        IF ( humidity_remote )  THEN
           ALLOCATE( qswst_remote(nys-1:nyn+1,nxl-1:nxr+1) )
@@ -1328,4 +1331,5 @@
 !--    Initialize quantities needed for the ocean model
        CALL init_ocean
+
     ELSE
 !
Index: palm/trunk/SOURCE/init_ocean.f90
===================================================================
--- palm/trunk/SOURCE/init_ocean.f90	(revision 377)
+++ palm/trunk/SOURCE/init_ocean.f90	(revision 388)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! Bugfix: First calculation of hyp(0) changed
+! Bugfix: Initial profiles of hydrostatic pressure and density are calculated
+! iteratively. First calculation of hyp(0) changed.
 !
 ! Former revisions:
@@ -32,5 +33,5 @@
     INTEGER ::  k, n
 
-    REAL    ::  sa_l, pt_l, rho_l
+    REAL    ::  sa_l, pt_l
 
     REAL, DIMENSION(nzb:nzt+1) ::  rho_init
@@ -45,4 +46,5 @@
 !-- 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
 
@@ -55,12 +57,7 @@
     hyp(0) = hyp(1) + rho_surface * g * dzu(1)
 
-    IF ( myid == 0 )  THEN
-       print*,'hydro pres using rho_surface'
-       DO  k = nzt+1, 0, -1
-          print*, 'k = ', k, ' hyp = ', hyp(k)
-       ENDDO
-       print*, ' '
-    ENDIF
-
+!
+!-- Second step: Iteratively calculate in situ density (based on presssure)
+!-- and pressure (based on in situ density)
     DO  n = 1, 5
 
@@ -85,12 +82,4 @@
        ENDDO
 
-       IF ( myid == 0 )  THEN
-          print*,'hydro pres / rho  n = ', n
-          DO  k = nzt+1, 0, -1
-             print*, 'k = ', k, ' hyp = ', hyp(k), ' rho = ', rho_init(k)
-          ENDDO
-          print*, ' '
-       ENDIF
-
     ENDDO
 
@@ -111,6 +100,6 @@
 
 !
-!-- Calculate the initial potential density, based on the initial
-!-- temperature and salinity profile
+!-- Calculate the 3d array of initial in situ and potential density,
+!-- based on the initial temperature and salinity profile
     CALL eqn_state_seawater
 
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 377)
+++ palm/trunk/SOURCE/modules.f90	(revision 388)
@@ -5,4 +5,5 @@
 ! Current revisions:
 ! -----------------
+! +prho, prho_1
 ! +bc_lr_cyc, bc_ns_cyc
 ! +output_for_t0 
@@ -169,12 +170,12 @@
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
-          e_1, e_2, e_3, kh_1, kh_2, km_1, km_2, p, pt_1, pt_2, pt_3, q_1,     &
-          q_2, q_3, ql_1, ql_2, rho_1, sa_1, sa_2, sa_3, u_1, u_2, u_3, v_1,   &
-          v_2, v_3, vpt_1, vpt_2, w_1, w_2, w_3
+          e_1, e_2, e_3, kh_1, kh_2, km_1, km_2, p, prho_1, pt_1, pt_2, pt_3,  &
+          q_1, q_2, q_3, ql_1, ql_2, rho_1, sa_1, sa_2, sa_3, u_1, u_2, u_3,   &
+          v_1, v_2, v_3, vpt_1, vpt_2, w_1, w_2, w_3
 
     REAL, DIMENSION(:,:,:), POINTER ::                                         &
-          e, e_m, e_p, kh, kh_m, km, km_m, pt, pt_m, pt_p, q, q_m, q_p, ql,    &
-          ql_c, rho, sa, sa_p, te_m, tpt_m, tq_m, tsa_m, tu_m, tv_m, tw_m, u,  &
-          u_m, u_p, v, v_m, v_p, vpt, vpt_m, w, w_m, w_p
+          e, e_m, e_p, kh, kh_m, km, km_m, prho, pt, pt_m, pt_p, q, q_m, q_p,  &
+          ql, ql_c, rho, sa, sa_p, te_m, tpt_m, tq_m, tsa_m, tu_m, tv_m, tw_m, &
+          u, u_m, u_p, v, v_m, v_p, vpt, vpt_m, w, w_m, w_p
 
     REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  rif_wall
Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 377)
+++ palm/trunk/SOURCE/production_e.f90	(revision 388)
@@ -4,4 +4,6 @@
 ! Actual revisions:
 ! -----------------
+! Bugfix: wrong sign in buoyancy production of ocean part in case of not using
+!         the reference density (only in 3D routine production_e)
 ! Bugfix to avoid zero division by km_neutral
 !
@@ -450,6 +452,6 @@
                    DO  j = nys, nyn
                       DO  k = nzb_s_inner(j,i)+1, nzt
-                         tend(k,j,i) = tend(k,j,i) +                    &
-                                       kh(k,j,i) * g / prho_reference * &
+                         tend(k,j,i) = tend(k,j,i) +                   &
+                                       kh(k,j,i) * g / rho_reference * &
                                        ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
                       ENDDO
@@ -487,5 +489,5 @@
                    DO  j = nys, nyn
                       DO  k = nzb_s_inner(j,i)+1, nzt
-                         tend(k,j,i) = tend(k,j,i) -                &
+                         tend(k,j,i) = tend(k,j,i) +                &
                                        kh(k,j,i) * g / rho(k,j,i) * &
                                        ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
@@ -946,5 +948,5 @@
 !--             bottom and top surface layer
                 DO  k = nzb_s_inner(j,i)+1, nzt
-                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / prho_reference * &
+                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / rho_reference * &
                                       ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
                 ENDDO
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 377)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 388)
@@ -4,4 +4,5 @@
 ! Actual revisions:
 ! -----------------
+! prho is used instead of rho in diffusion_e,
 ! external pressure gradient
 !
@@ -330,5 +331,5 @@
           CALL coriolis( i, j, 3 )
           IF ( ocean )  THEN
-             CALL buoyancy( i, j, rho, prho_reference, 3, 64 )
+             CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
           ELSE
              IF ( .NOT. humidity )  THEN
@@ -725,10 +726,10 @@
                 IF ( .NOT. humidity )  THEN
                    IF ( ocean )  THEN
-                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,   &
-                                        l_grid, rho, prho_reference, rif, tend, &
-                                        zu, zw )
+                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,  &
+                                        l_grid, prho, prho_reference, rif,     &
+                                        tend, zu, zw )
                    ELSE
-                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
-                                        l_grid, pt, pt_reference, rif, tend,  &
+                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,  &
+                                        l_grid, pt, pt_reference, rif, tend,   &
                                         zu, zw )
                    ENDIF
@@ -769,5 +770,5 @@
                       IF ( ocean )  THEN
                          CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
-                                           km, l_grid, rho, prho_reference,   &
+                                           km, l_grid, prho, prho_reference,  &
                                            rif, tend, zu, zw )
                       ELSE
@@ -1013,5 +1014,5 @@
           CALL coriolis( i, j, 3 )
           IF ( ocean )  THEN
-             CALL buoyancy( i, j, rho, prho_reference, 3, 64 )
+             CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
           ELSE
              IF ( .NOT. humidity )  THEN
@@ -1263,6 +1264,6 @@
                 IF ( .NOT. humidity )  THEN
                    IF ( ocean )  THEN
-                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
-                                        km, l_grid, rho, prho_reference,  &
+                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
+                                        km, l_grid, prho, prho_reference,  &
                                         rif, tend, zu, zw )
                    ELSE
@@ -1552,5 +1553,5 @@
     CALL coriolis( 3 )
     IF ( ocean )  THEN
-       CALL buoyancy( rho, prho_reference, 3, 64 )
+       CALL buoyancy( rho, rho_reference, 3, 64 )
     ELSE
        IF ( .NOT. humidity )  THEN
@@ -1967,6 +1968,6 @@
           IF ( .NOT. humidity )  THEN
              IF ( ocean )  THEN
-                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, rho, &
-                                  prho_reference, rif, tend, zu, zw )
+                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
+                                  prho, prho_reference, rif, tend, zu, zw )
              ELSE
                 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
@@ -2004,5 +2005,5 @@
                 IF ( ocean )  THEN
                    CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
-                                     rho, prho_reference, rif, tend, zu, zw )
+                                     prho, prho_reference, rif, tend, zu, zw )
                 ELSE
                    CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 377)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 388)
@@ -4,4 +4,5 @@
 ! Actual revisions:
 ! -----------------
+! Using prho instead of rho in diffusvities.
 ! Coupling with independent precursor runs.
 ! Bugfix: output of particle time series only if particle advection is switched
@@ -284,5 +285,5 @@
              IF ( .NOT. humidity ) THEN
                 IF ( ocean )  THEN
-                   CALL diffusivities( rho, prho_reference )
+                   CALL diffusivities( prho, prho_reference )
                 ELSE
                    CALL diffusivities( pt, pt_reference )
