Index: palm/trunk/SOURCE/average_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/average_3d_data.f90	(revision 2030)
+++ palm/trunk/SOURCE/average_3d_data.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! renamed variable rho to rho_ocean and rho_av to rho_ocean_av
 ! 
 ! Former revisions:
@@ -318,9 +318,9 @@
              ENDDO
 
-          CASE ( 'rho' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
-                      rho_av(k,j,i) = rho_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+          CASE ( 'rho_ocean' )
+             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
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 2030)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! renamed variable rho to rho_ocean
 ! 
 ! Former revisions:
@@ -299,5 +299,5 @@
 ! 1179 2013-06-14 05:57:58Z raasch
 ! checks and settings of buoyancy parameters and switches revised,
-! initial profile for rho added to hom (id=77)
+! initial profile for rho_ocean added to hom (id=77)
 !
 ! 1174 2013-05-31 10:28:08Z gryschka
@@ -2528,5 +2528,5 @@
              hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
 
-          CASE ( 'rho' )
+          CASE ( 'rho_ocean' )
              IF (  .NOT.  ocean ) THEN
                 message_string = 'data_output_pr = ' //                        &
@@ -3005,5 +3005,5 @@
              unit = 'kg/kg'
 
-          CASE ( 'rho' )
+          CASE ( 'rho_ocean' )
              IF (  .NOT.  ocean )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
Index: palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_2d.f90	(revision 2030)
+++ palm/trunk/SOURCE/data_output_2d.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! renamed variable rho to rho_ocean and rho_av to rho_ocean_av
 ! 
 ! Former revisions:
@@ -157,5 +157,5 @@
     USE arrays_3d,                                                             &
         ONLY:  dzw, e, nr, ol, p, pt, precipitation_amount, precipitation_rate,&
-               prr,q, qc, ql, ql_c, ql_v, ql_vp, qr, qsws, rho, s, sa, shf,    &
+               prr,q, qc, ql, ql_c, ql_v, ql_vp, qr, qsws, rho_ocean, s, sa, shf,    &
                ssws, tend, ts, u, us, v, vpt, w, z0, z0h, z0q, zu, zw
         
@@ -769,7 +769,7 @@
              CASE ( 'rho_xy', 'rho_xz', 'rho_yz' )
                 IF ( av == 0 )  THEN
-                   to_be_resorted => rho
-                ELSE
-                   to_be_resorted => rho_av
+                   to_be_resorted => rho_ocean
+                ELSE
+                   to_be_resorted => rho_ocean_av
                 ENDIF
 
Index: palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_3d.f90	(revision 2030)
+++ palm/trunk/SOURCE/data_output_3d.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! renamed variable rho to rho_ocean and rho_av to rho_ocean_av
 ! 
 ! Former revisions:
@@ -146,5 +146,5 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  e, nr, p, pt, prr, q, qc, ql, ql_c, ql_v, qr, rho, s, sa, tend, &
+        ONLY:  e, nr, p, pt, prr, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, sa, tend, &
                u, v, vpt, w
         
@@ -531,9 +531,9 @@
              ENDIF
 
