# This controlfile demonstrates the calculation of Jacobians under cloudy # conditions using the scattering module DOIT. # # It is based on TestDoitBatch, but has been slightly adapted in certain # RT parameters. The Jacobian calculation control parameters are set in # section 2 (around l.130-160). # # Author: Jana Mendrok # Arts2{ # include general settings INCLUDE "general/general.arts" INCLUDE "general/continua.arts" INCLUDE "general/agendas.arts" INCLUDE "general/agendasDOIT.arts" INCLUDE "general/planet_earth.arts" # 1.General ARTS calculation Setup ------------------------------- #----------------------------------------------------------------- # Preparing a bunch of required agendas Copy(abs_xsec_agenda, abs_xsec_agenda__noCIA) Copy( iy_main_agenda, iy_main_agenda__Emission ) Copy( iy_space_agenda, iy_space_agenda__CosmicBackground ) Copy( iy_surface_agenda, iy_surface_agenda__UseSurfaceRtprop ) Copy( ppath_agenda, ppath_agenda__FollowSensorLosPath ) Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath ) output_file_formatSetAscii VectorSet( f_grid, [190.e9] ) IndexSet ( stokes_dim, 1 ) AtmosphereSet1D StringSet( iy_unit, "PlanckBT" ) VectorSetConstant( surface_scalar_reflectivity, 1, 0.4) Copy( surface_rtprop_agenda, surface_rtprop_agenda__Specular_NoPol_ReflFix_SurfTFromt_field ) MatrixSet( sensor_los, [180.]) nrowsGet( nrows, sensor_los ) ncolsGet( ncols, sensor_los ) MatrixSetConstant( sensor_pos, nrows, ncols, 600e3 ) sensorOff ReadXML( batch_atm_fields_compact, "testdata/chevallierl91_all_extract.xml" ) batch_atm_fields_compactAddConstant( name="abs_species-O2", value=0.2095 ) batch_atm_fields_compactAddConstant( name="abs_species-N2", value=0.7808 ) abs_speciesSet( species=[ "H2O-PWR98", "O3", "O2-PWR93", "N2-SelfContStandardType" ] ) ReadXML( abs_lookup, "artscomponents/doitbatch/abs_lookupBatch.xml" ) abs_lookupAdapt Copy( propmat_clearsky_agenda, propmat_clearsky_agenda__LookUpTable ) ScatSpeciesInit scat_speciesSet(scat_tags=[ "LWC-MGD_LWC", "IWC-MGD_IWC", "RR-MP48" # "SR-F07TR" ]) ArrayOfStringCreate( scat_data_files ) ReadXML( scat_data_files, "scattering/H2O_liquid/SingleScatteringFile_allH2Oliquid.xml" ) ScatSpeciesScatAndMetaRead ( scat_data_files=scat_data_files ) ReadXML( scat_data_files, "scattering/H2O_ice/SingleScatteringFile_allH2Oice.xml" ) ScatSpeciesScatAndMetaRead ( scat_data_files=scat_data_files ) ReadXML( scat_data_files, "scattering/H2O_liquid/SingleScatteringFile_allH2Oliquid.xml" ) ScatSpeciesScatAndMetaRead ( scat_data_files=scat_data_files ) 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 ) #ScatElementsSelect( species="RR", sizeparam="diameter_volume_equ" ) #ScatElementsSelect( species="SR", sizeparam="diameter_volume_equ" ) scat_dataCalc doit_za_interpSet( doit_za_interp, atmosphere_dim, "linear" ) DOAngularGridsSet( N_za_grid=19, N_aa_grid=37 ) #----------------------------------------------------------------- # 2. Doit Jacobian specific settings ----------------------------- #----------------------------------------------------------------- # Convergence test agenda # ----------------------- # Prepare some container variables holding the Doit Jacobian control parameters VectorCreate( doit_epsilon_ref ) VectorCreate( doit_epsilon_perturb ) IndexCreate( n_iterations ) ArrayOfIndexCreate( cloudlimit_levels ) # convergence agenda for reference case AgendaCreate( doit_conv_test_agenda_refcase ) AgendaSet( doit_conv_test_agenda_refcase ){ doit_conv_flagAbsBT( epsilon=doit_epsilon_ref ) #Print(doit_iteration_counter,0) #WriteXMLIndexed(in=doit_i_field,file_index=doit_iteration_counter) } # convergence agenda for perturbation runs AgendaCreate( doit_conv_test_agenda_perturb ) AgendaSet( doit_conv_test_agenda_perturb ){ doit_conv_flagAbsBT( epsilon=doit_epsilon_perturb, max_iterations=n_iterations ) #Print(doit_iteration_counter,0) } # The main DOIT agenda # -------------------- AgendaSet( doit_mono_agenda ){ DoitScatteringDataPrepare Ignore( f_grid ) doit_i_field_monoIterate } ###################################################################################### # Set Jacobian parameters. # ----------------------- # Set epsilon for reference run of DOIT (give limits for all Stokes components) VectorSet( doit_epsilon_ref, [0.001] ) # Epsilon for perturbation runs of DOIT - here we want n_iterations as break # criterion, hence we set epsilon very small here VectorScale( doit_epsilon_perturb, doit_epsilon_ref, 1e-6 ) # Set fixed number of iteration for perturbation runs IndexSet( n_iterations, 23 ) # Set the extend of the cloudbox, i.e. of the Jacobians (ATTENTION: all non-zero # regions of scattering species fields need to be within cloudbox! (you can # switch off the fixed cloudbox limit setting and switch on automatic cloudbox # limits below. then setting here will have no effect). ArrayOfIndexSet( cloudlimit_levels, [0,32] ) ###################################################################################### # 3. The RT and Jacobian calculation ----------------------------- #----------------------------------------------------------------- ArrayOfArrayOfSingleScatteringDataCreate( scat_data_orig ) ArrayOfArrayOfScatteringMetaDataCreate( scat_meta_orig ) ArrayOfStringCreate( scat_species_orig ) Tensor4Create( pnd_field_orig ) AgendaSet( ybatch_calc_agenda ){ Touch( y_aux ) # Extract compact atmosphere from the batch of atmospheric states Extract( atm_fields_compact, batch_atm_fields_compact, ybatch_index ) # Split up compact atmosphere into atmospheric fields AtmFieldsFromCompact # Remove unrealistic (small) values from scat_species_*_*_fields particle_fieldCleanup (scat_species_mass_density_field, scat_species_mass_density_field, 1e-7) particle_fieldCleanup (scat_species_mass_flux_field, scat_species_mass_flux_field, 1e-7) # Get some surface properties from corresponding atmospheric fields Extract( z_surface, z_field, 0 ) # Set cloudbox limits # --- # here, we prefer the manual setting with a resetting to a constant range of # atm levels (to make all cases having roughly the same vertical extend) #cloudboxSetAutomatically cloudboxSetManually(p1=11e4, p2=1e-1, lat1=0, lat2=0, lon1=0, lon2=0 ) Copy(cloudbox_limits, cloudlimit_levels) ###################################################################################### # Case specific Jacobian settings # ------------------------------- # if you want to have DOIT jacobians, this need to be done, but BEFORE the # jacobianDoitAddSpecies! # it's used here also the (re)initialize the jacobian WSV for each batch # calculation (for example, retrieval grids might be different). jacobianOff # DOIT Jacobians to calculate #jacobianDoitAddSpecies( species="T", unit="abs", dx=1. ) #jacobianDoitAddSpecies( species="abs_species.H2O-PWR98", unit="abs", dx=1e-6 ) #jacobianDoitAddSpecies( species="scat_species.LWC.mass_density", unit="abs", dx=1e-6 ) #jacobianDoitAddSpecies( species="scat_species.IWC.mass_density", unit="abs", dx=1e-6 ) jacobianDoitAddSpecies( species="scat_species.RR.mass_flux", unit="abs", dx=1e-6 ) jacobianDoitClose #WriteXMLIndexed( in=jacobian_quantities, file_index=ybatch_index ) #WriteXMLIndexed( in=jacobian_indices, file_index=ybatch_index ) ###################################################################################### # Particle Number Density calculations pnd_fieldCalcFromscat_speciesFields # Consistency checks atmgeom_checkedCalc atmfields_checkedCalc( bad_partition_functions_ok = 1 ) cloudbox_checkedCalc scat_data_checkedCalc sensor_checkedCalc propmat_clearsky_agenda_checkedCalc # when doing no ScatSpeciesMerge (neither here, nor in # jacobianDoit!), comment out this complete block Copy( scat_data_orig, scat_data ) Copy( scat_meta_orig, scat_meta ) Copy( scat_species_orig, scat_species ) Copy( pnd_field_orig, pnd_field ) ScatSpeciesMerge cloudbox_checkedCalc scat_data_checkedCalc # Initialize Doit variables DoitInit # Calculate incoming radiation field on cloudbox boundaries DoitGetIncoming(rigorous=0) # Set first guess field: doit_i_fieldSetClearsky # Initialization convergence agenda for reference RT case (for Jacobians) Copy( doit_conv_test_agenda, doit_conv_test_agenda_refcase ) # Executes doit_mono_agenda for all frequencies DoitCalc #WriteXMLIndexed(in=doit_i_field,file_index=ybatch_index) # Re-setting agendas and WSV containers for perturbation runs Copy( doit_conv_test_agenda, doit_conv_test_agenda_perturb ) # Comment that out if ScatSpeciesMerge not used (neither here, nor in # jacobianDoit!) Copy( scat_data, scat_data_orig ) Copy( scat_meta, scat_meta_orig ) Copy( scat_species, scat_species_orig ) Copy( pnd_field, pnd_field_orig ) jacobianDoit( ScatSpeciesMerge_do=1 ) } # for batch cases: nelemGet( ybatch_n, batch_atm_fields_compact ) IndexSet(ybatch_start, 5) IndexSet(ybatch_n, 1) ybatchCalc #( robust=1 ) # Store results: # --- # for non-batch cases #WriteXML( "zascii", y ) #WriteXML( "zascii", jacobian ) #WriteXML( "ascii", jacobian_quantities ) #WriteXML( "ascii", jacobian_indices ) # for batch cases #WriteXML( "zascii", ybatch ) #WriteXML( "zascii", ybatch_jacobians ) # note: jacobian_quantities and jacobian_indices are not part of ybatchCalc # output and are (or can be) specific to each batch case. hence, they should be # written out within the ybatch agenda, e.g. using WriteXMLIndexed. # Verify results # --- ArrayOfVectorCreate( ybatch_ref ) ReadXML( ybatch_ref, "ybatchREFERENCE_DoitJacobians.xml" ) Compare( ybatch, ybatch_ref, 1e-3 ) ArrayOfMatrixCreate( ybatch_jacobians_ref ) ReadXML( ybatch_jacobians_ref, "ybatch_jacobiansREFERENCE_DoitJacobians.xml" ) Compare( ybatch_jacobians, ybatch_jacobians_ref, 1e-4 ) }