#DEFINITIONS:  -*-sh-*-

#This control file handles a DOIT batch calculation with ARTS internal
#pnd_field calculations. Six selected atmospheres from Chevallier91L
#data are used for atmospheric input. No sensor characteristics 
#are applied. Only two frequencies and two line of sights are used,
#for calculation time reasons.
#
# Author: Daniel Kreyling
# 

Arts2 {

INCLUDE "general/general.arts"
INCLUDE "general/continua.arts"
INCLUDE "general/agendas.arts"
INCLUDE "general/agendasDOIT.arts"
INCLUDE "general/planet_earth.arts"

# 1.General Settings:---------------------------------------------
#----------------------------------------------------------------- 

# Agenda for scalar gas absorption calculation
Copy(abs_xsec_agenda, abs_xsec_agenda__noCIA)

# (standard) emission calculation
Copy( iy_main_agenda, iy_main_agenda__Emission )

# cosmic background radiation
Copy( iy_space_agenda, iy_space_agenda__CosmicBackground )

# standard surface agenda (i.e., make use of surface_rtprop_agenda)
Copy( iy_surface_agenda, iy_surface_agenda__UseSurfaceRtprop )

# sensor-only path
Copy( ppath_agenda, ppath_agenda__FollowSensorLosPath )

# no refraction
Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath )


# Set out put file format
#------------------------
output_file_formatSetAscii

# Define f_grid
#--------------
VectorSet( f_grid, [9.0e10, 19.0e10])
#ReadXML(f_grid, "f_grid.xml")
#VectorNLinSpace (f_grid, 10, 120e9, 300e9)
#WriteXML (output_file_format, f_grid, "f_grid.xml")

#Set stokes dim
#--------------
IndexSet (stokes_dim, 1)

#def of atmosphere
#-----------------
IndexSet (atmosphere_dim, 1)

#No jacobian calculations
#-----------------
jacobianOff

# Modifiy the maximum propagation step, from the default(10.e3)
# to 250 m:---------------------------------------------------
NumericSet( ppath_lmax, 250 )

# Surface properties
#-------------------
# Set surface reflectivity (=1-emissivity)
# corresponds to emissivity=0.75
VectorSetConstant( surface_scalar_reflectivity, 1, 0.25 )

Copy( surface_rtprop_agenda,
      surface_rtprop_agenda__Specular_NoPol_ReflFix_SurfTFromt_surface )


# 2. Sensor:---------------------------------------------------------
#--------------------------------------------------------------------

# Definition of sensor position and LOS
# ------------------------------------
# Line of sight
MatrixSet( sensor_los,  [131; 179])
# Sensor position
nrowsGet( nrows, sensor_los )
ncolsGet( ncols, sensor_los )
MatrixSetConstant( sensor_pos, nrows, ncols, 850e3 )
# No sensor characteristics are specified
sensorOff


# 3. Read chevallier atmospheric profiles for batch calc--------------
#---------------------------------------------------------------------
ReadXML( batch_atm_fields_compact, "testdata/chevallierl91_all_extract.xml" )

# add constant profiles for O2 and N2
batch_atm_fields_compactAddConstant( name="abs_species-O2", value=0.2095 )
batch_atm_fields_compactAddConstant( name="abs_species-N2", value=0.7808 )


# 4. Absorption-------------------------------------------------
#---------------------------------------------------------------

abs_speciesSet( species=[ "H2O-PWR98",
                          "O3",
                          "O2-PWR93",
                          "N2-SelfContStandardType" ] )

# Creation of abs_lookup table
#-----------------------------
  # ATTENTION: The abs_lookup table used with this test control file was 
  #generated by the following code. It is adapted to the specified abs_species,
  #to the specified atmospheric data and frequencies ranging from 80e9 Hz to 
  #200e9 Hz. If changes to these inputs are applied the abs_lookup table must 
  #be recalculated.

  # Read HITRAN catalog (needed for O3):
  #abs_linesReadFromHitran2004( abs_lines,
  #                             "/storage3/data/catalogue/hitran/hitran2004/HITRAN04.par",
  #                             80e9,
  #                             200e9 )

  #abs_lines_per_speciesCreateFromLines

  #abs_lookupSetupBatch   	

  #abs_xsec_agenda_checkedCalc
  #abs_lookupCalc

  #WriteXML("ascii", abs_lookup, "abs_lookupBatch.xml")

