Changeset 4088

Timestamp:

Jul 11, 2019 1:57:56 PM (5 years ago)

Author:

Giersch

Message:

Pressure and density profile calculations revised using basic functions, comments improved, function for ideal gas law revised

Location:

palm/trunk/SOURCE

Files:

• basic_constants_and_equations_mod.f90 (modified) (5 diffs)
• init_3d_model.f90 (modified) (2 diffs)

palm/trunk/SOURCE/basic_constants_and_equations_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Comment of barometric formula improved, function for ideal gas law revised 
+! 
+! 4084 2019-07-10 17:09:11Z knoop
 ! Changed precomputed fractions to be variable based
 ! 
@@ -236 +239 @@
 !
 !--    compute density according to ideal gas law:
-       ideal_gas_law_rho_pt_0d = p / (r_d * (1.0_wp / exner_function_invers(p)) * t)
+       ideal_gas_law_rho_pt_0d = p / (r_d * exner_function(p) * t)
 
     END FUNCTION ideal_gas_law_rho_pt_0d
@@ -255 +258 @@
 !
 !--    compute density according to ideal gas law:
-       ideal_gas_law_rho_pt_1d = p / (r_d * (1.0_wp / exner_function_invers(p)) * t)
+       ideal_gas_law_rho_pt_1d = p / (r_d * exner_function(p) * t)
 
     END FUNCTION ideal_gas_law_rho_pt_1d
@@ -334 +337 @@
 ! Description:
 ! ------------
-!> Compute the barometric formula for scalar arguments.
+!> Compute the barometric formula for scalar arguments. The calculation is 
+!> based on the assumption of a polytropic atmosphere and neutral 
+!> stratification, where the temperature lapse rate is g/cp. 
 !------------------------------------------------------------------------------!
     FUNCTION barometric_formula_0d( z, t_0, p_0)
@@ -354 +359 @@
 ! Description:
 ! ------------
-!> Compute the barometric formula for 1-D array arguments.
+!> Compute the barometric formula for 1-D array arguments. The calculation is 
+!> based on the assumption of a polytropic atmosphere and neutral 
+!> stratification, where the temperature lapse rate is g/cp.
 !------------------------------------------------------------------------------!
     FUNCTION barometric_formula_1d( z, t_0, p_0)

palm/trunk/SOURCE/init_3d_model.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Pressure and density profile calculations revised using basic functions
+! 
+! 4048 2019-06-21 21:00:21Z knoop
 ! Further modularization of particle code components
 ! 
@@ -845 +848 @@
     ALLOCATE( drho_air_zw(nzb:nzt+1) )
 !
-!-- Density profile calculation for anelastic approximation
-    t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**( r_d / c_p )
-    IF ( TRIM( approximation ) == 'anelastic' ) THEN
+!-- Density profile calculation for anelastic and Boussinesq approximation
+!-- In case of a Boussinesq approximation, a constant density is calculated
+!-- mainly for output purposes. This density do not need to be considered 
+!-- in the model's system of equations.
+!    t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**( r_d / c_p )
+    IF ( TRIM( approximation ) == 'anelastic' )  THEN
        DO  k = nzb, nzt+1
-          p_hydrostatic(k)    = surface_pressure * 100.0_wp *                  &
-                                ( 1 - ( g * zu(k) ) / ( c_p * t_surface )      &
-                                )**( c_p / r_d )
-          rho_air(k)          = ( p_hydrostatic(k) *                           &
-                                  ( 100000.0_wp / p_hydrostatic(k)             &
-                                  )**( r_d / c_p )                             &
-                                ) / ( r_d * pt_init(k) )
+          p_hydrostatic(k) = barometric_formula(zu(k), pt_surface *            & 
+                             exner_function(surface_pressure * 100.0_wp),      &
+                             surface_pressure * 100.0_wp)
+          
+          rho_air(k) = ideal_gas_law_rho_pt(p_hydrostatic(k), pt_init(k))
        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)  )
+                            
     ELSE
        DO  k = nzb, nzt+1
-          p_hydrostatic(k)    = surface_pressure * 100.0_wp *                  &
-                                ( 1 - ( g * zu(nzb) ) / ( c_p * t_surface )    &
-                                )**( c_p / r_d )
-          rho_air(k)          = ( p_hydrostatic(k) *                           &
-                                  ( 100000.0_wp / p_hydrostatic(k)             &
-                                  )**( r_d / c_p )                             &
-                                ) / ( r_d * pt_init(nzb) )
+          p_hydrostatic(k) = barometric_formula(zu(nzb), pt_surface *          & 
+                             exner_function(surface_pressure * 100.0_wp),      &
+                             surface_pressure * 100.0_wp)
+
+          rho_air(k) = ideal_gas_law_rho_pt(p_hydrostatic(k), 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
 !