#DEFINITIONS:  -*-sh-*-
#
# filename: TestDISORT.arts
#
# Demonstration of a DISORT scattering calculation
#
# The setup follows the one used by TestDOITFromIndividualFields.arts (in doit/).
#
# Author: Jana Mendrok
# 

Arts2 {

INCLUDE "artscomponents/disort/indivfieldsatmo_setup.arts"

#####
# Reference output calculation
#  - For reference, a DOIT scattering calculation is used.
#  - Switch this block on, if changes were made to indivfieldsatmo_setup.arts,
#    hence a reference recalculation is required.
#  - Note 1: The DOIT reference solution will take  significantly longer 
#    computation time than DISORT.
#  - Note 2: A recalculation will overwrite the precalculated reference output.
#####
#Copy( surface_rtprop_agenda,
#      surface_rtprop_agenda__Blackbody_SurfTFromt_surface )
#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 )
#pnd_fieldCalcFromscat_speciesFields
##cloudbox_checkedCalc  
##ScatSpeciesMerge
##cloudbox_checkedCalc
#INCLUDE "general/agendasDOIT.arts"
#doit_za_interpSet( doit_za_interp, atmosphere_dim, "linear" )
#DOAngularGridsSet( doit_za_grid_size, scat_aa_grid, scat_za_grid,
#      37, 73, "" )
#AgendaSet( doit_conv_test_agenda ){
#  doit_conv_flagAbsBT( epsilon=[0.1],
#                       max_iterations=100,
#                       nonconv_return_nan=1) }
#atmfields_checkedCalc
#atmgeom_checkedCalc
#cloudbox_checkedCalc  
#sensor_checkedCalc
#DoitInit
#DoitGetIncoming( rigorous=0 )
#doit_i_fieldSetClearsky
#DoitCalc
#yCalc
#WriteXML( output_file_format, y, "yREFERENCEfromDOIT.y.xml" )
#####
# End of reference calculation block
#####



# Set surface skin temperature based on t_field
Extract( t_surface, t_field, 0 )
VectorCreate( tskin_vector )
VectorExtractFromMatrix( tskin_vector, t_surface, 0, "column" )
Extract( surface_skin_t, tskin_vector, 0 )


# Set cloudbox limits 
# ---------------------------------------------------------------------
#cloudboxSetFullATm
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 )

# Particle Number Density field calculation
# -----------------------------------------
pnd_fieldCalcFromscat_speciesFields



# 7. DISORT settings-----------------------------------------------------
#------------------------------------------------------------------------

# For DISORT we use the same angular grids as output.
DOAngularGridsSet( N_za_grid=38 )

# Calculate opticle properties of particles and add particle absorption
# and extiction to the gaseous properties to get total extinction and
# absorption:

AgendaSet( iy_cloudbox_agenda ){
  iyInterpCloudboxField
}


# 8. The RT calculation--------------------------------------------------
#------------------------------------------------------------------------

# Consistency checks
atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc  
sensor_checkedCalc
scat_data_checkedCalc

# Perform scattering calculation 
#DisortCalc #( nstreams=64 ) #( non_iso_inc=1 )
DisortCalc( new_optprop=1 )
#WriteXML( output_file_format, doit_i_field )
   
# Calculate RT from cloudbox boundary to the sensor
yCalc

#WriteXML( output_file_format, y )

# Verify results
VectorCreate(yREFERENCE)
ReadXML( yREFERENCE, "yREFERENCEfromDOIT.y.xml" )
Compare( y, yREFERENCE, 1. )

ReadXML( yREFERENCE, "yREFERENCE.y.xml" )
Compare( y, yREFERENCE, 1e-6 )

} # End of Main