Index: palm/trunk/SOURCE/buoyancy.f90
===================================================================
--- palm/trunk/SOURCE/buoyancy.f90	(revision 1178)
+++ palm/trunk/SOURCE/buoyancy.f90	(revision 1179)
@@ -19,6 +19,7 @@
 !
 ! Currrent revisions:
-! -----------------
-!
+! ------------------
+! steering of reference state revised (var_reference and pr removed from
+! parameter list), use_reference renamed use_single_reference_value
 !
 ! Former revisions:
@@ -108,5 +109,5 @@
 ! Call for all grid points
 !------------------------------------------------------------------------------!
-    SUBROUTINE buoyancy( var, var_reference, wind_component, pr )
+    SUBROUTINE buoyancy( var, wind_component )
 
        USE arrays_3d
@@ -114,10 +115,8 @@
        USE indices
        USE pegrid
-       USE statistics
 
        IMPLICIT NONE
 
-       INTEGER ::  i, j, k, pr, wind_component
-       REAL    ::  var_reference
+       INTEGER ::  i, j, k, wind_component
 #if defined( __nopointer )
        REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
@@ -130,29 +129,14 @@
 !
 !--       Normal case: horizontal surface
-          IF ( use_reference )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_s_inner(j,i)+1, nzt-1
-                      tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
-                                                            (                  &
-                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / var_reference + &
-                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference   &
-                                                            )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_s_inner(j,i)+1, nzt-1
-                      tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
-                                                            (                  &
-                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
-                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
-                                                            )
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb_s_inner(j,i)+1, nzt-1
+                   tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (  &
+                          ( var(k,j,i)   - ref_state(k) )   / ref_state(k) +   &
+                          ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1)   &
+                                                                            )
+                ENDDO
+             ENDDO
+          ENDDO
 
        ELSE
@@ -205,5 +189,5 @@
 ! Call for all grid points - accelerator version
 !------------------------------------------------------------------------------!
-    SUBROUTINE buoyancy_acc( var, var_reference, wind_component, pr )
+    SUBROUTINE buoyancy_acc( var, wind_component )
 
        USE arrays_3d
@@ -211,10 +195,8 @@
        USE indices
        USE pegrid
-       USE statistics
 
        IMPLICIT NONE
 
-       INTEGER ::  i, j, k, pr, wind_component
-       REAL    ::  var_reference
+       INTEGER ::  i, j, k, wind_component
 #if defined( __nopointer )
        REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
@@ -227,37 +209,18 @@
 !
 !--       Normal case: horizontal surface
-          IF ( use_reference )  THEN
-             !$acc kernels present( nzb_s_inner, hom, tend, var )
-             !$acc loop
-             DO  i = i_left, i_right
-                DO  j = j_south, j_north
-                   !$acc loop independent vector(32)
-                   DO  k = nzb_s_inner(j,i)+1, nzt-1
-                      tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
-                                                            (                  &
-                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / var_reference + &
-                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference   &
-                                                            )
-                   ENDDO
-                ENDDO
-             ENDDO
-             !$acc end kernels
-          ELSE
-             !$acc kernels present( nzb_s_inner, hom, tend, var )
-             !$acc loop
-             DO  i = i_left, i_right
-                DO  j = j_south, j_north
-                   !$acc loop independent vector(32)
-                   DO  k = nzb_s_inner(j,i)+1, nzt-1
-                      tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
-                                                            (                  &
-                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
-                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
-                                                            )
-                   ENDDO
-                ENDDO
-             ENDDO
-             !$acc end kernels
-          ENDIF
+          !$acc kernels present( nzb_s_inner, ref_state, tend, var )
+          !$acc loop
+          DO  i = i_left, i_right
+             DO  j = j_south, j_north
+                !$acc loop independent vector(32)
+                DO  k = nzb_s_inner(j,i)+1, nzt-1
+                   tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (  &
+                          ( var(k,j,i)   - ref_state(k) )   / ref_state(k) +   &
+                          ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1)   &
+                                                                            )
+                ENDDO
+             ENDDO
+          ENDDO
+          !$acc end kernels
 
        ELSE
@@ -313,5 +276,5 @@
 !------------------------------------------------------------------------------!
 !pgi$r opt=1
-    SUBROUTINE buoyancy_ij( i, j, var, var_reference, wind_component, pr )
+    SUBROUTINE buoyancy_ij( i, j, var, wind_component )
 
        USE arrays_3d