# Reading of abs_lookup table
#----------------------------
ReadXML( abs_lookup, "abs_lookupBatch.xml" )

# Check if lookup table is fitting input
abs_lookupAdapt

# absorption from LUT
Copy( propmat_clearsky_agenda, propmat_clearsky_agenda__LookUpTable )


# 6.Specification of cloud-----------------------------------------------
# -----------------------------------------------------------------------
ScatSpeciesInit

# Definition of particle species
#--------------------------------
scat_speciesSet(scat_species, [
                "LWC-H98_STCO",
                "IWC-MH97",
                "RR-MP48"
])

# Read SingleScatteringData and ScatteringMetaData
#-------------------------------------------------
ArrayOfStringCreate( scat_data_files )
ReadXML( scat_data_files, "testdata/SingleScatteringFile_H2Oliquid.xml" )
ScatSpeciesScatAndMetaRead ( scat_data_files=scat_data_files )
ReadXML( scat_data_files, "testdata/SingleScatteringFile_H2Oice.xml" )
ScatSpeciesScatAndMetaRead ( scat_data_files=scat_data_files )
ReadXML( scat_data_files, "testdata/SingleScatteringFile_H2Oliquid.xml" )
ScatSpeciesScatAndMetaRead ( scat_data_files=scat_data_files )

# Selection/filterting of particles according to *scat_species*
#--------------------------------------------------------------
ScatElementsSelect( species="LWC",
                    sizeparam="diameter_volume_equ", sizemin=0.1e-6, sizemax=2000e-6 )
ScatElementsSelect( species="IWC",
                    sizeparam="diameter_volume_equ", sizemin=0.1e-6, sizemax=2000e-6 )

# Extend scattering data temperature range (allow ice at positive atm temps)
#--------------------------------------------------------------
ScatSpeciesExtendTemperature( species="IWC", T_high=350. )

scat_dataCheck( scat_data=scat_data_raw )
scat_dataCalc
scat_dataCheck

# 7.AGENDAS--------------------------------------------------------------
#------------------------------------------------------------------------

# Check the include files to see the setting of Agendas and if needed,
# overwrite them by re-setting the agendas here.------------------------


# Select zenith angle interpolation method ('linear' or 'polynomial'):
# --------------------------------------------------------------------
# If "polynomial" is selected one has to use an optimized grid. Please     
# use *doit_za_grid_optCalc* to optimize the grid.
doit_za_interpSet( doit_za_interp, atmosphere_dim, "linear" )

# Sets the angular grids for DOIT calculation     
# --------------------------------------------------------------------
# For down- and up-looking geometries. 
DOAngularGridsSet( doit_za_grid_size, scat_aa_grid, scat_za_grid,
      19, 37, "" )

# For limb calculations it is very important to have a finer zenith angle
# grid resolution about 90 degrees sza.
# The finer grid will be used for the RT calculations, while the scattering
# integral is still solved for the original sza grid.
# Make sure to set the path to "doit_za_grid_opt.xml" 
#   DOAngularGridsSet( doit_za_grid_size, scat_aa_grid, scat_za_grid,
#                     19, 10, "PATH/doit_za_grid_opt.xml" )

AgendaSet( doit_conv_test_agenda ){
  # Give limits for all Stokes components in BT:
  doit_conv_flagAbsBT( epsilon=[0.1],
                       max_iterations=100,
                       nonconv_return_nan=1)
}

