code_original - Sean's code and sample output - 23/03/07

ancil - N96 HadGEM1 orography field - retrieved from NEC - 23/03/07

Missing data - REAL, PARAMETER  :: RMDI    = -9999.0 - 23/03/07

------------------------------------------------------------------------

Running the code - email from Sean, 23/03/07
------------------------------------------------------------------------

Ok, 3 files

        refrac_info.f90

is just a module containing various definitions.

        gpsro_driver.f90

is an example of a call to the 1d operator. I have
set the atmospheric parameters up in a simple way, so that it
should be pretty obvious what you have to extract from the UM
model output. Basically, you need the P, T, Q on the (full) model
levels, and the full model level heights.

        gpsro_op.f90

is the 1d bending angle operator. It is self contained,
so hopefully you should not have to go into it, assuming everything
works!


To compile and run

        f90 -c refrac_info.f90

this should create refrac_info.mod. Then type

        f90 gpsro_driver.f90 gpsro_op.f90 -o gpsro_op

Then

        gpsro_op

This will write a profile of bending angles to screen. My results are
given in alpha.out, so you can check its running ok.
-----------------------------------------------------------------------

UM height levels - from http://www-nwp.metoffice.com/~frcw/UM_Info/ND_levels/ND_levels.html  - 16/04/07
-----------------------------------------------------------------------------------------------------------
Blev (=pphead(52) ) contains Zsea  , the height of the level over the sea,
bhlev (=pphead(54)) contains C(k) = [ 1 - eta_value(k) / eta_value(first constant rho level) ] ^2,

so that the height for the level is  Z(i,j)= Zsea+ C(k) * orography(i,j) at each point (i,j)  in the grid.
 
------------------------------------------------------------------------------------------------------------

Tangent linear code - 03/07/07
--------------------------------------------------------------------------


I have attached the tangent linear code with a driver routine.
To compile and run the driver

f90  gpsro_drivertl.f90 gpsro_opTL.f90 gpsro_op.f90 -o gpsro_optl

You always need to compile gpsro_opTL.f90 with gpsro_op.f90 because they
share routines. The interface to the new code is

SUBROUTINE alpha_optl(nlev,   &   ! no. of model levs (=38)
                    nobs,   &   ! no. of bending angles in profile
                    roc,    &   ! radius of curv.
                    undul,  &   ! undulation (set to 0.0)
                    lat,    &   ! latitude of ob. location (degrees)
                    pres,   &   ! pressure on mod levels (hPa)
                    pres_prime, &
                    temp,   &   ! temp on model levels
                    temp_prime, &
                    q,      &   ! specific humidity (g/kg)
                    q_prime,    &
                    zg,     &   ! geopotential height of model levels(m)
                    a,      &   ! impact parameters
                    alpha,  &
                    alpha_prime)      ! bending angles

Basically, when you see a "_prime" that means a perturbation. The code
calculates the bending angles (alpha) as before with P,T and Q, but now 
also calculates the bending angle perturbations, alpha_prime, caused by 
the P, T and Q perturbations, p_prime, t_prime and q_prime, respectively.

The driver code calculates alpha_prime along with the differences 
between 2 calls to the normal operator for a simple profile. The results 
are very close.

The tangent linear code will enable you to look at the relative 
contributions of the p_prime, t_prime and q_prime to the alpha_prime.


------------------------------------------------------------------------------------------------------------

Further calculations - 15/02/08
--------------------------------------------------------------------------

1) set q to 0.0 just before refrac_levs is called and
look at the bending angle trends.

2) set q to 2000 value just before refrac_levs and
look at trends etc.

3) Look at refractivity trend on model levels. One
of the outputs from refrac_levs is refrac. This is
the refractivity on the model levels. (you could
also look at the q = 0 for refrac if you get time.)

4) Compare the Tdry and T trends on model levels.
You can calculate Tdry after calling refrac_levs

Tdry(:) =aval*pres(:)/refrac(:)

where aval is the refractivity constant defined in
refrac_info.

Look at error (Tdry(:)-temp(:)) and trend differences.