@@ -319,10 +282,8 @@
        USE indices
        USE pegrid
-       USE statistics
 
        IMPLICIT NONE
 
        INTEGER ::  i, j, k, pr, wind_component
-       REAL    ::  var_reference
 #if defined( __nopointer )
        REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
@@ -335,19 +296,10 @@
 !
 !--       Normal case: horizontal surface
-          IF ( use_reference )  THEN
-             DO  k = nzb_s_inner(j,i)+1, nzt-1
-                 tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (   &
-                         ( var(k,j,i)   - hom(k,1,pr,0)   ) / var_reference + &
-                         ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference   &
-                                                                          )
-             ENDDO
-          ELSE
-             DO  k = nzb_s_inner(j,i)+1, nzt-1
-                 tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (    &
-                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
-                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
-                                                                          )
-             ENDDO
-          ENDIF
+          DO  k = nzb_s_inner(j,i)+1, nzt-1
+              tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (     &
+                        ( var(k,j,i)   - ref_state(k)   ) / ref_state(k)   + &
+                        ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1)   &
+                                                                       )
+          ENDDO
 
        ELSE
@@ -395,7 +347,9 @@
 ! ------------
 ! Calculate the horizontally averaged vertical temperature profile (pr=4 in case
-! of potential temperature and 44 in case of virtual potential temperature).
-!------------------------------------------------------------------------------!
-
+! of potential temperature, 44 in case of virtual potential temperature, and 64
+! in case of density (ocean runs)).
+!------------------------------------------------------------------------------!
+
+       USE arrays_3d,  ONLY: ref_state
        USE control_parameters
        USE indices
@@ -457,4 +411,6 @@
        ENDIF
 
+       ref_state(:)  = hom(:,1,pr,0)   ! this is used in the buoyancy term
+
     END SUBROUTINE calc_mean_profile
 
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 1178)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 1179)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! checks and settings of buoyancy parameters and switches revised,
+! initial profile for rho added to hom (id=77)
 !
 ! Former revisions:
@@ -1323,7 +1324,5 @@
        large_scale_subsidence = .TRUE.
        CALL init_w_subsidence
-    END IF
- 
-             
+    ENDIF
 
 !
@@ -1333,11 +1332,25 @@
 
 !
-!-- Ocean runs always use reference values in the buoyancy term. Therefore
-!-- set the reference temperature equal to the surface temperature.
-    IF ( ocean  .AND.  pt_reference == 9999999.9 )  pt_reference = pt_surface
-
-!
-!-- Reference value has to be used in buoyancy terms
-    IF ( pt_reference /= 9999999.9 )  use_reference = .TRUE.
+!-- Check and set buoyancy related parameters and switches
+    IF ( reference_state == 'horizontal_average' )  THEN
+       CONTINUE
+    ELSEIF ( reference_state == 'initial_profile' )  THEN
+       use_initial_profile_as_reference = .TRUE.
+    ELSEIF ( reference_state == 'single_value' )  THEN
+       use_single_reference_value = .TRUE.
+       IF ( pt_reference == 9999999.9 )  pt_reference = pt_surface
+       vpt_reference = pt_reference * ( 1.0 + 0.61 * q_surface )
+    ELSE
+       message_string = 'illegal value for reference_state: "' // &
+                        TRIM( reference_state ) // '"'
+       CALL message( 'check_parameters', 'PA0056', 1, 2, 0, 6, 0 )
+    ENDIF
+
+!
+!-- Ocean runs always use reference values in the buoyancy term
+    IF ( ocean )  THEN
+       reference_state = 'single_value'
+       use_single_reference_value = .TRUE.
+    ENDIF
 
 !
@@ -2463,4 +2476,10 @@
                 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    ! because zu(nzb) is negative
+                   data_output_pr(i)     = data_output_pr(i)(2:)
+                ENDIF
              ENDIF
 
Index: palm/trunk/SOURCE/diffusion_e.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_e.f90	(revision 1178)
+++ palm/trunk/SOURCE/diffusion_e.f90	(revision 1179)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! use_reference renamed use_single_reference_value
 !
 ! Former revisions:
@@ -134,5 +134,5 @@
 !--    This if clause must be outside the k-loop because otherwise
 !--    runtime errors occur with -C hopt on NEC
-       IF ( use_reference )  THEN
+       IF ( use_single_reference_value )  THEN
 
           DO  i = nxl, nxr