-          CASE ( 'rho' )
-             IF ( av == 0 )  THEN
-                to_be_resorted => rho
-             ELSE
-                to_be_resorted => rho_av
+          CASE ( 'rho_ocean' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rho_ocean
+             ELSE
+                to_be_resorted => rho_ocean_av
              ENDIF
 
Index: palm/trunk/SOURCE/data_output_mask.f90
===================================================================
--- palm/trunk/SOURCE/data_output_mask.f90	(revision 2030)
+++ palm/trunk/SOURCE/data_output_mask.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! renamed variable rho to rho_ocean and rho_av to rho_ocean_av
 ! 
 ! Former revisions:
@@ -102,10 +102,10 @@
 #if defined( __netcdf )
     USE arrays_3d,                                                             &
-        ONLY:  e, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho, s, sa, tend, u,   &
+        ONLY:  e, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, sa, tend, u,   &
                v, vpt, w
     
     USE averaging,                                                             &
         ONLY:  e_av, lpt_av, nr_av, p_av, pc_av, pr_av, pt_av, q_av, qc_av,    &
-               ql_av, ql_c_av, ql_v_av, ql_vp_av, qv_av, qr_av, rho_av, s_av,  &
+               ql_av, ql_c_av, ql_v_av, ql_vp_av, qv_av, qr_av, rho_ocean_av, s_av,  &
                sa_av, u_av, v_av, vpt_av, w_av 
     
@@ -436,9 +436,9 @@
              ENDIF
 
-          CASE ( 'rho' )
-             IF ( av == 0 )  THEN
-                to_be_resorted => rho
-             ELSE
-                to_be_resorted => rho_av
+          CASE ( 'rho_ocean' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rho_ocean
+             ELSE
+                to_be_resorted => rho_ocean_av
              ENDIF
 
Index: palm/trunk/SOURCE/diffusivities.f90
===================================================================
--- palm/trunk/SOURCE/diffusivities.f90	(revision 2030)
+++ palm/trunk/SOURCE/diffusivities.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! renamed variable rho to rho_ocean
 ! 
 ! Former revisions:
@@ -160,5 +160,5 @@
 !
 !--             Determine the mixing length
-                dvar_dz = atmos_ocean_sign * &  ! inverse effect of pt/rho gradient
+                dvar_dz = atmos_ocean_sign * &  ! inverse effect of pt/rho_ocean gradient
                           ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
                 IF ( dvar_dz > 0.0_wp ) THEN
Index: palm/trunk/SOURCE/eqn_state_seawater.f90
===================================================================
--- palm/trunk/SOURCE/eqn_state_seawater.f90	(revision 2030)
+++ palm/trunk/SOURCE/eqn_state_seawater.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! renamed variable rho to rho_ocean
 ! 
 ! Former revisions:
@@ -115,5 +115,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  hyp, prho, pt_p, rho, sa_p
+           ONLY:  hyp, prho, pt_p, rho_ocean, sa_p
        USE indices,                                                            &
            ONLY:  nxl, nxr, nyn, nys, nzb_s_inner, nzt
@@ -181,5 +181,5 @@
 !
 !--             In-situ density
-                rho(k,j,i) = pnom / pden
+                rho_ocean(k,j,i) = pnom / pden
 
              ENDDO
@@ -188,6 +188,6 @@
              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)
+             rho_ocean(nzt+1,j,i)             = rho_ocean(nzt,j,i)
+             rho_ocean(nzb_s_inner(j,i),j,i)  = rho_ocean(nzb_s_inner(j,i)+1,j,i)
 
           ENDDO
@@ -205,5 +205,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  hyp, prho, pt_p, rho, sa_p
+           ONLY:  hyp, prho, pt_p, rho_ocean, sa_p
            
        USE indices,                                                            &
@@ -255,5 +255,5 @@
 !
 !--       In-situ density
-          rho(k,j,i) = pnom / pden
+          rho_ocean(k,j,i) = pnom / pden
 
 
@@ -264,6 +264,6 @@
        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)
+       rho_ocean(nzt+1,j,i)             = rho_ocean(nzt,j,i)
+       rho_ocean(nzb_s_inner(j,i),j,i)  = rho_ocean(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 2030)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! renamed variable rho to rho_ocean
 ! 
 ! Former revisions:
@@ -212,5 +212,5 @@
     USE arrays_3d,                                                             &
         ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, ol, p, prho, prr, pt, q, qc, ql,&
-               qr, qs, qsws, qswst, rho, s, sa, ss, ssws, sswst, saswsb,       &
+               qr, qs, qsws, qswst, rho_ocean, s, sa, ss, ssws, sswst, saswsb,       &
                saswst, shf, td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q,        &
                time_vert, ts, tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, &
@@ -926,5 +926,5 @@
                                         rmask(j,i,sr)
                    ENDIF
-                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) *            &
+                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho_ocean(k,j,i) *            &
                                                        rmask(j,i,sr)
                    sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) *           &
@@ -1400,5 +1400,5 @@
        hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
        hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
