# computing pressure from PP – in #9: CCLM

### 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.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,

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