Index: palm/trunk/SOURCE/poismg_mod.f90
===================================================================
--- palm/trunk/SOURCE/poismg_mod.f90	(revision 2022)
+++ palm/trunk/SOURCE/poismg_mod.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Anelastic approximation implemented (stll error in optimized multigrid)
 ! 
 ! Former revisions:
@@ -262,5 +262,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  f1_mg, f2_mg, f3_mg
+           ONLY:  f1_mg, f2_mg, f3_mg, rho_air_mg
 
        USE control_parameters,                                                 &
@@ -309,7 +309,7 @@
                 kp1 = k-ind_even_odd
                 r(k,j,i) = f_mg(k,j,i)                                         &
-                      - ddx2_mg(l) *                                           &
+                      - rho_air_mg(k,l) * ddx2_mg(l) *                         &
                       ( p_mg(k,j,i+1) +  p_mg(k,j,i-1)  )                      &
-                      - ddy2_mg(l) *                                           &
+                      - rho_air_mg(k,l) * ddy2_mg(l) *                         &
                       ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )                       &
                       - f2_mg_b(k,l) * p_mg(kp1,j,i)                           &
@@ -322,7 +322,7 @@
                 kp1 = k+ind_even_odd+1
                 r(k,j,i) = f_mg(k,j,i)                                         &
-                      - ddx2_mg(l) *                                           &
+                      - rho_air_mg(k,l) * ddx2_mg(l) *                         &
                       ( p_mg(k,j,i+1) +  p_mg(k,j,i-1)  )                      &
-                      - ddy2_mg(l) *                                           &
+                      - rho_air_mg(k,l) * ddy2_mg(l) *                         &
                       ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )                       &
                       - f2_mg_b(k,l) * p_mg(kp1,j,i)                           &
@@ -693,5 +693,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  f1_mg, f2_mg, f3_mg
+           ONLY:  f1_mg, f2_mg, f3_mg, rho_air_mg
 
        USE control_parameters,                                                 &
@@ -755,7 +755,7 @@
                          kp1 = k-ind_even_odd
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -774,7 +774,7 @@
                          kp1 = k-ind_even_odd
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -793,7 +793,7 @@
                          kp1 = k+ind_even_odd+1
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -812,7 +812,7 @@
                          kp1 = k+ind_even_odd+1
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -843,7 +843,7 @@
                          j   = jj
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -852,7 +852,7 @@
                          j = jj+2
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -869,7 +869,7 @@
                          j   = jj
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -878,7 +878,7 @@
                          j = jj+2
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -895,7 +895,7 @@
                          j   = jj
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -904,7 +904,7 @@
                          j = jj+2
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -921,7 +921,7 @@
                          j   = jj
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
@@ -930,7 +930,7 @@
                          j = jj+2
                          p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
-                                 ddx2_mg(l) *                                  &
+                                 rho_air_mg(k,l) * ddx2_mg(l) *                &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
-                               + ddy2_mg(l) *                                  &
+                               + rho_air_mg(k,l) * ddy2_mg(l) *                &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