-       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
+       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho_ocean )/dz
        hom(:,1,58,sr) = sums(:,58)     ! u"pt"
        hom(:,1,59,sr) = sums(:,59)     ! u*pt*
@@ -1407,10 +1407,10 @@
        hom(:,1,62,sr) = sums(:,62)     ! v*pt*
        hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
-       hom(:,1,64,sr) = sums(:,64)     ! rho
+       hom(:,1,64,sr) = sums(:,64)     ! rho_ocean
        hom(:,1,65,sr) = sums(:,65)     ! w"sa"
        hom(:,1,66,sr) = sums(:,66)     ! w*sa*
        hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
        hom(:,1,68,sr) = sums(:,68)     ! w*p*
-       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
+       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho_ocean
        hom(:,1,70,sr) = sums(:,70)     ! q*2
        hom(:,1,71,sr) = sums(:,71)     ! prho
@@ -1706,5 +1706,5 @@
     USE arrays_3d,                                                             &
         ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
-               qsws, qswst, rho, s, sa, saswsb, saswst, shf, ss, ssws, sswst,  &
+               qsws, qswst, rho_ocean, s, sa, saswsb, saswst, shf, ss, ssws, sswst,  &
                td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts,      &
                tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, vswst, w,      &
@@ -2786,5 +2786,5 @@
           ENDIF
 
-          !$acc parallel loop gang present( rflags_invers, rho, prho, rmask, sums_l ) create( s1, s2 )
+          !$acc parallel loop gang present( rflags_invers, rho_ocean, prho, rmask, sums_l ) create( s1, s2 )
           DO  k = nzb, nzt_diff
              s1 = 0
@@ -2793,5 +2793,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   s1 = s1 + rho(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                   s1 = s1 + rho_ocean(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
                    s2 = s2 + prho(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
                 ENDDO
@@ -3424,5 +3424,5 @@
        hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
        hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
-       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
+       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho_ocean )/dz
        hom(:,1,58,sr) = sums(:,58)     ! u"pt"
        hom(:,1,59,sr) = sums(:,59)     ! u*pt*
@@ -3431,10 +3431,10 @@
        hom(:,1,62,sr) = sums(:,62)     ! v*pt*
        hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
-       hom(:,1,64,sr) = sums(:,64)     ! rho
+       hom(:,1,64,sr) = sums(:,64)     ! rho_ocean
        hom(:,1,65,sr) = sums(:,65)     ! w"sa"
        hom(:,1,66,sr) = sums(:,66)     ! w*sa*
        hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
        hom(:,1,68,sr) = sums(:,68)     ! w*p*
-       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
+       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho_ocean
        hom(:,1,70,sr) = sums(:,70)     ! q*2
        hom(:,1,71,sr) = sums(:,71)     ! prho
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 2030)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! renamed variable rho to rho_ocean
 ! 
 ! Former revisions:
@@ -623,5 +623,5 @@
 #if defined( __nopointer )
        ALLOCATE( prho(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
-                 rho(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),                          &
@@ -634,5 +634,5 @@
                  sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
        prho => prho_1
-       rho  => rho_1  ! routines calc_mean_profile and diffusion_e require
+       rho_ocean  => rho_1  ! routines calc_mean_profile and diffusion_e require
                       ! density to be apointer
 #endif
Index: palm/trunk/SOURCE/init_masks.f90
===================================================================
--- palm/trunk/SOURCE/init_masks.f90	(revision 2030)
+++ palm/trunk/SOURCE/init_masks.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! renamed variable rho to rho_ocean
 ! 
 ! Former revisions:
@@ -350,5 +350,5 @@
                 unit = 'kg/kg'
 
-             CASE ( 'rho' )
+             CASE ( 'rho_ocean' )
                 IF ( .NOT. ocean )  THEN
                    WRITE ( message_string, * ) 'output of "', TRIM( var ), &
Index: palm/trunk/SOURCE/init_ocean.f90
===================================================================
--- palm/trunk/SOURCE/init_ocean.f90	(revision 2030)
+++ palm/trunk/SOURCE/init_ocean.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! renamed variable rho_init to rho_ocean_init
 ! 
 ! Former revisions:
@@ -90,5 +90,5 @@
     REAL(wp)     ::  sa_l !<
 
-    REAL(wp), DIMENSION(nzb:nzt+1) ::  rho_init !<
+    REAL(wp), DIMENSION(nzb:nzt+1) ::  rho_ocean_init !<
 
     ALLOCATE( hyp(nzb:nzt+1) )
@@ -105,5 +105,5 @@
 
     hyp(nzt)      = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1)
-    rho_init(nzt) = rho_surface
+    rho_ocean_init(nzt) = rho_surface
 
     DO  k = nzt-1, 1, -1
@@ -124,7 +124,7 @@
           pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) )
 
