Changeset 2252

Timestamp:

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

Author:

knoop

Message:

air density now depending on surface_pressure even in boussinesq mode

Location:

palm/trunk/SOURCE

Files:

• flow_statistics.f90 (modified) (5 diffs)
• init_3d_model.f90 (modified) (3 diffs)

palm/trunk/SOURCE/flow_statistics.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! perturbation pressure now depending on flux_output_mode
+! 
+! 2233 2017-05-30 18:08:54Z suehring
 !
 ! 2232 2017-05-30 17:47:52Z suehring
@@ -637 +640 @@
                 sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)  &
                                                               * flag
-                sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)  * flag
+                sums_l(k,40,tn) = sums_l(k,40,tn) + ( p(k,j,i)                 &
+                                         / momentumflux_output_conversion(k) ) &
+                                                              * flag
 
                 sums_l(k,33,tn) = sums_l(k,33,tn) + &
@@ -1415 +1420 @@
 
                    sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
-                                               * ( p(k,j,i) + p(k+1,j,i) ) * flag
+                                       * ( ( p(k,j,i) + p(k+1,j,i) )           &
+                                         / momentumflux_output_conversion(k) ) &
+                                       * flag
 
                 ENDDO
@@ -1750 +1757 @@
        hom(:,1,70,sr) = sums(:,70)     ! q*2
        hom(:,1,71,sr) = sums(:,71)     ! prho
-       hom(:,1,72,sr) = hyp * 1E-4_wp  ! hyp in dbar
+       hom(:,1,72,sr) = hyp * 1E-2_wp  ! hyp in hPa
        hom(:,1,73,sr) = sums(:,73)     ! nr
        hom(:,1,74,sr) = sums(:,74)     ! qr
@@ -1931 +1938 @@
        IF ( hom(k_surface_level,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )&
        THEN
-          hom(nzb+8,1,pr_palm,sr) = &
+          hom(nzb+8,1,pr_palm,sr) =                                            &
              ( g / hom(k_surface_level+1,1,4,sr) *                             &
-             ( hom(k_surface_level,1,18,sr) / heatflux_output_conversion(nzb) )&
+             ( hom(k_surface_level,1,18,sr) /                                  &
+             ( heatflux_output_conversion(nzb) * rho_air(nzb) ) )              &
              * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
        ELSE

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