Index: palm/trunk/SOURCE/buoyancy.f90
===================================================================
--- palm/trunk/SOURCE/buoyancy.f90	(revision 1172)
+++ 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
 