-          rho_init(k) = eqn_state_seawater_func( hyp(k), pt_l, sa_l )
+          rho_ocean_init(k) = eqn_state_seawater_func( hyp(k), pt_l, sa_l )
 
-          rho_reference = rho_reference + rho_init(k) * dzu(k+1)
+          rho_reference = rho_reference + rho_ocean_init(k) * dzu(k+1)
 
        ENDDO
@@ -133,5 +133,5 @@
 
        DO  k = nzt-1, 0, -1
-          hyp(k) = hyp(k+1) + g * 0.5_wp * ( rho_init(k) + rho_init(k+1 ) ) * &
+          hyp(k) = hyp(k+1) + g * 0.5_wp * ( rho_ocean_init(k) + rho_ocean_init(k+1 ) ) * &
                               dzu(k+1)
        ENDDO
@@ -161,5 +161,5 @@
 !
 !-- Store initial density profile
-    hom(:,1,77,:)  = SPREAD( rho_init(:), 2, statistic_regions+1 )
+    hom(:,1,77,:)  = SPREAD( rho_ocean_init(:), 2, statistic_regions+1 )
 
 !
@@ -168,5 +168,5 @@
        ref_state(:) = prho_reference
     ELSE
-       ref_state(:) = rho_init(:)
+       ref_state(:) = rho_ocean_init(:)
     ENDIF
 
Index: palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2030)
+++ palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! renamed variable rho to rho_ocean
 ! 
 ! Former revisions:
@@ -257,5 +257,5 @@
                 residual_moisture = 9999999.9_wp,       & !< NAMELIST m_res
                 rho_cp,                                 & !< rho_surface * cp
-                rho_lv,                                 & !< rho * l_v
+                rho_lv,                                 & !< rho_ocean * l_v
                 rd_d_rv,                                & !< r_d / r_v
                 saturation_moisture = 9999999.9_wp,     & !< NAMELIST m_sat
Index: palm/trunk/SOURCE/microphysics_mod.f90
===================================================================
--- palm/trunk/SOURCE/microphysics_mod.f90	(revision 2030)
+++ palm/trunk/SOURCE/microphysics_mod.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! renamed variable rho to rho_ocean
 ! 
 ! Former revisions:
@@ -82,5 +82,5 @@
 !
 ! 1241 2013-10-30 11:36:58Z heinze
-! hyp and rho have to be calculated at each time step if data from external
+! hyp and rho_ocean have to be calculated at each time step if data from external
 ! file LSF_DATA are used
 !
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 2030)
+++ palm/trunk/SOURCE/modules.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! renamed variable rho to rho_ocean and rho_av to rho_ocean_av
 ! 
 ! Former revisions:
@@ -473,5 +473,5 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
           e, e_p, nr, nr_p, p, prho, pt, pt_p, q, q_p, qc, ql, ql_c, ql_v,     &
-          ql_vp, qr, qr_p, rho, s, s_p, sa, sa_p, te_m, tnr_m, tpt_m, tq_m,    &
+          ql_vp, qr, qr_p, rho_ocean, s, s_p, sa, sa_p, te_m, tnr_m, tpt_m, tq_m,    &
           tqr_m, ts_m, tsa_m, tu_m, tv_m, tw_m, u, u_p, v, v_p, vpt, w, w_p
 #else