@@ -319,5 +319,5 @@
 !--    This if clause must be outside the k-loop because otherwise
 !--    runtime errors occur with -C hopt on NEC
-       IF ( use_reference )  THEN
+       IF ( use_single_reference_value )  THEN
 
           !$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
@@ -506,5 +506,5 @@
                     ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
           IF ( dvar_dz > 0.0 ) THEN
-             IF ( use_reference )  THEN
+             IF ( use_single_reference_value )  THEN
                 l_stable = 0.76 * SQRT( e(k,j,i) ) / &
                                   SQRT( g / var_reference * dvar_dz ) + 1E-5
Index: palm/trunk/SOURCE/diffusivities.f90
===================================================================
--- palm/trunk/SOURCE/diffusivities.f90	(revision 1178)
+++ palm/trunk/SOURCE/diffusivities.f90	(revision 1179)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! use_reference renamed use_single_reference_value
 !
 ! Former revisions:
@@ -132,5 +132,5 @@
                           ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
                 IF ( dvar_dz > 0.0 ) THEN
-                   IF ( use_reference )  THEN
+                   IF ( use_single_reference_value )  THEN
                       l_stable = 0.76 * sqrt_e / &
                                  SQRT( g / var_reference * dvar_dz ) + 1E-5
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 1178)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 1179)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! comment for profile 77 added
 !
 ! Former revisions:
@@ -1148,4 +1148,5 @@
        hom(:,1,75,sr) = sums(:,75)     ! qc
        hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
+                                       ! 77 is initial density profile
 
        hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 1178)
+++ palm/trunk/SOURCE/header.f90	(revision 1179)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! output of reference_state, use_reference renamed use_single_reference_value
 !
 ! Former revisions:
@@ -1230,9 +1230,10 @@
 !-- Other quantities
     WRITE ( io, 411 )  g
-    IF ( use_reference )  THEN
+    WRITE ( io, 412 )  TRIM( reference_state )
+    IF ( use_single_reference_value )  THEN
        IF ( ocean )  THEN
-          WRITE ( io, 412 )  prho_reference
-       ELSE
-          WRITE ( io, 413 )  pt_reference
+          WRITE ( io, 413 )  prho_reference
+       ELSE
+          WRITE ( io, 414 )  pt_reference
        ENDIF
     ENDIF
@@ -1900,6 +1901,7 @@
             '                            f*    = ',F9.6,' 1/s')
 411 FORMAT (/'    Gravity             :   g     = ',F4.1,' m/s**2')
-412 FORMAT (/'    Reference density in buoyancy terms: ',F8.3,' kg/m**3')
-413 FORMAT (/'    Reference temperature in buoyancy terms: ',F8.4,' K')
+412 FORMAT (/'    Reference state used in buoyancy terms: ',A)
+413 FORMAT ('       Reference density in buoyancy terms: ',F8.3,' kg/m**3')
+414 FORMAT ('       Reference temperature in buoyancy terms: ',F8.4,' K')
 415 FORMAT (/'    Cloud physics parameters:'/ &
              '    ------------------------'/)
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 1178)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 1179)
@@ -23,5 +23,5 @@
 ! Current revisions:
 ! ------------------
-!
+! allocate and set ref_state to be used in buoyancy terms
 !
 ! Former revisions:
@@ -279,5 +279,6 @@
               sums_divnew_l(0:statistic_regions),                           &
               sums_divold_l(0:statistic_regions) )
-    ALLOCATE( dp_smooth_factor(nzb:nzt), rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt) )
+    ALLOCATE( dp_smooth_factor(nzb:nzt), rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt), &
+              ref_state(nzb:nzt+1) )
     ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                 &
               ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),               &
@@ -1039,4 +1040,20 @@
        ENDIF
 
+!
+!--    Set the reference state to be used in the buoyancy terms (for ocean runs
+!--    the reference state will be set (overwritten) in init_ocean)
+       IF ( use_single_reference_value )  THEN
+          IF ( .NOT. humidity )  THEN
+             ref_state(:) = pt_reference
+          ELSE
+             ref_state(:) = vpt_reference
+          ENDIF
+       ELSE
+          IF ( .NOT. humidity )  THEN
+             ref_state(:) = pt_init(:)
+          ELSE
+             ref_state(:) = vpt(:,nys,nxl)
+          ENDIF
+       ENDIF
 
 !
