#
# Testing functionality of Regridding routines GriddedFieldLatLonRegrid and
#  GriddedFieldPRegrid
#
# Jana Mendrok 2013-02-26


Arts2 {

StringCreate( atmcase )
StringSet( atmcase, "planets/Earth/Fascod/tropical/tropical" )
StringCreate( surfpath )
StringSet( surfpath, "planets/Earth/ERA40/" )
StringCreate( zsurfname )
StringSet( zsurfname, "SurfaceAltitude_ERA40_1.0Degree" )
StringCreate( surfmaskname )
StringSet( surfmaskname, "LandSeaMask_ERA40_1.0Degree" )
StringCreate( Bname )
StringSet( Bname, "planets/Earth/IGRF/IGRF11_2010_200km-5deg-5deg" )
StringCreate( windname )
StringSet( windname, "planets/Mars/MPS/Mars.Ls0.day.dust-high/Mars.Ls0.day.dust-high.sol-avg/Mars.Ls0.day.dust-high.sol-avg" )

StringCreate( casefull )
StringCreate( caseext )
StringCreate( slat )
StringSet( slat, "lat:" )
StringCreate( slon )
StringSet( slon, "lon:" )

VectorNLogSpace( p_grid, 11, 1000e2, 0.1 )
VectorCreate( lat_regrid )
VectorLinSpace( lat_regrid, 45, -20, 10 )
VectorCreate( lon_regrid )
VectorLinSpace( lon_regrid, -3.5, 12.5, 0.25 )
NumericCreate( lattrue )
NumericCreate( lontrue )
IndexCreate( ncases )

IndexCreate( interp_order )
IndexSet( interp_order, 1 )

GriddedField2Create( z_surface_raw )
GriddedField2Create( lsmask_raw )
MatrixCreate( lsmask )
GriddedField3Create( B_field_raw )
GriddedField3Create( wind_field_raw )

MatrixCreate( mtmp )
Tensor3Create( t3tmp )

abs_speciesSet( species=["CH4", "H2O", "O2", "N2"] )
# get atm scenario raw data
AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw,
            abs_species, atmcase )

# for testing and demonstration purposes, we manually expand 1D atmosphere first
# to 3D then regrid back to 1D
GriddedFieldLatLonExpand( t_field_raw, t_field_raw )
GriddedFieldLatLonExpand( z_field_raw, z_field_raw )
GriddedFieldLatLonExpand( vmr_field_raw, vmr_field_raw )

# reading the surface altitude field
Copy( casefull, surfpath )
Append( casefull, zsurfname )
ReadXML( z_surface_raw, casefull )

# reading surface mask
Copy( casefull, surfpath )
Append( casefull, surfmaskname )
ReadXML( lsmask_raw, casefull )

# reading B-field w-component
Copy( casefull, Bname )
StringSet( caseext, ".B_w.xml" )
Append( casefull, caseext )
ReadXML( B_field_raw, casefull )

# reading wind-field w-component
Copy( casefull, windname )
StringSet( caseext, ".wind_w.xml" )
Append( casefull, caseext )
ReadXML( wind_field_raw, casefull )


AgendaCreate( forloop_agenda_lon )
AgendaSet( forloop_agenda_lon ){
  Extract( lontrue, lon_regrid, forloop_index )
  VectorSetConstant( lon_true, 1, lontrue )
  #Print( slon, 0 )
  #Print( lon_true, 0 )

  GriddedFieldLatLonRegrid( t_field_raw, lat_true, lon_true, t_field_raw )
  GriddedFieldPRegrid( t_field_raw, p_grid, t_field_raw, interp_order, 0 )
  FieldFromGriddedField( t_field, p_grid, lat_grid, lon_grid, t_field_raw )
  Extract( mtmp, t_field, 5 )
  #Print( mtmp, 0 )

  GriddedFieldLatLonRegrid( z_field_raw, lat_true, lon_true, z_field_raw )
  GriddedFieldPRegrid( z_field_raw, p_grid, z_field_raw, interp_order, 0 )
  FieldFromGriddedField( z_field, p_grid, lat_grid, lon_grid, z_field_raw )
  Extract( mtmp, z_field, 5 )
  #Print( mtmp, 0 )

  GriddedFieldLatLonRegrid( vmr_field_raw, lat_true, lon_true, vmr_field_raw )
  GriddedFieldPRegrid( vmr_field_raw, p_grid, vmr_field_raw, interp_order, 1 )
  FieldFromGriddedField( vmr_field, p_grid, lat_grid, lon_grid, vmr_field_raw )
  Extract( t3tmp, vmr_field, 0 )
  Extract( mtmp, t3tmp, 5 )
  #Print( mtmp, 0 )

  GriddedFieldLatLonRegrid( z_surface_raw, lat_true, lon_true, z_surface_raw )
  FieldFromGriddedField( z_surface, p_grid, lat_grid, lon_grid, z_surface_raw )
  #Print( z_surface, 0 )

  GriddedFieldLatLonRegrid( lsmask_raw, lat_true, lon_true, lsmask_raw )
  FieldFromGriddedField( lsmask, p_grid, lat_grid, lon_grid, lsmask_raw )
  #Print( lsmask, 0 )

  GriddedFieldPRegrid( B_field_raw, p_grid, B_field_raw, interp_order, 1 )
  GriddedFieldLatLonRegrid( B_field_raw, lat_true, lon_true, B_field_raw, interp_order )
  FieldFromGriddedField( mag_w_field, p_grid, lat_grid, lon_grid, B_field_raw )
  Extract( mtmp, mag_w_field, 5 )
  #Print( mtmp, 0 )

  # we can also do the expanding here...
  GriddedFieldPRegrid( wind_field_raw, p_grid, wind_field_raw, interp_order, 1 )
  GriddedFieldLatLonExpand(	wind_field_raw, wind_field_raw )
  GriddedFieldLatLonRegrid( wind_field_raw, lat_true, lon_true, wind_field_raw, interp_order )
  FieldFromGriddedField( wind_w_field, p_grid, lat_grid, lon_grid, wind_field_raw )
  Extract( mtmp, wind_w_field, 5 )
  #Print( mtmp, 0 )
}

AgendaCreate( forloop_agenda_lat )
AgendaSet( forloop_agenda_lat ){
  Extract( lattrue, lat_regrid, forloop_index )
  VectorSetConstant( lat_true, 1, lattrue )
  #Print( slat, 0 )
  #Print( lat_true, 0 )

  nelemGet( ncases, lon_regrid )
  IndexStepDown( ncases, ncases )
  Copy( forloop_agenda, forloop_agenda_lon )
  ForLoop( forloop_agenda, 0, ncases, 1 )
}

AtmosphereSet1D
nelemGet( ncases, lat_regrid )
IndexStepDown( ncases, ncases )
Copy( forloop_agenda, forloop_agenda_lat )
ForLoop( forloop_agenda, 0, ncases, 1 )

}