@@ -484,5 +484,5 @@
     REAL(wp), DIMENSION(:,:,:), POINTER ::                                     &
           e, e_p, nr, nr_p, prho, pt, pt_p, q, q_p, qc, ql, ql_c, qr, qr_p,    &
-          rho, s, s_p, sa, sa_p, te_m, tnr_m, tpt_m, tq_m, tqr_m, ts_m,        &
+          rho_ocean, s, s_p, sa, sa_p, te_m, tnr_m, tpt_m, tq_m, tqr_m, ts_m,        &
           tsa_m, tu_m, tv_m, tw_m, u, u_p, v, v_p, vpt, w, w_p
 #endif
@@ -519,5 +519,5 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
           e_av, lpt_av, nr_av, p_av, pc_av, pr_av, prr_av, pt_av, q_av, qc_av, &
-          ql_av, ql_c_av, ql_v_av, ql_vp_av, qr_av, qv_av, rho_av, s_av, sa_av,&
+          ql_av, ql_c_av, ql_v_av, ql_vp_av, qr_av, qv_av, rho_ocean_av, s_av, sa_av,&
           u_av, v_av, vpt_av, w_av
   
Index: palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 2030)
+++ palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! renamed variable rho to rho_ocean
 ! 
 ! Former revisions:
@@ -778,5 +778,5 @@
                 CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',    &
                        'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',        &
-                       'rho', 's', 'sa', 'vpt'  )
+                       'rho_ocean', 's', 'sa', 'vpt'  )
 
                    grid_x = 'x'
@@ -1299,5 +1299,5 @@
 !--             Most variables are defined on the scalar grid
                 CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',    &
-                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho', &
+                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho_ocean', &
                        's', 'sa', 'vpt' )
 
@@ -4403,5 +4403,5 @@
 !--             Most variables are defined on the zu levels
                 CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',   &
-                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho',&
+                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho_ocean',&
                        's', 'sa', 'u', 'v', 'vpt' )
 
Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 2030)
+++ palm/trunk/SOURCE/production_e.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! renamed variable rho to rho_ocean
 ! 
 ! Former revisions:
@@ -136,5 +136,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzw, dd2zu, kh, km, pt, q, ql, qsws, qswst, rho, shf,       &
+           ONLY:  ddzw, dd2zu, kh, km, pt, q, ql, qsws, qswst, rho_ocean, shf,       &
                   tend, tswst, u, v, vpt, w
 
@@ -537,5 +537,5 @@
                             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) ) * &
+                                          ( rho_ocean(k+1,j,i) - rho_ocean(k-1,j,i) ) * &
                                           dd2zu(k)
                          ENDDO
@@ -576,6 +576,6 @@
                          DO  k = nzb_s_inner(j,i)+1, nzt
                             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) ) * &
+                                          kh(k,j,i) * g / rho_ocean(k,j,i) *      &
+                                          ( rho_ocean(k+1,j,i) - rho_ocean(k-1,j,i) ) * &
                                           dd2zu(k)
                          ENDDO
@@ -742,5 +742,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzw, dd2zu, kh, km, pt, q, ql, qsws, qswst, rho, shf,       &
+           ONLY:  ddzw, dd2zu, kh, km, pt, q, ql, qsws, qswst, rho_ocean, shf,       &
                   tend, tswst, u, v, vpt, w
 
@@ -806,5 +806,5 @@
 !--    Calculate TKE production by shear
        !$acc kernels present( ddzw, dd2zu, kh, km, nzb_diff_s_inner, nzb_diff_s_outer ) &
-       !$acc         present( nzb_s_inner, pt, q, ql, qsws, qswst, rho )                &
+       !$acc         present( nzb_s_inner, pt, q, ql, qsws, qswst, rho_ocean )                &
        !$acc         present( shf, tend, tswst, u, v, vpt, w, wall_e_x, wall_e_y )      &
        !$acc         copyin( u_0, v_0 )
@@ -1166,5 +1166,5 @@
                                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) ) * &