Index: palm/trunk/SOURCE/init_ocean.f90
===================================================================
--- palm/trunk/SOURCE/init_ocean.f90	(revision 1178)
+++ palm/trunk/SOURCE/init_ocean.f90	(revision 1179)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Initial density profile is stored in array hom
 !
 ! Former revisions:
@@ -47,7 +47,8 @@
     USE control_parameters
     USE eqn_state_seawater_mod
-    USE pegrid
     USE grid_variables
     USE indices
+    USE pegrid
+    USE statistics
 
     IMPLICIT NONE
@@ -126,4 +127,16 @@
     CALL eqn_state_seawater
 
+!
+!-- Store initial density profile
+    hom(:,1,77,:)  = SPREAD( rho_init(:), 2, statistic_regions+1 )
+
+!
+!-- Set the reference state to be used in the buoyancy terms
+    IF ( use_single_reference_value )  THEN
+       ref_state(:) = prho_reference
+    ELSE
+       ref_state(:) = rho_init(:)
+    ENDIF
+
 
  END SUBROUTINE init_ocean
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 1178)
+++ palm/trunk/SOURCE/modules.f90	(revision 1179)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-!
+! +reference_state, ref_state, use_initial_profile_as_reference, vpt_reference,
+! use_reference renamed use_single_reference_value
 !
 ! Former revisions:
@@ -420,6 +421,6 @@
           dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor, lad, l_grid,      &
           nc_1d, nr_1d, ptdf_x, ptdf_y, pt_1d, pt_init, q_1d, q_init, qc_1d,   &
-          qr_1d, rdf, rdf_sc, sa_init, ug, u_init, u_nzb_p1_for_vfc, vg,       &
-          v_init, v_nzb_p1_for_vfc, w_subs, zu, zw
+          qr_1d, rdf, rdf_sc, ref_state, sa_init, ug, u_init,                  &
+          u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, w_subs, zu, zw
 
     REAL, DIMENSION(:,:), ALLOCATABLE ::                                       &
@@ -656,4 +657,5 @@
                              mixing_length_1d = 'as_in_3d_model', &
                              random_generator = 'numerical-recipes', &
+                             reference_state = 'initial_profile', &
                              return_addres, return_username, &
                              timestep_scheme = 'runge-kutta-3'                             
@@ -770,5 +772,7 @@
                 terminate_run = .FALSE., turbulence = .FALSE., &
                 turbulent_inflow = .FALSE., use_cmax = .TRUE., &
-                use_prescribed_profile_data = .FALSE., use_reference = .FALSE.,&
+                use_initial_profile_as_reference = .FALSE., &
+                use_prescribed_profile_data = .FALSE., &
+                use_single_reference_value = .FALSE., &
                 use_surface_fluxes = .FALSE., use_top_fluxes = .FALSE., &
                 use_ug_for_galilei_tr = .TRUE., use_upstream_for_tke = .FALSE.,&
@@ -846,5 +850,5 @@
              top_momentumflux_v = 9999999.9, top_salinityflux = 9999999.9, &
              ug_surface = 0.0, u_bulk = 0.0, u_gtrans = 0.0, &
-             vg_surface = 0.0, &
+             vg_surface = 0.0, vpt_reference = 9999999.9, &
              v_bulk = 0.0, v_gtrans = 0.0, wall_adjustment_factor = 1.8, &
              z_max_do2d = -1.0, z0h_factor = 1.0
Index: palm/trunk/SOURCE/palm.f90
===================================================================
--- palm/trunk/SOURCE/palm.f90	(revision 1178)
+++ palm/trunk/SOURCE/palm.f90	(revision 1179)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! ref_state added to copyin-list
 !
 ! Former revisions:
@@ -241,5 +241,5 @@
 !-- host values
     !$acc  data copyin( d, diss, e, e_p, kh, km, p, pt, pt_p, q, ql, tend, te_m, tpt_m, tu_m, tv_m, tw_m, u, u_p, v, vpt, v_p, w, w_p )          &
