#
# 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 Mars and specifically tests
#
#  - for the 72 Mars scenarios: 4seasons x 2daytimes x 3dustloads x 3solaractivities
#  - electron densities (given: 1D cases, 1 per night scenario, 4 per day scenario)
#      - all cases in 1D with p_grid from altitude grid
#      - a single case per scenario expanded to 3D (assuming the other cases
#        behave in the same way) with p_grid from altitude grid
#  - wind (given: one 1D case per scenario with all 3 components)
#      - 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_mars.arts"


# 1-dimensional atmosphere
AtmosphereSet1D

# We have to set the absorption lookup table interpolation order, 
# since we are using wind fields. atmfields_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 )
IndexCreate( interp_order )

IndexSet( interp_order, 1 )


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

# Arrays with (sub)case names
ArrayOfStringCreate( seasoncasearray )
ArrayOfStringSet( seasoncasearray, ["Mars.Ls0", "Mars.Ls90", "Mars.Ls180", "Mars.Ls270"] )
ArrayOfStringCreate( timecasearray )
ArrayOfStringSet( timecasearray, [".day", ".night"] )
ArrayOfStringCreate( dustcasearray )
ArrayOfStringSet( dustcasearray, [".dust-high", ".dust-low", ".dust-medium"] )
ArrayOfStringCreate( solarcasearray )
ArrayOfStringSet( solarcasearray, [".sol-avg", ".sol-max", ".sol-min"] )

ArrayOfStringCreate( necasearray )
ArrayOfStringCreate( nightnearray )
ArrayOfStringSet( nightnearray, [".SZA.120-180.Ne"] )
ArrayOfStringCreate( daynearray )
ArrayOfStringSet( daynearray,
 [".SZA.0-30.Ne", ".SZA.30-50.Ne", ".SZA.50-70.Ne", ".SZA.70-90.Ne"] )

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


# we go with several nested foorloop through the different cases.
#  All those cases have identical abs species to process.
#  Order of agenda definitions has to be inverse from their execution (as we can
#  only copy an agenda AFTER we have defined it).

AgendaCreate( forloop_agenda_Ne )
AgendaSet( forloop_agenda_Ne ){
  # construct the full case name
  Copy( casefull, basename )
  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_solar )
AgendaSet( forloop_agenda_solar ){
  # construct atmcase name IV (Mars.LsXX.YY.dust-ZZ.sol-WW)
  Extract( casefull, solarcasearray, forloop_index )
  Append( atmcase, casefull )
  Append( basename, atmcase )
  StringSet( caseext, "/" )
  Append( basename, caseext )
  Append( basename, atmcase )
  #Print( basename, 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, basename )
  p_gridFromZRaw(	p_grid, z_field_raw, 0 )
  AtmFieldsCalc

  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, basename )
  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 w-component
  Copy( casefull, basename )
  StringSet( caseext, ".wind_w" )
  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_w_field, p_grid, lat_grid, lon_grid, gf3tmp )

  #sencond the u-component
  Copy( casefull, basename )
  StringSet( caseext, ".wind_u" )
  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_u_field, p_grid, lat_grid, lon_grid, gf3tmp )

  #last the v-component
  Copy( casefull, basename )
  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
}

AgendaCreate( forloop_agenda_dust )
AgendaSet( forloop_agenda_dust ){
  # construct atmcase name III (Mars.LsXX.YY.dust-ZZ)
  Extract( casefull, dustcasearray, forloop_index )
  Append( atmcase, casefull )
  # keep the casestring till dust and make upper-level folder name
  Append( basename, atmcase )
  StringSet( caseext, "/" )
  Append( basename, caseext )
  Copy( forloop_agenda, forloop_agenda_solar )
  nelemGet( ncases, solarcasearray )
  IndexStepDown( ncases, ncases )
  ForLoop( forloop_agenda, 0, ncases, 1 )
}

AgendaCreate( forloop_agenda_season )
AgendaSet( forloop_agenda_season ){
  #since we have different (number of) cases (Ne) depending on day/night, we
  # can't use a loop to process them, but need to set and call each separately.
  # hence we replace the timeloop (which would loops over 2 instances) by
  # explicit calls.

  #first: day
  # construct atmcase name I (Mars.LsXX)
  Extract( atmcase, seasoncasearray, forloop_index )
  # construct atmcase name II (Mars.LsXX.d/n)
  Extract( casefull, timecasearray, 0 )
  Append( atmcase, casefull )
  # set to use the correct Ne case array
  Copy( necasearray, daynearray )
  Copy( forloop_agenda, forloop_agenda_dust )
  nelemGet( ncases, dustcasearray )
  IndexStepDown( ncases, ncases )
  ForLoop( forloop_agenda, 0, ncases, 1 )

  #second: night
  # construct atmcase name I (Mars.LsXX)
  Extract( atmcase, seasoncasearray, forloop_index )
  # construct atmcase name II (Mars.LsXX.d/n)
  Extract( casefull, timecasearray, 1 )
  Append( atmcase, casefull )
  # set to use the correct Ne case array
  Copy( necasearray, nightnearray )
  Copy( forloop_agenda, forloop_agenda_dust )
  nelemGet( ncases, dustcasearray )
  IndexStepDown( ncases, ncases )
  ForLoop( forloop_agenda, 0, ncases, 1 )
}


nelemGet( ncases, seasoncasearray )
IndexStepDown( ncases, ncases )
Copy( forloop_agenda, forloop_agenda_season )
ForLoop( forloop_agenda, 0, ncases, 1  )

}