+                                             ( rho_ocean(k+1,j,i) - rho_ocean(k-1,j,i) ) * &
                                              dd2zu(k)
                             ENDIF
@@ -1212,6 +1212,6 @@
                             IF ( k > nzb_s_inner(j,i) )  THEN
                                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) ) * &
+                                             kh(k,j,i) * g / rho_ocean(k,j,i) *      &
+                                             ( rho_ocean(k+1,j,i) - rho_ocean(k-1,j,i) ) * &
                                              dd2zu(k)
                             ENDIF
@@ -1410,5 +1410,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzw, dd2zu, kh, km, pt, q, ql, qsws, qswst, rho, shf,       &
+           ONLY:  ddzw, dd2zu, kh, km, pt, q, ql, qsws, qswst, rho_ocean, shf,       &
                   tend, tswst, u, v, vpt, w
 
@@ -1766,5 +1766,5 @@
                       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)
+                                    ( rho_ocean(k+1,j,i) - rho_ocean(k-1,j,i) ) * dd2zu(k)
                    ENDDO
 
@@ -1797,6 +1797,6 @@
                    DO  k = nzb_s_inner(j,i)+1, nzt
                       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)
+                                    kh(k,j,i) * g / rho_ocean(k,j,i) * &
+                                    ( rho_ocean(k+1,j,i) - rho_ocean(k-1,j,i) ) * dd2zu(k)
                    ENDDO
 
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 2030)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! renamed variable rho to rho_ocean
 ! 
 ! Former revisions:
@@ -215,5 +215,5 @@
                nr, nr_p, nrsws, nrswst, pt, ptdf_x, ptdf_y, pt_init, pt_p,     &
                prho, q, q_init, q_p, qsws, qswst, qr, qr_p, qrsws, qrswst, rdf,&
-               rdf_sc, ref_state, rho, s, s_init, s_p, sa, sa_init, sa_p,      &
+               rdf_sc, ref_state, rho_ocean, s, s_init, s_p, sa, sa_init, sa_p,      &
                saswsb, saswst, shf, ssws, sswst, tend,                         &
                te_m, tnr_m, tpt_m, tq_m, tqr_m, ts_m, tsa_m, tswst, tu_m, tv_m,&
@@ -562,5 +562,5 @@
           IF ( .NOT. neutral )  THEN
              IF ( ocean )  THEN
-                CALL buoyancy( i, j, rho, 3 )
+                CALL buoyancy( i, j, rho_ocean, 3 )
              ELSE
                 IF ( .NOT. humidity )  THEN
@@ -1267,5 +1267,5 @@
     IF ( .NOT. neutral )  THEN
        IF ( ocean )  THEN
-          CALL buoyancy( rho, 3 )
+          CALL buoyancy( rho_ocean, 3 )
        ELSE
           IF ( .NOT. humidity )  THEN
@@ -2205,5 +2205,5 @@
     IF ( .NOT. neutral )  THEN
        IF ( ocean )  THEN
-          CALL buoyancy( rho, 3 )
+          CALL buoyancy( rho_ocean, 3 )
        ELSE
           IF ( .NOT. humidity )  THEN
Index: palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/read_3d_binary.f90	(revision 2030)
+++ palm/trunk/SOURCE/read_3d_binary.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! renamed variable rho_av to rho_ocean_av
 ! 
 ! Former revisions:
@@ -712,10 +712,10 @@
                    IF ( k == 1 )  READ ( 13 )  seq_random_array
 
-                CASE ( 'rho_av' )
-                   IF ( .NOT. ALLOCATED( rho_av ) )  THEN
-                      ALLOCATE( rho_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   rho_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                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)
 
Index: palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 2030)
+++ palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! renamed variable rho to rho_ocean and rho_av to rho_ocean_av
 ! 
 ! Former revisions:
@@ -130,5 +130,5 @@
     USE arrays_3d,                                                             &
         ONLY:  dzw, e, nr, ol, p, pt, precipitation_rate, q, qc, ql, ql_c,     &
