#
# Testing functionality (meeting format requirements, etc.) of supplemental
#  atmospheric scenario data.
#
# General test setup: reading in raw data (including a basic atmosphere),
#  expanding (where necessary), regridding (incl. extracting), executing
#  standard pre-RT calc internal test method atmfields_checkedCalc.
#
#
# This case is for Venus and specifically tests
#
#  - for the five Venus scenarions: Venus.spicav.night, Venus.spicav.night_cold,
#     Venus.vira.night, Venus.vira.day, Venus.vira.day_highlat
#  - electron densities (given: 1D cases, 2 per night scneario, 5 per day scenario)
#      - all cases in 1D with p_grid from altitude grid
#      - a single case expanded to 3D (assuming the other cases behave in the
#        same way) with p_grid from altitude grid
#  - wind (given: 1D cases with both horizontal components, 1 case for NS
#     component and 3 cases for EW component per scenario)
#      - all expanded to global 3D case with p_grid from altitude grid
#
# Jana Mendrok 2013-02-26

Arts2 {
INCLUDE "general/general.arts"
INCLUDE "general/agendas.arts"
INCLUDE "general/planet_venus.arts"


# 1-dimensional atmosphere
AtmosphereSet1D

# We have to set the absorption lookup table interpolation order, 
# since we are using wind fields. atmferlds_checkedCalc will otherwise throw an error.
IndexSet( abs_f_interp_order, 1 )

Tensor3Create( edensity_field )
GriddedField3Create( gf3tmp )
StringCreate( caseext )
StringCreate( casefull )
StringCreate( atmcase )
IndexCreate( ncases )
ArrayOfStringCreate( atmcasearray )
ArrayOfStringCreate( necasearray )
IndexCreate( interp_order )

IndexSet( interp_order, 1 )

# set basic case folder
StringCreate( basename )
StringSet( basename, "planets/Venus/MPS/" )

# Array with case names
ArrayOfStringCreate( nightcasearray )
ArrayOfStringSet( nightcasearray,
  ["Venus.spicav.night", "Venus.spicav.night_cold", "Venus.vira.night"] )
ArrayOfStringCreate( nightnearray )
ArrayOfStringSet( nightnearray,
  [".SZA.90-100.Ne", ".SZA.100-120.Ne"] )

ArrayOfStringCreate( daycasearray )
ArrayOfStringSet( daycasearray,
  ["Venus.vira.day", "Venus.vira.day_highlat"] )
ArrayOfStringCreate( daynearray )
ArrayOfStringSet( daynearray,
  [".SZA.0-30.Ne", ".SZA.30-50.Ne", ".SZA.50-70.Ne", ".SZA.70-80.Ne", ".SZA.80-90.Ne"] )

ArrayOfStringCreate( windcasearray )
ArrayOfStringSet( windcasearray,
  ["_max", "_mid", "_min"] )

# the 3D geo grid to test
VectorCreate( lat_grid3D )
VectorCreate( lon_grid3D )
VectorLinSpace( lat_grid3D, -90, 90, 18 )
VectorLinSpace( lon_grid3D, -20, 350, 18 )

#####
# CASES A-C (night, 2 Ne profiles) and D-E (day, 5 Ne profiles)
#####

#we go with a foorloop through the different cases - night/day in an outer loop
# and the different Ne cases within in an inner loop.


AgendaCreate( forloop_agenda_Ne )
AgendaSet( forloop_agenda_Ne ){
  # construct the full case name
  Copy( casefull, atmcase )
  Extract( caseext, necasearray, forloop_index )
  Append( casefull, caseext )
  Print( casefull, 0 )
  # readin in raw data
  ReadXML( gf3tmp, casefull )

  #this is 1D data and we're doing 1D. but we need to regrid the raw data to
  # the calculation grid(s). for supplemental atm data (Ne, B, winds) this
  # requires manual regridding (in contrast to basic atm data, where this is
  # handled by AtmFieldsCalc.

  # so, first: regrid to p_grid (as we are in 1D, we don't need latlon regridding)
  GriddedFieldPRegrid( gf3tmp, p_grid, gf3tmp, interp_order, 1 )

  # eventually: extract the data Tensor from the regridded GriddedField
  FieldFromGriddedField( edensity_field, p_grid, lat_grid, lon_grid, gf3tmp )

  atmfields_checkedCalc
  atmgeom_checkedCalc
  #WriteXML( "ascii", p_grid )
  #WriteXML( "ascii", z_field )
  #WriteXML( "ascii", t_field )
  #WriteXMLIndexed( "ascii", forloop_index, edensity_field )
}


AgendaCreate( forloop_agenda_wind )
AgendaSet( forloop_agenda_wind ){
  # construct the full case name
  Copy( atmcase, casefull )
  Extract( caseext, windcasearray, forloop_index )
  Append( atmcase, caseext )
  Print( atmcase, 0 )
  # readin in raw data
  ReadXML( gf3tmp, atmcase )

  # first: regrid to p_grid
  GriddedFieldPRegrid( gf3tmp, p_grid, gf3tmp, interp_order, 1 )
  # second: make 3D fields from 1D, then regrid to lat/lon_grid
  GriddedFieldLatLonExpand(	gf3tmp, gf3tmp )
  GriddedFieldLatLonRegrid( gf3tmp, lat_grid, lon_grid, gf3tmp, interp_order )
  # last, extract the data Tensor from the regridded GriddedField
  FieldFromGriddedField( wind_v_field, p_grid, lat_grid, lon_grid, gf3tmp )

  atmfields_checkedCalc
  atmgeom_checkedCalc
}


AgendaCreate( forloop_agenda_scen )
AgendaSet( forloop_agenda_scen ){
  # construct atmscen name
  Copy( atmcase, basename )
  Extract( casefull, atmcasearray, forloop_index )
  Append( atmcase, casefull )
  StringSet( caseext, "/" )
  Append( atmcase, caseext )
  Append( atmcase, casefull )
  Print( atmcase, 0 )

  # we manually select a minumim set of basic atm data (main atm constituents)
  abs_speciesSet( species=["CO2"] )
  AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw,
              abs_species, atmcase )
  p_gridFromZRaw(	p_grid, z_field_raw, 0 )
  AtmFieldsCalc( vmr_zeropadding=1 )

  Extract( z_surface, z_field, 0 )

  # now get the Ne data. with several cases per scenario, we use another forloop
  Copy( forloop_agenda, forloop_agenda_Ne )
  nelemGet( ncases, necasearray )
  IndexStepDown( ncases, ncases )
  ForLoop( forloop_agenda, 0, ncases, 1  )

  # now changing to 3D for the winds (but do one Ne case per scenario as well)
  AtmosphereSet3D
  Copy( lat_grid, lat_grid3D )
  Copy( lon_grid, lon_grid3D )

  # blowing up basic atmosphere (and surface)
  AtmFieldsExpand1D
  Extract( z_surface, z_field, 0 )

  # getting and preprocessing Ne data in 3D (single case per scenario)
  Copy( casefull, atmcase )
  Extract( caseext, necasearray, 0 )
  Append( casefull, caseext )
  Print( casefull, 0 )
  # reading in raw data
  ReadXML( gf3tmp, casefull )
  GriddedFieldLatLonExpand(	gf3tmp, gf3tmp )
  GriddedFieldPRegrid( gf3tmp, p_grid, gf3tmp, interp_order, 1 )
  GriddedFieldLatLonRegrid( gf3tmp, lat_grid, lon_grid, gf3tmp, interp_order )
  FieldFromGriddedField( edensity_field, p_grid, lat_grid, lon_grid, gf3tmp )
  atmfields_checkedCalc
  atmgeom_checkedCalc

  # and eventually the winds - here we the two horizontal components (u&v), and
  # the u-component (zonal) in 3 different cases.

  #first the v-component (as this is not changing within one scenario)
  Copy( casefull, atmcase )
  StringSet( caseext, ".wind_v" )
  Append( casefull, caseext )
  ReadXML( gf3tmp, casefull )
  Print( casefull, 0 )
  # first: regrid to p_grid
  GriddedFieldPRegrid( gf3tmp, p_grid, gf3tmp, interp_order, 1 )
  # second: make 3D fields from 1D, then regrid to lat/lon_grid
  GriddedFieldLatLonExpand(	gf3tmp, gf3tmp )
  GriddedFieldLatLonRegrid( gf3tmp, lat_grid, lon_grid, gf3tmp, interp_order )
  # last, extract the data Tensor from the regridded GriddedField
  FieldFromGriddedField( wind_v_field, p_grid, lat_grid, lon_grid, gf3tmp )
  atmfields_checkedCalc
  atmgeom_checkedCalc

  #now the u-component. for the 3 cases per scenario, we again use a forloop.  
  Copy( casefull, atmcase )
  StringSet( caseext, ".wind_u" )
  Append( casefull, caseext )
  Copy( forloop_agenda, forloop_agenda_wind )
  nelemGet( ncases, windcasearray )
  IndexStepDown( ncases, ncases )
  ForLoop( forloop_agenda, 0, ncases, 1  )

}


Copy( forloop_agenda, forloop_agenda_scen )

Copy( atmcasearray, nightcasearray )
Copy( necasearray, nightnearray )
nelemGet( ncases, atmcasearray )
IndexStepDown( ncases, ncases )
ForLoop( forloop_agenda, 0, ncases, 1  )

Copy( atmcasearray, daycasearray )
Copy( necasearray, daynearray )
nelemGet( ncases, atmcasearray )
IndexStepDown( ncases, ncases )
ForLoop( forloop_agenda, 0, ncases, 1  )

}