### computing pressure from PP

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.0000from 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?

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.

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

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

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

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

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

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

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

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

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

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

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

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

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

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.0000from 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?

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 + PPI must compute the reference pressure on full levels right?

Is this so far correct?

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 )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.

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