StringCreate( info )
# Batch Agenda-----------------------------------------------------------
#------------------------------------------------------------------------
AgendaSet( ybatch_calc_agenda ){
  
  # Extract the atmospheric profiles for current atmosphere:
  Extract(
    atm_fields_compact, 
    batch_atm_fields_compact, 
    ybatch_index
  )

  # Split up *atm_fields_compact* to
  # generate p_grid, t_field, z_field, massdensity_field, vmr_field
  # Additionally moving TOA down to 10hPa altitude.
  AtmFieldsFromCompact( p_min=1e3 )

  # Unrealistic (small) values in scat_species_*_*_fields set to zero
  particle_fieldCleanup (scat_species_mass_density_field,
                             scat_species_mass_density_field,
                             1e-15)
  particle_fieldCleanup (scat_species_mass_flux_field,
                             scat_species_mass_flux_field,
                             1e-15)
  #WriteXMLIndexed (output_file_format, ybatch_index, scat_species_mass_density_field)

  # Get some surface properties from corresponding atmospheric fields
  Extract( z_surface, z_field, 0 )
  Extract( t_surface, t_field, 0 )

  # Set cloudbox limits automatically from scattering particle profiles 
  #(pressure units)
  # ---------------------------------------------------------------------
  cloudboxSetAutomatically( particle_field=scat_species_mass_density_field )
  cloudboxSetAutomatically( particle_field=scat_species_mass_flux_field,
                            cloudbox_limits_old=cloudbox_limits )
  cloudboxSetAutomatically( particle_field=scat_species_number_density_field,
                            cloudbox_limits_old=cloudbox_limits )
  cloudboxSetAutomatically( particle_field=scat_species_mean_mass_field,
                            cloudbox_limits_old=cloudbox_limits )

  # Alternative: set cloudbox manually
  # In this case the extent of the cloudbox should cover almost the 
  # whole troposphere, to make sure, that all scattering particles of all 
  # atmospheres are encompassed by it.
  #  cloudboxSetManually(cloudbox_on, cloudbox_limits, atmosphere_dim, 
  #      p_grid, lat_grid, lon_grid, 101300, 5000, 0, 0, 0, 0 )

  # Particle Number Density calculations
  # ------------------------------------
  pnd_fieldCalcFromscat_speciesFields
  
  #Write indexed pnd_fields to sub-folder for each atmosphere
  #WriteXMLIndexed (output_file_format, ybatch_index, pnd_field)

  # Consistency checks
  atmfields_checkedCalc( bad_partition_functions_ok = 1 )
  atmgeom_checkedCalc
  cloudbox_checkedCalc
  scat_data_checkedCalc
  sensor_checkedCalc

  # Initialize Doit variables
  DoitInit

  # Calculate incoming radiation field on cloudbox boundaries
  DoitGetIncoming

  # Set first guess field
  doit_i_fieldSetClearsky

  # Executes doit_mono_agenda for all frequencies
  DoitCalc

  # Here we choose our *y* output unit 
  StringSet( iy_unit, "PlanckBT" )

  # Calculate complete measurement vector 
  yCalc
}


# Set number of batch cases:
nelemGet( ybatch_n, batch_atm_fields_compact )
#IndexSet(ybatch_start, 0)
#IndexSet(ybatch_n, 4)

#===========start batch calc=====================
# Execute the batch calculations:
# This test control file can be run multi-threaded, since it was approved
# that none of the jobs fails. 
# If settings are changed, especially if the input atmospheres are altered 
# or exchanged, the robust option in *ybatchCalc* should be used.
propmat_clearsky_agenda_checkedCalc

ybatchCalc

# Call *ybatchCalc* robust:
# Set robust flag to 1. If one individual job fails, ARTS will continue
# with the next batch job.
# ybatchCalc(ybatch, ybatch_jacobians, ybatch_start, ybatch_n, 
#    ybatch_calc_agenda, 1 )

WriteXML( output_file_format, ybatch)

# Verify results
ArrayOfVectorCreate( ybatch_ref )
ReadXML( ybatch_ref, "TestDOITBatch.ybatch.ref.xml" )

Compare( ybatch, ybatch_ref, 0.01,
         "Total BT should be close to the reference values")

#==================stop==========================

} # End of Main