= Revising the way the height of atmospheric variables are specified in the driver =
Jan Polcher June 2012


== Current situation ==

The driver reads the height of the atmospheric forcing variables (Temperature, humidity and wind) from the run.def. These are the variables '''HEIGHT_LEV1''' (for T and Q) and '''HEIGHT_LEVW''' (for wind). This is done in routine forcing_grid in file readdim2.f90.

The scalar values obtained are then copied onto the global gird (i.e. iim x jjm in size).

When calling forcing_read these variables can be updated over the global grid. This is only the case when we are in a "WATCHOUT" situation. That is when when using output from an LMDZ_OR to drive ORCHIDEE off-line. 

Later, in the temporal loop of dim2_driver, if the level of the wind and temperature are not the same, the wind is interpolated vertically to the level of T and Q.

== Problems caused by this state of affairs ==

There are 3 problems with this approach :

  1. We impose the same value for the heigh of the atmospheric forcing variables to all points of the globe (except in the case of WATCHOUT).

  2. As the height of the forcing variable is in the run.def the user can change forcing data and forget to provide the new height of the variables. This can only be solved if the height of the forcing data is in the file containing the data itself.

  3. The height of the variables cannot change over time. But this is the case if the atmospheric variable come from a model which uses some type of pressure coordinate.

As the fluxes computed are strongly dependent on the height of the atmospheric variable a small error here as a large impact.

== Strategy for adapting the driver ==

The strategy to adapt the driver needs to ensure that what works today continues to work. It also needs to limit the work required to adapt existing netCDF forcing files.

We propose to add to the netCDF files a scalar information on the height of the 2 types of atmospheric variables (scalar : T, Q and vector : wind). This can be done easily with the following type of script :

Creat a cdf file (levels.cdf) of the following type :
{{{
netcdf filename {
variables:
        float HybSigA ;
              HybSigA:long_name = "hybrid level at layer midpoints" ;
              HybSigA:title="hyam (mlev=hyam+hybm*aps)";
        float HybSigB ;
              HybSigB:long_name = "hybrid level at layer midpoints" ;
              HybSigB:title="hybm (mlev=hyam+hybm*aps)";
data:
 HybSigA = 0.0 ;
 HybSigB = 0.998815059661865 ;
}
}}}

Then have a small script when generates the netCDF file (levels.nc) and adds the variables to an existing forcing file :
{{{
ncgen -o levels.nc levels.cdf
ncks -A -v HybSigA,HybSigB levels.nc My_ORCHIDEE_Forcing_2012.n
}}}

This scalar information can either be of 3 types :

  1. Sigma level height : This allows then to compute the actual height given the surface pressure.

  2. Hybrid sigma levels (as in the example above)  : Also given surface pressure we can compute the height.

  3. height in meters : This is the usual way we get forcing data.

Following a convention readdim2.f90 will check for the presence of either one of the 3 above listed cases. 

If one of the cases is found in the forcing file, forcing_just_read will compute the height directly after reading surface pressure.

If none of the above cases is found, the variables  '''HEIGHT_LEV1''' and '''HEIGHT_LEVW''' found in run.def will be used to fill the global fields.

== Proposed convention ==

We will test for the following variables in the forcing file :

=== Sigma and Sigma_uv ===

Expected is one real value which allows to compute the height of the first level with the following formula :

zlev_vec(i,j) = rau(i,j) * cte_grav * (psurf(i,j) - Sigma * psurf(i,j))

zlevuv_vec() = rau(i,j) * cte_grav * (psurf(i,j) - Sigma_uv * psurf(i,j))

Should only Sigma exist, Sigma=Sigma_uv will be assumed.

=== HybSigA, HybSigB, HybSigA_uv and HybSigB_uv ===

Expected is one real value for each variable. They allow to compute the height of the first level with the following formula :

zlev_vec(i,j) = rau(i,j) * cte_grav * (psurf(i,j) - (HybSigA + HybSigB * psurf(i,j)))

zlevuv_vec(i,j) = rau(i,j) * cte_grav * (psurf(i,j) - (HybSigA_uv + HybSigB_uv * psurf(i,j)))

If the UV values do not exist we assume that they are the same as for T and Q.

Sample CDF file is provided : [https://forge.ipsl.fr/orchidee/attachment/wiki/Branches/Driver_Atm_Lev/levels_Hybrid.cdf]

=== Levels and Levels_uv === 

This is the case when the levels are provied as 2D time evolving fields. The height is supposed to be in meters.

The following assignation will be done for the values read at each forcing time step :

zlev_vec(i,j) = Levles(i,j)

If levels_uv is present :

   zlevuv_vec(i,j) = levels_uv(i,j)

Else :

   zlevuv_vec(i,j) = levels(i,j)

=== Height_Lev1 and Height_Levuv ===

This comes back to the same case as when these variables are provided through the run.def file.

The following assignations will be done :

zlev_vec(i,j) = Height_Lev1

If Height_Levuv is present :

   zlevuv_vec(i,j) = Height_Levuv

Else :

   zlevuv_vec(i,j) = Height_Lev1

Sample CDF file is provided : [https://forge.ipsl.fr/orchidee/attachment/wiki/Branches/Driver_Atm_Lev/levels_Height.cdf]

== Proposed modifications in the code ==

=== dim2_driver.f90 ===

simplify the code by keeping only the following variables :
  * zlev_vec(:,:)
  * zlevuv_vec(:,:)
Their content will be provided by the call to '''forcing_read'''.
This makes redundant variables lev and levuv and they can be deleted from the arguments of '''forcing_grid'''.


The vertical interpolation would be done in all cases. Should zlev_vec and zlevuv_vec have the same value the interpolation does not change the values of u and v. This simplifies the code and makes the variable lower_wind redundant.

=== readdim2.f90, forcing_info ===

Will call forcing_vertical (before the call to init_data_para) to do most of the work  !

=== readdim2.f90, forcing_vertical ===

  * This routine will read '''HEIGHT_LEV1''' and '''HEIGHT_LEVW''' ... and use them if needed.

  * will explore the netCDF file to see what variables are available

  * Set the flags according to what needs to be done in order to compute zlev and zlevuv.

  * read the scalar variables giving the height.

The information gathered here will be saved in shared variables of the module.
 
=== readdim2.f90, forcing_grid ===

This routine will be simplified as it will not deal any more with the vertical layer. The levels will be considered as a time variable quantity.

=== readdim2.f90, forcing_read ===

weathgen_main will be expanded to also return the content of zlev and zlevuv.

=== readdim2.f90, forcing_read_interpol ===

The temporal interpolation of zlevuv needs to be added.

=== readdim2.f90, weathgen_main ===

=== readdim2.f90, forcing_just_read ===

Once surface pressure is read we compute or fill zlev and zlevuv

This will be done on the zoomed variables.

The "levels" variable in the is_watchout case will remain as it is. This should be homogenized but only be done in a second stage.

Comments in code need to be corrected.