computing pressure from PP – in #9: CCLM

in #9: CCLM

Hi,

I have my 3d (+time = 4d) variable PP (deviation from reference pressure).
Now I simple want to compute my pressure on the full levels.

from YUSPECIF :

Vertical coordinate: type 2 (Gal-Chen hybrid) k Z(k) P0(k) 1 25000.0000 28.1726 2 23800.0000 33.9524 […..] 60 10.0000 998.8149 61 0.0000 1000.0000 Boundary for change from z to terrain-following: vcflat = 11430.0000 Half-level index where levels become flat: k = 15 Exponential profile with asymptotic isothermal stratosphere Values of the Reference Atmosphere: Pressure at Sea-Level: 1000.0000 Temperature at Sea-Level: 15.0000 Temp. diff. sea-level-stratosph.: 75.0000 Scale height for temp. decrease: 10000.0000

from Output constant file (lffd2015123118c.nc)

vcoord:long_name = “Height-based hybrid Gal-Chen coordinate” ; vcoord:units = “Pa” ; vcoord:ivctype = 2 ; vcoord:irefatm = 2 ; vcoord:p0sl = 100000. ; vcoord:t0sl = 288.149993896484 ; vcoord:dt0lp = 42. ; vcoord:vcflat = 11430. ; vcoord:delta_t = 75. ; vcoord:h_scal = 10000. ;

As I understand it, I need to compute the reference pressure for every grid point and then just add PP?
I must compute the reference pressure on full levels right?

Step 1: H = ( HLL [:,:,1:60] + HLL [:,:,2:61] ) * 0.5 Step 2: P_ref = specialFunction(H) = “pressure in Height H“ Step 3: P = P_ref + PP

Is this so far correct?
In that case what is this “specialFunction”? (I know it is basically barometric height formula… but which version?)
I tried the formula from the COSMO Documentation (Part I: Dynamics and Numerics) page 28/29, equation (3.85)
p_0(z) = … for both cases beta = 0 and beta = 42
but if I test it with the values from YUSPECIF I get lower or higher values than there.

