Ignore:
Timestamp:
Jul 11, 2019 1:57:56 PM (2 years ago)
Author:
Giersch
Message:

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/init_3d_model.f90

    r4048 r4088  
    2525! -----------------
    2626! $Id$
     27! Pressure and density profile calculations revised using basic functions
     28!
     29! 4048 2019-06-21 21:00:21Z knoop
    2730! Further modularization of particle code components
    2831!
     
    845848    ALLOCATE( drho_air_zw(nzb:nzt+1) )
    846849!
    847 !-- Density profile calculation for anelastic approximation
    848     t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**( r_d / c_p )
    849     IF ( TRIM( approximation ) == 'anelastic' ) THEN
     850!-- Density profile calculation for anelastic and Boussinesq approximation
     851!-- In case of a Boussinesq approximation, a constant density is calculated
     852!-- mainly for output purposes. This density do not need to be considered
     853!-- in the model's system of equations.
     854!    t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**( r_d / c_p )
     855    IF ( TRIM( approximation ) == 'anelastic' )  THEN
    850856       DO  k = nzb, nzt+1
    851           p_hydrostatic(k)    = surface_pressure * 100.0_wp *                  &
    852                                 ( 1 - ( g * zu(k) ) / ( c_p * t_surface )      &
    853                                 )**( c_p / r_d )
    854           rho_air(k)          = ( p_hydrostatic(k) *                           &
    855                                   ( 100000.0_wp / p_hydrostatic(k)             &
    856                                   )**( r_d / c_p )                             &
    857                                 ) / ( r_d * pt_init(k) )
     857          p_hydrostatic(k) = barometric_formula(zu(k), pt_surface *            &
     858                             exner_function(surface_pressure * 100.0_wp),      &
     859                             surface_pressure * 100.0_wp)
     860         
     861          rho_air(k) = ideal_gas_law_rho_pt(p_hydrostatic(k), pt_init(k))
    858862       ENDDO
     863       
    859864       DO  k = nzb, nzt
    860865          rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
    861866       ENDDO
     867       
    862868       rho_air_zw(nzt+1)  = rho_air_zw(nzt)                                    &
    863869                            + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
     870                           
    864871    ELSE
    865872       DO  k = nzb, nzt+1
    866           p_hydrostatic(k)    = surface_pressure * 100.0_wp *                  &
    867                                 ( 1 - ( g * zu(nzb) ) / ( c_p * t_surface )    &
    868                                 )**( c_p / r_d )
    869           rho_air(k)          = ( p_hydrostatic(k) *                           &
    870                                   ( 100000.0_wp / p_hydrostatic(k)             &
    871                                   )**( r_d / c_p )                             &
    872                                 ) / ( r_d * pt_init(nzb) )
     873          p_hydrostatic(k) = barometric_formula(zu(nzb), pt_surface *          &
     874                             exner_function(surface_pressure * 100.0_wp),      &
     875                             surface_pressure * 100.0_wp)
     876
     877          rho_air(k) = ideal_gas_law_rho_pt(p_hydrostatic(k), pt_init(nzb))
    873878       ENDDO
     879       
    874880       DO  k = nzb, nzt
    875881          rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
    876882       ENDDO
     883       
    877884       rho_air_zw(nzt+1)  = rho_air_zw(nzt)                                    &
    878885                            + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
     886                           
    879887    ENDIF
    880888!
Note: See TracChangeset for help on using the changeset viewer.