Changeset 2037 for palm/trunk/SOURCE/sor.f90

Timestamp:

Oct 26, 2016 11:15:40 AM (7 years ago)

Author:

knoop

Message:

Anelastic approximation implemented

File:

• palm/trunk/SOURCE/sor.f90 (modified) (7 diffs)

palm/trunk/SOURCE/sor.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -59 +59 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE sor( d, ddzu, ddzw, p )
- 
+
+    USE arrays_3d,                                                             &
+        ONLY:  rho_air, rho_air_zw
 
     USE grid_variables,                                                        &
@@ -101 +103 @@
 !-- Compute pre-factors.
     DO  k = 1, nz
-         f2(k) = ddzu(k+1) * ddzw(k)
-         f3(k) = ddzu(k)   * ddzw(k)
-         f1(k) = 2.0_wp * ( ddx2 + ddy2 ) + f2(k) + f3(k)
+         f2(k) = ddzu(k+1) * ddzw(k) * rho_air_zw(k)
+         f3(k) = ddzu(k)   * ddzw(k) * rho_air_zw(k-1)
+         f1(k) = 2.0_wp * ( ddx2 + ddy2 ) * rho_air(k) + f2(k) + f3(k)
     ENDDO
 
@@ -131 +133 @@
              DO  k = nzb+1, nzt
                 p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
-                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
-                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
-                               f2(k) * p(k+1,j,i)                 +    &
-                               f3(k) * p(k-1,j,i)                 -    &
-                               d(k,j,i)                           -    &
-                               f1(k) * p(k,j,i)           )
+                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +   &
+                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +   &
+                           f2(k) * p(k+1,j,i)                              +   &
+                           f3(k) * p(k-1,j,i)                              -   &
+                           d(k,j,i)                                        -   &
+                           f1(k) * p(k,j,i)           )
              ENDDO
           ENDDO
@@ -144 +146 @@
           DO  j = nys1, nyn, 2
              DO  k = nzb+1, nzt
-                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
-                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
-                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
-                               f2(k) * p(k+1,j,i)                 +    &
-                               f3(k) * p(k-1,j,i)                 -    &
-                               d(k,j,i)                           -    &
-                               f1(k) * p(k,j,i)           )
+                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (                    &
+                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +   &
+                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +   &
+                           f2(k) * p(k+1,j,i)                              +   &
+                           f3(k) * p(k-1,j,i)                              -   &
+                           d(k,j,i)                                        -   &
+                           f1(k) * p(k,j,i)           )
              ENDDO
           ENDDO
@@ -176 +178 @@
              DO  k = nzb+1, nzt
                 p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
-                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
-                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
-                               f2(k) * p(k+1,j,i)                 +    &
-                               f3(k) * p(k-1,j,i)                 -    &
-                               d(k,j,i)                           -    &
-                               f1(k) * p(k,j,i)           )
+                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +   &
+                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +   &
+                           f2(k) * p(k+1,j,i)                              +   &
+                           f3(k) * p(k-1,j,i)                              -   &
+                           d(k,j,i)                                        -   &
+                           f1(k) * p(k,j,i)           )
              ENDDO
           ENDDO
@@ -190 +192 @@
              DO  k = nzb+1, nzt
                 p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
-                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
-                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
-                               f2(k) * p(k+1,j,i)                 +    &
-                               f3(k) * p(k-1,j,i)                 -    &
-                               d(k,j,i)                           -    &
-                               f1(k) * p(k,j,i)           )
+                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +   &
+                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +   &
+                           f2(k) * p(k+1,j,i)                              +   &
+                           f3(k) * p(k-1,j,i)                              -   &
+                           d(k,j,i)                                        -   &
+                           f1(k) * p(k,j,i)           )
              ENDDO
           ENDDO