-    !$acc       copyin( tric, dzu, ddzu, ddzw, dd2zu, l_grid, l_wall, ptdf_x, ptdf_y, pt_init, rdf, rdf_sc, ug, u_init, vg, v_init, zu, zw )   &
+    !$acc       copyin( tric, dzu, ddzu, ddzw, dd2zu, l_grid, l_wall, ptdf_x, ptdf_y, pt_init, rdf, rdf_sc, ref_state, ug, u_init, vg, v_init, zu, zw )   &
     !$acc       copyin( hom, qs, qsws, qswst, rif, rif_wall, shf, ts, tswst, us, usws, uswst, vsws, vswst, z0, z0h )      &
     !$acc       copyin( fxm, fxp, fym, fyp, fwxm, fwxp, fwym, fwyp, nzb_diff_s_inner, nzb_diff_s_outer, nzb_diff_u )       &
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 1178)
+++ palm/trunk/SOURCE/parin.f90	(revision 1179)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! +reference_state in inipar
 !
 ! Former revisions:
@@ -234,6 +234,5 @@
              mg_switch_to_pe0_level, mixing_length_1d, momentum_advec, &
              nc_const, netcdf_precision, neutral, ngsrb, &
-             nsor, &
-             nsor_ini, nx, ny, nz, ocean, omega, omega_sor, &
+             nsor, nsor_ini, nx, ny, nz, ocean, omega, omega_sor, &
              passive_scalar, pch_index, phi, plant_canopy, prandtl_layer, &
              prandtl_number, precipitation, psolver, pt_damping_factor, &
@@ -242,7 +241,7 @@
              pt_vertical_gradient_level, q_surface, q_surface_initial_change, &
              q_vertical_gradient, q_vertical_gradient_level, &
-             radiation, &
-             random_generator, random_heatflux, rayleigh_damping_factor, &
-             rayleigh_damping_height, recycling_width, residual_limit, &
+             radiation, random_generator, random_heatflux, &
+             rayleigh_damping_factor, rayleigh_damping_height, recycling_width,&
+             reference_state, residual_limit, &
              rif_max, rif_min, roughness_length, sa_surface, &
              sa_vertical_gradient, sa_vertical_gradient_level, scalar_advec, &
Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 1178)
+++ palm/trunk/SOURCE/production_e.f90	(revision 1179)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! use_reference renamed use_single_reference_value
 !
 ! Former revisions:
@@ -497,5 +497,5 @@
              IF ( .NOT. humidity )  THEN
 
-                IF ( use_reference )  THEN
+                IF ( use_single_reference_value )  THEN
 
                    IF ( ocean )  THEN
@@ -1099,5 +1099,5 @@
           IF ( .NOT. humidity )  THEN
 
-             IF ( use_reference )  THEN
+             IF ( use_single_reference_value )  THEN
 
                 IF ( ocean )  THEN
@@ -1677,5 +1677,5 @@
           IF ( .NOT. humidity )  THEN
 
-             IF ( use_reference )  THEN
+             IF ( use_single_reference_value )  THEN
 
                 IF ( ocean )  THEN
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 1178)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 1179)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! two arguments removed from routine buoyancy, ref_state updated on device
 !
 ! Former revisions:
@@ -267,5 +267,5 @@
              CALL coriolis( i, j, 1 )
              IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
-                CALL buoyancy( i, j, pt, pt_reference, 1, 4 )
+                CALL buoyancy( i, j, pt, 1 )
              ENDIF
 
@@ -380,10 +380,10 @@
           IF ( .NOT. neutral )  THEN
              IF ( ocean )  THEN
-                CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
+                CALL buoyancy( i, j, rho, 3 )
              ELSE
                 IF ( .NOT. humidity )  THEN
-                   CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
+                   CALL buoyancy( i, j, pt, 3 )
                 ELSE
-                   CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
+                   CALL buoyancy( i, j, vpt, 3 )
                 ENDIF
              ENDIF
@@ -825,5 +825,5 @@
     CALL coriolis( 1 )
     IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
-       CALL buoyancy( pt, pt_reference, 1, 4 )
+       CALL buoyancy( pt, 1 )
     ENDIF
 
@@ -974,10 +974,10 @@
     IF ( .NOT. neutral )  THEN
        IF ( ocean )  THEN
-          CALL buoyancy( rho, rho_reference, 3, 64 )
+          CALL buoyancy( rho, 3 )
        ELSE
           IF ( .NOT. humidity )  THEN
-             CALL buoyancy( pt, pt_reference, 3, 4 )
+             CALL buoyancy( pt, 3 )
           ELSE
-             CALL buoyancy( vpt, pt_reference, 3, 44 )
+             CALL buoyancy( vpt, 3 )
           ENDIF
        ENDIF