-               ql_v, qr, qsws, rho, s, sa, shf, ssws, ts, u, us, v, vpt, w, z0,&
+               ql_v, qr, qsws, rho_ocean, s, sa, shf, ssws, ts, u, us, v, vpt, w, z0,&
                z0h, z0q
 
@@ -136,5 +136,5 @@
         ONLY:  e_av, lpt_av, lwp_av, nr_av, ol_av, p_av, pc_av, pr_av, prr_av, &
                precipitation_rate_av, pt_av, q_av, qc_av, ql_av, ql_c_av,      &
-               ql_v_av, ql_vp_av, qr_av, qsws_av, qv_av, rho_av, s_av, sa_av,  &
+               ql_v_av, ql_vp_av, qr_av, qsws_av, qv_av, rho_ocean_av, s_av, sa_av,  &
                shf_av, ssws_av, ts_av, u_av, us_av, v_av, vpt_av, w_av, z0_av, &
                z0h_av, z0q_av
@@ -323,9 +323,9 @@
                 qv_av = 0.0_wp
 
-             CASE ( 'rho' )
-                IF ( .NOT. ALLOCATED( rho_av ) )  THEN
-                   ALLOCATE( rho_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                rho_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' )
@@ -657,9 +657,9 @@
              ENDDO
 
-          CASE ( 'rho' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
-                      rho_av(k,j,i) = rho_av(k,j,i) + rho(k,j,i)
+          CASE ( 'rho_ocean' )
+             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
Index: palm/trunk/SOURCE/surface_coupler.f90
===================================================================
--- palm/trunk/SOURCE/surface_coupler.f90	(revision 2030)
+++ palm/trunk/SOURCE/surface_coupler.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! renamed variable rho to rho_ocean
 ! 
 ! Former revisions:
@@ -80,5 +80,5 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  pt, shf, qsws, qswst_remote, rho, sa, saswst, total_2d_a,       &
+        ONLY:  pt, shf, qsws, qswst_remote, rho_ocean, sa, saswst, total_2d_a,       &
                total_2d_o, tswst, u, usws, uswst, v, vsws, vswst
 
@@ -421,5 +421,5 @@
 !--       S'w' = -S * evaporation / ( rho_water * ( 1 - S ) )
           saswst = -1.0_wp * sa(nzt,:,:) * 0.001_wp * qswst_remote /  &
-                    ( rho(nzt,:,:) * ( 1.0_wp - sa(nzt,:,:) * 0.001_wp ) )
+                    ( rho_ocean(nzt,:,:) * ( 1.0_wp - sa(nzt,:,:) * 0.001_wp ) )
        ENDIF
 
@@ -428,10 +428,10 @@
 !--    (constants are the specific heat capacities for air and water)
 !--    now tswst is the ocean top heat flux
-       tswst = tswst / rho(nzt,:,:) * cp / cpw
+       tswst = tswst / rho_ocean(nzt,:,:) * cp / cpw
 
 !
 !--    Adjust the momentum fluxes with respect to ocean density
-       uswst = uswst / rho(nzt,:,:)
-       vswst = vswst / rho(nzt,:,:)
+       uswst = uswst / rho_ocean(nzt,:,:)
+       vswst = vswst / rho_ocean(nzt,:,:)
 
     ENDIF
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 2030)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! renamed variable rho to rho_ocean
 ! 
 ! Former revisions:
@@ -242,5 +242,5 @@
     USE arrays_3d,                                                             &
         ONLY:  diss, dzu, e, e_p, nr_p, prho, pt, pt_p, pt_init, q_init, q,    &
-               ql, ql_c, ql_v, ql_vp, qr_p, q_p, ref_state, rho, s, s_p, sa_p, &
+               ql, ql_c, ql_v, ql_vp, qr_p, q_p, ref_state, rho_ocean, s, s_p, sa_p, &
                tend, u, u_p, v, vpt, v_p, w, w_p
 
@@ -482,5 +482,5 @@
              ENDIF
              IF ( ocean )  THEN
-                CALL calc_mean_profile( rho, 64 )
+                CALL calc_mean_profile( rho_ocean, 64 )
                 ref_state(:)  = hom(:,1,64,0)
              ENDIF
@@ -552,5 +552,5 @@
                 IF ( ocean )  THEN
                    CALL exchange_horiz( sa_p, nbgp )
-                   CALL exchange_horiz( rho, nbgp )
+                   CALL exchange_horiz( rho_ocean, nbgp )
                    CALL exchange_horiz( prho, nbgp )
                 ENDIF
@@ -620,5 +620,5 @@
                 IF ( ocean )  THEN
                    CALL exchange_horiz( sa_p, nbgp )
-                   CALL exchange_horiz( rho, nbgp )
+                   CALL exchange_horiz( rho_ocean, nbgp )
                   CALL exchange_horiz( prho, nbgp )
                 ENDIF
@@ -710,5 +710,5 @@
              IF ( ocean )  THEN
                 CALL exchange_horiz( sa_p, nbgp )
-                CALL exchange_horiz( rho, nbgp )
+                CALL exchange_horiz( rho_ocean, nbgp )
                 CALL exchange_horiz( prho, nbgp )
              ENDIF
Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 2030)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 2031)
@@ -21,5 +21,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! renamed variable rho to rho_ocean
 ! 
 ! Former revisions:
