Index: palm/trunk/SOURCE/buoyancy.f90
===================================================================
--- palm/trunk/SOURCE/buoyancy.f90	(revision 96)
+++ palm/trunk/SOURCE/buoyancy.f90	(revision 97)
@@ -4,4 +4,7 @@
 ! Actual revisions:
 ! -----------------
+! Routine reneralized to be used with temperature AND density:
+! argument theta renamed var, new argument var_reference,
+! use_pt_reference renamed use_reference,
 ! calc_mean_pt_profile renamed calc_mean_profile
 !
@@ -46,5 +49,5 @@
 ! Call for all grid points
 !------------------------------------------------------------------------------!
-    SUBROUTINE buoyancy( theta, wind_component, pr )
+    SUBROUTINE buoyancy( var, var_reference, wind_component, pr )
 
        USE arrays_3d
@@ -57,5 +60,6 @@
 
        INTEGER ::  i, j, k, pr, wind_component
-       REAL, DIMENSION(:,:,:), POINTER  ::  theta
+       REAL    ::  var_reference
+       REAL, DIMENSION(:,:,:), POINTER  ::  var
 
 
@@ -63,11 +67,12 @@
 !
 !--       Normal case: horizontal surface
-          IF ( use_pt_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) + g * 0.5 * (                 &
-                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / pt_reference + &
-                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / pt_reference   &
+          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
@@ -78,7 +83,8 @@
                 DO  j = nys, nyn
                    DO  k = nzb_s_inner(j,i)+1, nzt-1
-                      tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                  &
-                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
-                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
+                      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
@@ -136,5 +142,5 @@
 ! Call for grid point i,j
 !------------------------------------------------------------------------------!
-    SUBROUTINE buoyancy_ij( i, j, theta, wind_component, pr )
+    SUBROUTINE buoyancy_ij( i, j, var, var_reference, wind_component, pr )
 
        USE arrays_3d
@@ -147,5 +153,6 @@
 
        INTEGER ::  i, j, k, pr, wind_component
-       REAL, DIMENSION(:,:,:), POINTER  ::  theta
+       REAL    ::  var_reference
+       REAL, DIMENSION(:,:,:), POINTER  ::  var
 
 
@@ -153,17 +160,17 @@
 !
 !--       Normal case: horizontal surface
-          IF ( use_pt_reference )  THEN
-             DO  k = nzb_s_inner(j,i)+1, nzt-1
-                 tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                      &
-                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / pt_reference + &
-                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / pt_reference   &
-                                                       )
-             ENDDO
-          ELSE
-             DO  k = nzb_s_inner(j,i)+1, nzt-1
-                 tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                       &
-                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
-                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
-                                                       )
+          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
