Index: palm/trunk/SOURCE/diffusion_v.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_v.f90	(revision 2001)
+++ palm/trunk/SOURCE/diffusion_v.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -111,5 +111,6 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w
+           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w,                  &
+                  drho_air, rho_air_zw
        
        USE control_parameters,                                                 &
@@ -212,9 +213,9 @@
                       & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
                       &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
-                      &            )                                           &
+                      &            ) * rho_air_zw(k)                           &
                       &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
                       &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
-                      &            )                                           &
-                      &   ) * ddzw(k)
+                      &            ) * rho_air_zw(k-1)                         &
+                      &   ) * ddzw(k) * drho_air(k)
              ENDDO
 
@@ -237,9 +238,9 @@
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
-                      &   ) * ddzw(k)                                          &
-                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
-                      &   + vsws(j,i)                                          &
-                      &   ) * ddzw(k)
+                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
+                      &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
+                      &            ) * rho_air_zw(k)                           &
+                      &   - ( -vsws(j,i) )                                     &
+                      &   ) * ddzw(k) * drho_air(k)
              ENDIF
 
@@ -255,9 +256,9 @@
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
-                      &   ) * ddzw(k)                                          &
-                      & + ( -vswst(j,i)                                        &
-                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
-                      &   ) * ddzw(k)
+                      & + ( ( -vswst(j,i) )                                    &
+                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
+                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
+                      &            ) * rho_air_zw(k-1)                         &
+                      &   ) * ddzw(k) * drho_air(k)
              ENDIF
 
@@ -276,5 +277,6 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w
+           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w,                  &
+                  drho_air, rho_air_zw
        
        USE control_parameters,                                                 &
@@ -385,9 +387,9 @@
                          & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)&
                          &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy      &
-                         &            )                                        &
+                         &            ) * rho_air_zw(k)                        &
                          &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)&
                          &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy    &
-                         &            )                                        &
-                         &   ) * ddzw(k)
+                         &            ) * rho_air_zw(k-1)                      &
+                         &   ) * ddzw(k) * drho_air(k)
                 ENDIF
              ENDDO
@@ -418,9 +420,9 @@
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
-                      &   ) * ddzw(k)                                          &
-                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
-                      &   + vsws(j,i)                                          &
-                      &   ) * ddzw(k)
+                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
+                      &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
+                      &            ) * rho_air_zw(k)                           &
+                      &   - ( -vsws(j,i) )                                     &
+                      &   ) * ddzw(k) * drho_air(k)
              ENDDO
           ENDDO
@@ -444,9 +446,9 @@
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
-                      &   ) * ddzw(k)                                          &
-                      & + ( -vswst(j,i)                                        &
-                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
-                      &   ) * ddzw(k)
+                      & + ( ( -vswst(j,i) )                                    &
+                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
+                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
+                      &            ) * rho_air_zw(k-1)                         &
+                      &   ) * ddzw(k) * drho_air(k)
              ENDDO
           ENDDO
@@ -466,5 +468,6 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w
+           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w,                  &
+                  drho_air, rho_air_zw
        
        USE control_parameters,                                                 &
@@ -556,9 +559,9 @@
                       & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
                       &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
-                      &            )                                           &
+                      &            ) * rho_air_zw(k)                           &
                       &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
                       &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
-                      &            )                                           &
-                      &   ) * ddzw(k)
+                      &            ) * rho_air_zw(k-1)                         &
+                      &   ) * ddzw(k) * drho_air(k)
        ENDDO
 
@@ -580,9 +583,9 @@
 
           tend(k,j,i) = tend(k,j,i)                                            &
-                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
-                      &   ) * ddzw(k)                                          &
-                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
-                      &   + vsws(j,i)                                          &
-                      &   ) * ddzw(k)
+                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
+                      &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
+                      &            ) * rho_air_zw(k)                           &
+                      &   - ( -vsws(j,i) )                                     &
+                      &   ) * ddzw(k) * drho_air(k)
        ENDIF
 
@@ -594,13 +597,12 @@
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
-          kmzm = 0.25_wp * &
-                 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
+          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
 
           tend(k,j,i) = tend(k,j,i)                                            &
-                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
-                      &   ) * ddzw(k)                                          &
-                      & + ( -vswst(j,i)                                        &
-                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
-                      &   ) * ddzw(k)
+                      & + ( ( -vswst(j,i) )                                    &
+                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
+                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
+                      &            ) * rho_air_zw(k-1)                         &
+                      &   ) * ddzw(k) * drho_air(k)
        ENDIF
 