@@ -419,5 +419,5 @@
     INTEGER(iwp), PARAMETER                        :: icsurf   = 5         !< Surface skin layer heat capacity (J mâ2 Kâ1 )
     INTEGER(iwp), PARAMETER                        :: ithick   = 6         !< thickness of the surface (wall, roof, land)  ( m )
-    INTEGER(iwp), PARAMETER                        :: irhoC    = 7         !< volumetric heat capacity rho*C of the material ( J mâ3 Kâ1 )
+    INTEGER(iwp), PARAMETER                        :: irhoC    = 7         !< volumetric heat capacity rho_ocean*C of the material ( J mâ3 Kâ1 )
     INTEGER(iwp), PARAMETER                        :: ilambdah = 8         !< thermal conductivity Î»H of the wall (W mâ1 Kâ1 )
     CHARACTER(12), DIMENSION(:), ALLOCATABLE       :: surface_type_names   !< names of wall types (used only for reports)
@@ -2886,5 +2886,5 @@
         IF ( plant_canopy )  THEN
             pchf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp &
-                        / (cp * hyp(nzub:nzut) * dx*dy*dz) !< equals to 1 / (rho * c_p * Vbox * T)
+                        / (cp * hyp(nzub:nzut) * dx*dy*dz) !< equals to 1 / (rho_ocean * c_p * Vbox * T)
         ENDIF
 
@@ -3863,5 +3863,5 @@
             ENDIF
             
-!--         volumetric heat capacity rho*C of the wall ( J mâ3 Kâ1 )
+!--         volumetric heat capacity rho_ocean*C of the wall ( J mâ3 Kâ1 )
             rho_c_wall(:,l) = surface_params(irhoC, ip)
             
@@ -3930,5 +3930,5 @@
             pt1  = pt(k,j,i)
 
-!--         calculate rho * cp coefficient at surface layer
+!--         calculate rho_ocean * cp coefficient at surface layer
             rho_cp  = cp * hyp(k) / ( r_d * pt1 * exn(k) )
 
Index: palm/trunk/SOURCE/write_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/write_3d_binary.f90	(revision 2030)
+++ palm/trunk/SOURCE/write_3d_binary.f90	(revision 2031)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! renamed variable rho_av to rho_ocean_av
 ! 
 ! Former revisions:
@@ -256,6 +256,6 @@
     ENDIF    
     IF ( ocean )  THEN
-       IF ( ALLOCATED( rho_av ) )  THEN
-          WRITE ( 14 )  'rho_av              ';  WRITE ( 14 )  rho_av
+       IF ( ALLOCATED( rho_ocean_av ) )  THEN
+          WRITE ( 14 )  'rho_ocean_av              ';  WRITE ( 14 )  rho_ocean_av
        ENDIF
        WRITE ( 14 )  'sa                  ';  WRITE ( 14 )  sa
