Changeset 2252 for palm/trunk/SOURCE/init_3d_model.f90

Timestamp:

Jun 7, 2017 9:35:37 AM (8 years ago)

Author:

knoop

Message:

air density now depending on surface_pressure even in boussinesq mode

File:

• palm/trunk/SOURCE/init_3d_model.f90 (modified) (3 diffs)

palm/trunk/SOURCE/init_3d_model.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! rho_air now depending on surface_pressure even in Boussinesq mode
+! 
+! 2233 2017-05-30 18:08:54Z suehring
 !
 ! 2232 2017-05-30 17:47:52Z suehring
@@ -657 +660 @@
 !
 !-- Density profile calculation for anelastic approximation
+    t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**( r_d / cp )
     IF ( TRIM( approximation ) == 'anelastic' ) THEN
-       t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**( r_d / cp )
        DO  k = nzb, nzt+1
           p_hydrostatic(k)    = surface_pressure * 100.0_wp *                  &
@@ -674 +677 @@
                             + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
     ELSE
-       rho_air     = 1.0_wp
-       rho_air_zw  = 1.0_wp
+       DO  k = nzb, nzt+1
+          p_hydrostatic(k)    = surface_pressure * 100.0_wp *                  &
+                                ( 1 - ( g * zu(nzb) ) / ( cp * t_surface )       &
+                                )**( cp / r_d )
+          rho_air(k)          = ( p_hydrostatic(k) *                           &
+                                  ( 100000.0_wp / p_hydrostatic(k)             &
+                                  )**( r_d / cp )                              &
+                                ) / ( r_d * pt_init(nzb) )
+       ENDDO
+       DO  k = nzb, nzt
+          rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
+       ENDDO
+       rho_air_zw(nzt+1)  = rho_air_zw(nzt)                                    &
+                            + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
     ENDIF