@@ -1450,5 +1450,5 @@
 !-- u-velocity component
 !++ Statistics still not ported to accelerators
-    !$acc update device( hom )
+    !$acc update device( hom, ref_state )
     CALL cpu_log( log_point(5), 'u-equation', 'start' )
 
@@ -1466,5 +1466,5 @@
     CALL coriolis_acc( 1 )
     IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
-       CALL buoyancy( pt, pt_reference, 1, 4 )
+       CALL buoyancy( pt, 1 )
     ENDIF
 
@@ -1595,10 +1595,10 @@
     IF ( .NOT. neutral )  THEN
        IF ( ocean )  THEN
-          CALL buoyancy( rho, rho_reference, 3, 64 )
+          CALL buoyancy( rho, 3 )
        ELSE
           IF ( .NOT. humidity )  THEN
-             CALL buoyancy_acc( pt, pt_reference, 3, 4 )
+             CALL buoyancy_acc( pt, 3 )
           ELSE
-             CALL buoyancy( vpt, pt_reference, 3, 44 )
+             CALL buoyancy( vpt, 3 )
           ENDIF
        ENDIF
Index: palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- palm/trunk/SOURCE/read_var_list.f90	(revision 1178)
+++ palm/trunk/SOURCE/read_var_list.f90	(revision 1179)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! +ref_state
 !
 ! Former revisions:
@@ -188,5 +188,5 @@
 !-- Make version number check first
     READ ( 13 )  version_on_file
-    binary_version = '3.9'
+    binary_version = '3.9a'
     IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
        WRITE( message_string, * ) 'version mismatch concerning control ', &
@@ -551,4 +551,8 @@
           CASE ( 'recycling_width' )
              READ ( 13 )  recycling_width
+          CASE ( 'reference_state' )
+             READ ( 13 )  reference_state
+          CASE ( 'ref_state' )
+             READ ( 13 )  ref_state
           CASE ( 'residual_limit' )
              READ ( 13 )  residual_limit
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 1178)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 1179)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-!
+! mean profiles for reference state are only calculated if required,
+! small bugfix for background communication
 !
 ! Former revisions:
@@ -230,7 +231,15 @@
 !--       Calculate those variables needed in the tendency terms which need
 !--       global communication
-          IF ( .NOT. neutral )  CALL calc_mean_profile( pt, 4 )
-          IF ( ocean         )  CALL calc_mean_profile( rho, 64 )
-          IF ( humidity      )  CALL calc_mean_profile( vpt, 44 )
+          IF ( .NOT. use_single_reference_value  .AND. &
+               .NOT. use_initial_profile_as_reference )  THEN
+!
+!--          Horizontally averaged profiles to be used as reference state in
+!--          buoyancy terms (WARNING: only the respective last call of
+!--          calc_mean_profile defines the reference state!)
+             IF ( .NOT. neutral )  CALL calc_mean_profile( pt, 4 )
+             IF ( ocean         )  CALL calc_mean_profile( rho, 64 )
+             IF ( humidity      )  CALL calc_mean_profile( vpt, 44 )
+          ENDIF
+
           IF ( .NOT. constant_diffusion )  CALL production_e_init
           IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
@@ -438,5 +447,6 @@
              IF (humidity  .OR.  passive_scalar)  THEN
                 CALL exchange_horiz( q_p, nbgp )
-                IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                IF ( cloud_physics .AND. icloud_scheme == 0  .AND.  &
+                     precipitation )  THEN
                    CALL exchange_horiz( qr_p, nbgp )
                    CALL exchange_horiz( nr_p, nbgp )
Index: palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- palm/trunk/SOURCE/write_var_list.f90	(revision 1178)
+++ palm/trunk/SOURCE/write_var_list.f90	(revision 1179)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! +reference_state, ref_state
 !
 ! Former revisions:
@@ -169,5 +169,5 @@
 
 
-    binary_version = '3.9'
+    binary_version = '3.9a'
 
     WRITE ( 14 )  binary_version
@@ -468,4 +468,8 @@
     WRITE ( 14 )  'recycling_width               '
     WRITE ( 14 )  recycling_width
+    WRITE ( 14 )  'reference_state               '
+    WRITE ( 14 )  reference_state
+    WRITE ( 14 )  'ref_state                     '
+    WRITE ( 14 )  ref_state
     WRITE ( 14 )  'residual_limit                '
     WRITE ( 14 )  residual_limit