nc = nc_open(file) vcoord = ncvar_get(nc,“vcoord”) p0sl = ncatt_get(nc,“vcoord”,“p0sl”)$value t0sl = ncatt_get(nc,“vcoord”,“t0sl”)$value dt0lp = ncatt_get(nc,“vcoord”,“dt0lp”)$value pref1 = p0sl * exp( -t0sl/dt0lp * (1-sqrt(1 – 2 * dt0lp * 9.81 * vcoord / 287.058 / t0sl^2) ) pref2 = p0sl * exp( -9.81 * vcoord / 287.058 / t0sl )

for 25000m
pref1 gives 13 hPa
pref2 gives 51 hPa
while YUSPECIF says 28 hPa

Any suggestions?
Rolf

P.S.: the attribute:units of the variable “vcoord” says “Pa” but its obviously(?) meter

edit: I attached the files

  @rolfzentek in #fcc580e

Hi,

I have my 3d (+time = 4d) variable PP (deviation from reference pressure).
Now I simple want to compute my pressure on the full levels.

from YUSPECIF :

Vertical coordinate: type 2 (Gal-Chen hybrid) k Z(k) P0(k) 1 25000.0000 28.1726 2 23800.0000 33.9524 […..] 60 10.0000 998.8149 61 0.0000 1000.0000 Boundary for change from z to terrain-following: vcflat = 11430.0000 Half-level index where levels become flat: k = 15 Exponential profile with asymptotic isothermal stratosphere Values of the Reference Atmosphere: Pressure at Sea-Level: 1000.0000 Temperature at Sea-Level: 15.0000 Temp. diff. sea-level-stratosph.: 75.0000 Scale height for temp. decrease: 10000.0000

from Output constant file (lffd2015123118c.nc)

vcoord:long_name = “Height-based hybrid Gal-Chen coordinate” ; vcoord:units = “Pa” ; vcoord:ivctype = 2 ; vcoord:irefatm = 2 ; vcoord:p0sl = 100000. ; vcoord:t0sl = 288.149993896484 ; vcoord:dt0lp = 42. ; vcoord:vcflat = 11430. ; vcoord:delta_t = 75. ; vcoord:h_scal = 10000. ;

As I understand it, I need to compute the reference pressure for every grid point and then just add PP?
I must compute the reference pressure on full levels right?

Step 1: H = ( HLL [:,:,1:60] + HLL [:,:,2:61] ) * 0.5 Step 2: P_ref = specialFunction(H) = “pressure in Height H“ Step 3: P = P_ref + PP

Is this so far correct?
In that case what is this “specialFunction”? (I know it is basically barometric height formula… but which version?)
I tried the formula from the COSMO Documentation (Part I: Dynamics and Numerics) page 28/29, equation (3.85)
p_0(z) = … for both cases beta = 0 and beta = 42
but if I test it with the values from YUSPECIF I get lower or higher values than there.

nc = nc_open(file) vcoord = ncvar_get(nc,“vcoord”) p0sl = ncatt_get(nc,“vcoord”,“p0sl”)$value t0sl = ncatt_get(nc,“vcoord”,“t0sl”)$value dt0lp = ncatt_get(nc,“vcoord”,“dt0lp”)$value pref1 = p0sl * exp( -t0sl/dt0lp * (1-sqrt(1 – 2 * dt0lp * 9.81 * vcoord / 287.058 / t0sl^2) ) pref2 = p0sl * exp( -9.81 * vcoord / 287.058 / t0sl )

for 25000m
pref1 gives 13 hPa
pref2 gives 51 hPa
while YUSPECIF says 28 hPa

Any suggestions?
Rolf

P.S.: the attribute:units of the variable “vcoord” says “Pa” but its obviously(?) meter

edit: I attached the files

computing pressure from PP

RZ wrote

Hi,

I have my 3d (+time = 4d) variable PP (deviation from reference pressure).
Now I simple want to compute my pressure on the full levels.

from YUSPECIF :

Vertical coordinate: type 2 (Gal-Chen hybrid) k Z(k) P0(k) 1 25000.0000 28.1726 2 23800.0000 33.9524 […..] 60 10.0000 998.8149 61 0.0000 1000.0000 Boundary for change from z to terrain-following: vcflat = 11430.0000 Half-level index where levels become flat: k = 15 Exponential profile with asymptotic isothermal stratosphere Values of the Reference Atmosphere: Pressure at Sea-Level: 1000.0000 Temperature at Sea-Level: 15.0000 Temp. diff. sea-level-stratosph.: 75.0000 Scale height for temp. decrease: 10000.0000

from Output constant file (lffd2015123118c.nc)

vcoord:long_name = “Height-based hybrid Gal-Chen coordinate” ; vcoord:units = “Pa” ; vcoord:ivctype = 2 ; vcoord:irefatm = 2 ; vcoord:p0sl = 100000. ; vcoord:t0sl = 288.149993896484 ; vcoord:dt0lp = 42. ; vcoord:vcflat = 11430. ; vcoord:delta_t = 75. ; vcoord:h_scal = 10000. ;

As I understand it, I need to compute the reference pressure for every grid point and then just add PP?
I must compute the reference pressure on full levels right?

Step 1: H = ( HLL [:,:,1:60] + HLL [:,:,2:61] ) * 0.5 Step 2: P_ref = specialFunction(H) = “pressure in Height H“ Step 3: P = P_ref + PP

Is this so far correct?
In that case what is this “specialFunction”? (I know it is basically barometric height formula… but which version?)
I tried the formula from the COSMO Documentation (Part I: Dynamics and Numerics) page 28/29, equation (3.85)
p_0(z) = … for both cases beta = 0 and beta = 42
but if I test it with the values from YUSPECIF I get lower or higher values than there.

nc = nc_open(file) vcoord = ncvar_get(nc,“vcoord”) p0sl = ncatt_get(nc,“vcoord”,“p0sl”)$value t0sl = ncatt_get(nc,“vcoord”,“t0sl”)$value dt0lp = ncatt_get(nc,“vcoord”,“dt0lp”)$value pref1 = p0sl * exp( -t0sl/dt0lp * (1-sqrt(1 – 2 * dt0lp * 9.81 * vcoord / 287.058 / t0sl^2) ) pref2 = p0sl * exp( -9.81 * vcoord / 287.058 / t0sl )

for 25000m
pref1 gives 13 hPa
pref2 gives 51 hPa
while YUSPECIF says 28 hPa

Any suggestions?
Rolf

P.S.: the attribute:units of the variable “vcoord” says “Pa” but its obviously(?) meter

edit: I attached the files

View in channel

Hi Rolf,

I send you an excel table for the case of 40 levels, which was prepared by H-J and myself in order to analyse the inconsistencies between the grid definitions in COSMO . It shows all related quantities of the reference atmosphere and the same values as in YUSPECIF .
You may adapt it for 60 levels and compare the formula with yours. This should give the correct values.

I’m sorry, but I do not have the time at the moment to check it on my own.

Regards, Andreas

Rolf Zentek wrote:
> Hi,
>
> I have my 3d (+time = 4d) variable PP (deviation from reference pressure).
> Now I simple want to compute my pressure on the full levels.
>
> from YUSPECIF:
>
> Vertical coordinate: type 2 (Gal-Chen hybrid)
> k Z(k) P0(k)
> 1 25000.0000 28.1726
> 2 23800.0000 33.9524
> […..]
> 60 10.0000 998.8149
> 61 0.0000 1000.0000
> Boundary for change from z to terrain-following: vcflat = 11430.0000
> Half-level index where levels become flat: k = 15
> Exponential profile with asymptotic isothermal stratosphere
> Values of the Reference Atmosphere:
> Pressure at Sea-Level: 1000.0000
> Temperature at Sea-Level: 15.0000
> Temp. diff. sea-level-stratosph.: 75.0000
> Scale height for temp. decrease: 10000.0000
>
> from Output constant file (lffd2015123118c.nc)
>
> vcoord:long_name = “Height-based hybrid Gal-Chen coordinate” ;
> vcoord:units = “Pa” ;
> vcoord:ivctype = 2 ;
> vcoord:irefatm = 2 ;
> vcoord:p0sl = 100000. ;
> vcoord:t0sl = 288.149993896484 ;
> vcoord:dt0lp = 42. ;
> vcoord:vcflat = 11430. ;
> vcoord:delta_t = 75. ;
> vcoord:h_scal = 10000. ;
>
> As I understand it, I need to compute the reference pressure for every grid point and then just add PP?
> I must compute the reference pressure on full levels right?
>
> Step 1: H = ( HLL[:,:,1:60] + HLL[:,:,2:61] ) * 0.5
> Step 2: P_ref = specialFunction(H) = “pressure in Height H“
> Step 3: P = P_ref + PP
>
> Is this so far correct?
> In that case what is this “specialFunction”? (I know it is basically barometric height formula… but which version?)
> I tried the formula from the COSMO Documentation (Part I: Dynamics and Numerics) page 28/29, equation (3.85)
> p_0(z) = … for both cases beta = 0 and beta = 42
> but if I test it with the values from YUSPECIF I get lower or higher values than there.
>
> nc = nc_open(file)
> vcoord = ncvar_get(nc,“vcoord”)
> p0sl = ncatt_get(nc,“vcoord”,“p0sl”)$value
> t0sl = ncatt_get(nc,“vcoord”,“t0sl”)$value
> dt0lp = ncatt_get(nc,“vcoord”,“dt0lp”)$value
> pref1 = p0sl * exp( -t0sl/dt0lp * (1-sqrt(1 – 2 * dt0lp * 9.81 * vcoord / 287.058 / t0sl^2) )
> pref2 = p0sl * exp( -9.81 * vcoord / 287.058 / t0sl )
>
> for 25000m
> pref1 gives 13 hPa
> pref2 gives 51 hPa
> while YUSPECIF says 28 hPa
>
> Any suggestions?
> Rolf
>
> P.S.: the attribute:units of the variable “vcoord” says “Pa” but its obviously(?) meter
>
> edit: I attached the files

  @andreaswill in #c06172b

Hi Rolf,

I send you an excel table for the case of 40 levels, which was prepared by H-J and myself in order to analyse the inconsistencies between the grid definitions in COSMO . It shows all related quantities of the reference atmosphere and the same values as in YUSPECIF .
You may adapt it for 60 levels and compare the formula with yours. This should give the correct values.

I’m sorry, but I do not have the time at the moment to check it on my own.

Regards, Andreas

Rolf Zentek wrote:
> Hi,
>
> I have my 3d (+time = 4d) variable PP (deviation from reference pressure).
> Now I simple want to compute my pressure on the full levels.
>
> from YUSPECIF:
>
> Vertical coordinate: type 2 (Gal-Chen hybrid)
> k Z(k) P0(k)
> 1 25000.0000 28.1726
> 2 23800.0000 33.9524
> […..]
> 60 10.0000 998.8149
> 61 0.0000 1000.0000
> Boundary for change from z to terrain-following: vcflat = 11430.0000
> Half-level index where levels become flat: k = 15
> Exponential profile with asymptotic isothermal stratosphere
> Values of the Reference Atmosphere:
> Pressure at Sea-Level: 1000.0000
> Temperature at Sea-Level: 15.0000
> Temp. diff. sea-level-stratosph.: 75.0000
> Scale height for temp. decrease: 10000.0000
>
> from Output constant file (lffd2015123118c.nc)
>
> vcoord:long_name = “Height-based hybrid Gal-Chen coordinate” ;
> vcoord:units = “Pa” ;
> vcoord:ivctype = 2 ;
> vcoord:irefatm = 2 ;
> vcoord:p0sl = 100000. ;
> vcoord:t0sl = 288.149993896484 ;
> vcoord:dt0lp = 42. ;
> vcoord:vcflat = 11430. ;
> vcoord:delta_t = 75. ;
> vcoord:h_scal = 10000. ;
>
> As I understand it, I need to compute the reference pressure for every grid point and then just add PP?
> I must compute the reference pressure on full levels right?
>
> Step 1: H = ( HLL[:,:,1:60] + HLL[:,:,2:61] ) * 0.5
> Step 2: P_ref = specialFunction(H) = “pressure in Height H“
> Step 3: P = P_ref + PP
>
> Is this so far correct?
> In that case what is this “specialFunction”? (I know it is basically barometric height formula… but which version?)
> I tried the formula from the COSMO Documentation (Part I: Dynamics and Numerics) page 28/29, equation (3.85)
> p_0(z) = … for both cases beta = 0 and beta = 42
> but if I test it with the values from YUSPECIF I get lower or higher values than there.
>
> nc = nc_open(file)
> vcoord = ncvar_get(nc,“vcoord”)
> p0sl = ncatt_get(nc,“vcoord”,“p0sl”)$value
> t0sl = ncatt_get(nc,“vcoord”,“t0sl”)$value
> dt0lp = ncatt_get(nc,“vcoord”,“dt0lp”)$value
> pref1 = p0sl * exp( -t0sl/dt0lp * (1-sqrt(1 – 2 * dt0lp * 9.81 * vcoord / 287.058 / t0sl^2) )
> pref2 = p0sl * exp( -9.81 * vcoord / 287.058 / t0sl )
>
> for 25000m
> pref1 gives 13 hPa
> pref2 gives 51 hPa
> while YUSPECIF says 28 hPa
>
> Any suggestions?
> Rolf
>
> P.S.: the attribute:units of the variable “vcoord” says “Pa” but its obviously(?) meter
>
> edit: I attached the files

AW replied

Hi Rolf,

I send you an excel table for the case of 40 levels, which was prepared by H-J and myself in order to analyse the inconsistencies between the grid definitions in COSMO . It shows all related quantities of the reference atmosphere and the same values as in YUSPECIF .
You may adapt it for 60 levels and compare the formula with yours. This should give the correct values.

I’m sorry, but I do not have the time at the moment to check it on my own.

Regards, Andreas

Rolf Zentek wrote:
> Hi,
>
> I have my 3d (+time = 4d) variable PP (deviation from reference pressure).
> Now I simple want to compute my pressure on the full levels.
>
> from YUSPECIF:
>
> Vertical coordinate: type 2 (Gal-Chen hybrid)
> k Z(k) P0(k)
> 1 25000.0000 28.1726
> 2 23800.0000 33.9524
> […..]
> 60 10.0000 998.8149
> 61 0.0000 1000.0000
> Boundary for change from z to terrain-following: vcflat = 11430.0000
> Half-level index where levels become flat: k = 15
> Exponential profile with asymptotic isothermal stratosphere
> Values of the Reference Atmosphere:
> Pressure at Sea-Level: 1000.0000
> Temperature at Sea-Level: 15.0000
> Temp. diff. sea-level-stratosph.: 75.0000
> Scale height for temp. decrease: 10000.0000
>
> from Output constant file (lffd2015123118c.nc)
>
> vcoord:long_name = “Height-based hybrid Gal-Chen coordinate” ;
> vcoord:units = “Pa” ;
> vcoord:ivctype = 2 ;
> vcoord:irefatm = 2 ;
> vcoord:p0sl = 100000. ;
> vcoord:t0sl = 288.149993896484 ;
> vcoord:dt0lp = 42. ;
> vcoord:vcflat = 11430. ;
> vcoord:delta_t = 75. ;
> vcoord:h_scal = 10000. ;
>
> As I understand it, I need to compute the reference pressure for every grid point and then just add PP?
> I must compute the reference pressure on full levels right?
>
> Step 1: H = ( HLL[:,:,1:60] + HLL[:,:,2:61] ) * 0.5
> Step 2: P_ref = specialFunction(H) = “pressure in Height H“
> Step 3: P = P_ref + PP
>
> Is this so far correct?
> In that case what is this “specialFunction”? (I know it is basically barometric height formula… but which version?)
> I tried the formula from the COSMO Documentation (Part I: Dynamics and Numerics) page 28/29, equation (3.85)
> p_0(z) = … for both cases beta = 0 and beta = 42
> but if I test it with the values from YUSPECIF I get lower or higher values than there.
>
> nc = nc_open(file)
> vcoord = ncvar_get(nc,“vcoord”)
> p0sl = ncatt_get(nc,“vcoord”,“p0sl”)$value
> t0sl = ncatt_get(nc,“vcoord”,“t0sl”)$value
> dt0lp = ncatt_get(nc,“vcoord”,“dt0lp”)$value
> pref1 = p0sl * exp( -t0sl/dt0lp * (1-sqrt(1 – 2 * dt0lp * 9.81 * vcoord / 287.058 / t0sl^2) )
> pref2 = p0sl * exp( -9.81 * vcoord / 287.058 / t0sl )
>
> for 25000m
> pref1 gives 13 hPa
> pref2 gives 51 hPa
> while YUSPECIF says 28 hPa
>
> Any suggestions?
> Rolf
>
> P.S.: the attribute:units of the variable “vcoord” says “Pa” but its obviously(?) meter
>
> edit: I attached the files

View in channel

Hey Andreas,

I took the first sheet (ivctype_2_Ref_Atmos_alt)
and just changed the Value in B5 to 25000 resulting in 13 hPa for J5 (p0_HalfLevel) which uses the same formula as pref1 from my post.
But its unequal to the Value from YUSPECIF .

regards
Rolf

  @rolfzentek in #ddc1938

Hey Andreas,

I took the first sheet (ivctype_2_Ref_Atmos_alt)
and just changed the Value in B5 to 25000 resulting in 13 hPa for J5 (p0_HalfLevel) which uses the same formula as pref1 from my post.
But its unequal to the Value from YUSPECIF .

regards
Rolf

RZ replied

Hey Andreas,

I took the first sheet (ivctype_2_Ref_Atmos_alt)
and just changed the Value in B5 to 25000 resulting in 13 hPa for J5 (p0_HalfLevel) which uses the same formula as pref1 from my post.
But its unequal to the Value from YUSPECIF .

regards
Rolf

View in channel

Hey,

does anyone know if vgrid_refatm_utils.f90 is the right file to look at?

The SUBROUTINE reference_atmosphere seems to include both formulas I cited from the documentation.

The SUBROUTINE reference_atmosphere_2 seems to include another formulas, which computes approximatly(*) the same values as my YUSPECIF .

(*) I suppose due to different values of constants?!

Here in R-style: 

    nc = nc_open(file)
    vcoord = ncvar_get(nc,"vcoord")
    p0sl = ncatt_get(nc,"vcoord","p0sl")$value
    t0sl = ncatt_get(nc,"vcoord","t0sl")$value
    dt0lp = ncatt_get(nc,"vcoord","dt0lp")$value
    delta_t = ncatt_get(nc,"vcoord","delta_t")$value
    h_scal = ncatt_get(nc,"vcoord","h_scal")$value
    g = 9.80665#9.807
    R = 287.05
    pref1 = p0sl * exp( -t0sl/dt0lp * (1-sqrt(1 - 2 * dt0lp * g * vcoord / R / t0sl^2) ) )
    pref2 = p0sl * exp( -g * vcoord / R / t0sl )
    t00 = t0sl - delta_t
    pref3 = p0sl * exp( -g*h_scal/ R/t00 * log( (exp(vcoord/h_scal)*t00 + delta_t)/(t00 + delta_t)) )

Can anyone confirm, that for my configuration I need to use pref3 ? (or that indeed reference_atmosphere_2 was used)

Rolf

  @rolfzentek in #c73b43f

Hey,

does anyone know if vgrid_refatm_utils.f90 is the right file to look at?

The SUBROUTINE reference_atmosphere seems to include both formulas I cited from the documentation.

The SUBROUTINE reference_atmosphere_2 seems to include another formulas, which computes approximatly(*) the same values as my YUSPECIF .

(*) I suppose due to different values of constants?!

Here in R-style: 

    nc = nc_open(file)
    vcoord = ncvar_get(nc,"vcoord")
    p0sl = ncatt_get(nc,"vcoord","p0sl")$value
    t0sl = ncatt_get(nc,"vcoord","t0sl")$value
    dt0lp = ncatt_get(nc,"vcoord","dt0lp")$value
    delta_t = ncatt_get(nc,"vcoord","delta_t")$value
    h_scal = ncatt_get(nc,"vcoord","h_scal")$value
    g = 9.80665#9.807
    R = 287.05
    pref1 = p0sl * exp( -t0sl/dt0lp * (1-sqrt(1 - 2 * dt0lp * g * vcoord / R / t0sl^2) ) )
    pref2 = p0sl * exp( -g * vcoord / R / t0sl )
    t00 = t0sl - delta_t
    pref3 = p0sl * exp( -g*h_scal/ R/t00 * log( (exp(vcoord/h_scal)*t00 + delta_t)/(t00 + delta_t)) )

Can anyone confirm, that for my configuration I need to use pref3 ? (or that indeed reference_atmosphere_2 was used)

Rolf

RZ replied

Hey,

does anyone know if vgrid_refatm_utils.f90 is the right file to look at?

The SUBROUTINE reference_atmosphere seems to include both formulas I cited from the documentation.

The SUBROUTINE reference_atmosphere_2 seems to include another formulas, which computes approximatly(*) the same values as my YUSPECIF .

(*) I suppose due to different values of constants?!

Here in R-style: 

    nc = nc_open(file)
    vcoord = ncvar_get(nc,"vcoord")
    p0sl = ncatt_get(nc,"vcoord","p0sl")$value
    t0sl = ncatt_get(nc,"vcoord","t0sl")$value
    dt0lp = ncatt_get(nc,"vcoord","dt0lp")$value
    delta_t = ncatt_get(nc,"vcoord","delta_t")$value
    h_scal = ncatt_get(nc,"vcoord","h_scal")$value
    g = 9.80665#9.807
    R = 287.05
    pref1 = p0sl * exp( -t0sl/dt0lp * (1-sqrt(1 - 2 * dt0lp * g * vcoord / R / t0sl^2) ) )
    pref2 = p0sl * exp( -g * vcoord / R / t0sl )
    t00 = t0sl - delta_t
    pref3 = p0sl * exp( -g*h_scal/ R/t00 * log( (exp(vcoord/h_scal)*t00 + delta_t)/(t00 + delta_t)) )

Can anyone confirm, that for my configuration I need to use pref3 ? (or that indeed reference_atmosphere_2 was used)

Rolf

View in channel

update / solution:

For everyone having the same problem: my INT2LM setting is irefatm = 2

From looking into the source code it seems to me that that:
reference_atmosphere_2 and thus pref3 is used because of my INT2LM setting irefatm = 2

I think the formula from the excel sheet and from the documentation are for the case irefatm = 1

From the documentation I get the same as from the source code: 

The new reference atmosphere is based on the temperature profile t0(z) = (t0sl-delta_t) + delta_t*exp(-z/h_scal), where z = hhl(k) is the height of a model grid point.

From the source code vgrid_refatm_utils.f90 / reference_atmosphere_2 : 
The base-state pressure at main levels is computed from the analytically integrated hydrostatic equation, assuming that the height at full levels (zhfl) is the arithmetic mean of the adjacent half levels.

which I interprete to be the equation pref3 with (hhl[,,1:60]+hhl[,,2:61])/2 instead of vcoord
When I did the analytical integration of the hydrostatic equation myself I got a different formula than pref3, but they give the same values :D so in short:

My opionion: If irefatm = 2 then p_reference(z) = p0sl * exp( -g*h_scal/ R/t00 * log( (exp(z/h_scal)*t00 + delta_t)/(t00 + delta_t)) )
with the values from my earlier posts

  @rolfzentek in #2b9ff43

update / solution:

For everyone having the same problem: my INT2LM setting is irefatm = 2

From looking into the source code it seems to me that that:
reference_atmosphere_2 and thus pref3 is used because of my INT2LM setting irefatm = 2

I think the formula from the excel sheet and from the documentation are for the case irefatm = 1

From the documentation I get the same as from the source code: 

The new reference atmosphere is based on the temperature profile t0(z) = (t0sl-delta_t) + delta_t*exp(-z/h_scal), where z = hhl(k) is the height of a model grid point.

From the source code vgrid_refatm_utils.f90 / reference_atmosphere_2 : 
The base-state pressure at main levels is computed from the analytically integrated hydrostatic equation, assuming that the height at full levels (zhfl) is the arithmetic mean of the adjacent half levels.

which I interprete to be the equation pref3 with (hhl[,,1:60]+hhl[,,2:61])/2 instead of vcoord
When I did the analytical integration of the hydrostatic equation myself I got a different formula than pref3, but they give the same values :D so in short:

My opionion: If irefatm = 2 then p_reference(z) = p0sl * exp( -g*h_scal/ R/t00 * log( (exp(z/h_scal)*t00 + delta_t)/(t00 + delta_t)) )
with the values from my earlier posts

RZ replied

update / solution:

For everyone having the same problem: my INT2LM setting is irefatm = 2

From looking into the source code it seems to me that that:
reference_atmosphere_2 and thus pref3 is used because of my INT2LM setting irefatm = 2

I think the formula from the excel sheet and from the documentation are for the case irefatm = 1

From the documentation I get the same as from the source code: 

The new reference atmosphere is based on the temperature profile t0(z) = (t0sl-delta_t) + delta_t*exp(-z/h_scal), where z = hhl(k) is the height of a model grid point.

From the source code vgrid_refatm_utils.f90 / reference_atmosphere_2 : 
The base-state pressure at main levels is computed from the analytically integrated hydrostatic equation, assuming that the height at full levels (zhfl) is the arithmetic mean of the adjacent half levels.

which I interprete to be the equation pref3 with (hhl[,,1:60]+hhl[,,2:61])/2 instead of vcoord
When I did the analytical integration of the hydrostatic equation myself I got a different formula than pref3, but they give the same values :D so in short:

My opionion: If irefatm = 2 then p_reference(z) = p0sl * exp( -g*h_scal/ R/t00 * log( (exp(z/h_scal)*t00 + delta_t)/(t00 + delta_t)) )
with the values from my earlier posts

View in channel