#DEFINITIONS: -*-sh-*- # # This is a test of weighting function calculations under cloudy conditions # using iyHybrid with background field from RT4 calculations. # # The setup is as far as possible equivalent to TestDoitJacobians.arts in the # same folder, but using yCalc with iyHybrid in the iy_main_agenda in place of # the jacobianDoit WSM. # # 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 ------------------------------- #----------------------------------------------------------------- VectorCreate( yref ) # 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 ) AgendaCreate( iy_main_agenda__Hybrid ) AgendaSet( iy_main_agenda__Hybrid ){ Ignore( nlte_field ) Ignore( iy_id ) ppathCalc( cloudbox_on = 0 ) iyHybrid Touch(iy_aux) } 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.] ) #; 150.]) nrowsGet( nrows, sensor_los ) ncolsGet( ncols, sensor_los ) MatrixSetConstant( sensor_pos, nrows, ncols, 600e3 ) #VectorNLinSpace( sensor_time, nrows, 0, 1 ) VectorSetConstant( sensor_time, nrows, 0. ) sensorOff # some stuff for temperature perturbations NumericSet( p_hse, 1000e2 ) NumericSet( z_hse_accuracy, 1 ) VectorSet( lat_true, [15] ) VectorSet( lon_true, [123] ) #ReadXML( batch_atm_fields_compact, "chevallierl91_all_extract_RWC.xml" ) 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-MH97", "RR-W16" # "SR-F07TR" ]) # Microphysics # scat species #0 ArrayOfStringSet( pnd_agenda_input_names, [ "IWC" ] ) # parameters of this PSD. # can be multiple (T not to # appear here; is passed # anyways). Append( pnd_agenda_array_input_names, pnd_agenda_input_names ) ArrayOfAgendaAppend( pnd_agenda_array ){ ScatSpeciesSizeMassInfo( species_index=agenda_array_index, x_unit="dveq" ) Copy( psd_size_grid, scat_species_x ) Copy( pnd_size_grid, scat_species_x ) psdMH97 #( t_min=258 ) # t-limits modified for consistency with # pnd_fieldCalcFromscat_speciesFields interface pndFromPsd( scat_index=agenda_array_index ) #pndAdjustFromScatMeta( scat_index=agenda_array_index, pnd_agenda_input_tag="RWC" ) } # scat species #1 ArrayOfStringSet( pnd_agenda_input_names, [ "RWC" ] ) # parameters of this PSD. # can be multiple (T not to # appear here; is passed # anyways). Append( pnd_agenda_array_input_names, pnd_agenda_input_names ) ArrayOfAgendaAppend( pnd_agenda_array ){ ScatSpeciesSizeMassInfo( species_index=agenda_array_index, x_unit="dveq" ) Copy( psd_size_grid, scat_species_x ) Copy( pnd_size_grid, scat_species_x ) psdW16( t_min=258 ) # t-limits modified for consistency with # pnd_fieldCalcFromscat_speciesFields interface # to W16, which has no t-limits at all (but scat_data # allow T down to 210K at most. hence, t_min=0 doesn't # work either.) # semi-sensible setting that is still consistent with # pnd_fieldCalcFromscat_speciesFields: # allowing supercooled rain drops down to -15C pndFromPsd( scat_index=agenda_array_index ) #pndAdjustFromScatMeta( scat_index=agenda_array_index, pnd_agenda_input_tag="RWC" ) } 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_dataCheck( scat_data=scat_data_raw ) scat_dataCalc #scat_dataCheck( sca_mat_threshold=8e-3 ) scat_data_checkedCalc doit_za_interpSet( doit_za_interp, atmosphere_dim, "linear" ) #----------------------------------------------------------------- # 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 (DOIT) # ----------------------- # 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] ) ###################################################################################### StringCreate( info ) StringCreate( scat_name ) Tensor3Create( scat_field ) ArrayOfTensor4Create( dpnd_dummy ) # 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 # for now AtmFieldsFromCompact extracts into scat_species_*_*_field(s) # hence, to get the data into particle_bulkprop_field, we need to do a bit of # massaging... # (what exactly needs to be done depends a bit on what is in the compact atm # data, and therefore need to be manually adjusted for your specific atm # compact data :-/ ) Extract( scat_field, scat_species_mass_density_field, 0 ) Append( particle_bulkprop_field, scat_field ) StringSet( scat_name, "IWC" ) Append( particle_bulkprop_names, scat_name ) #Extract( scat_field, scat_species_mass_density_field, 1 ) Extract( scat_field, scat_species_mass_flux_field, 1 ) Append( particle_bulkprop_field, scat_field ) StringSet( scat_name, "RWC" ) Append( particle_bulkprop_names, scat_name ) # 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) particle_fieldCleanup (particle_bulkprop_field, particle_bulkprop_field, 1e-7) # Get some surface properties from corresponding atmospheric fields Extract( z_surface, z_field, 0 ) ##### # iyHybrid ##### # iyHybrid calc (here instead of after jacDoit for debugging (don't want to # wait jacDoit to finish each time...). Also, doing it here we don't have to # deal with troubles from ScatSpeciesMerge. jacobianInit jacobianAddScatSpecies( g1=p_grid, g2=lat_grid, g3=lon_grid, species="RR-W16", quantity="RWC" ) jacobianAddScatSpecies( g1=p_grid, g2=lat_grid, g3=lon_grid, species="IWC-MH97", quantity="IWC" ) #jacobianAddWind( g1=p_grid, g2=lat_grid, g3=lon_grid, # component="v" ) #jacobianAddTemperature( g1=p_grid, g2=lat_grid, g3=lon_grid, # hse="off", method="analytical", dt=0.1 ) ### not working simultaneously with analytical jacs so far (assert fail) #jacobianAddTemperature( g1=p_grid, g2=lat_grid, g3=lon_grid, # hse="off", method="perturbation", dt=0.1 ) #jacobianAddWind( g1=p_grid, g2=lat_grid, g3=lon_grid, # component="w" ) jacobianAddAbsSpecies( g1=p_grid, g2=lat_grid, g3=lon_grid, species="H2O-PWR98", unit="vmr" ) ### not working simultaneously with analytical jacs so far (assert fail) #jacobianAddAbsSpecies( g1=p_grid, g2=lat_grid, g3=lon_grid, # species="O3", # method="perturbation", unit="rel", dx=0.01 ) ### not working so far (assert fail) #jacobianAddFreqShift( df=50e3 ) #jacobianAddFreqStretch( df=50e3 ) #jacobianAddPolyfit( poly_order=1 ) #jacobianAddSinefit( period_lengths=[200e6,40e6] ) jacobianClose cloudboxSetFullAtm #pnd_fieldCalcFromscat_speciesFields pnd_fieldCalcFromParticleBulkProps #WriteXML( in=t_field ) #WriteXML( in=z_field ) #WriteXML( in=vmr_field ) #WriteXML( in=particle_bulkprop_field, filename="BulkProp.particle_bulkprop_field.xml" ) #WriteXML( in=pnd_field, filename="BulkProp.pnd_field.xml" ) #WriteXML( in=dpnd_field_dx, filename="BulkProp.dpnd_field_dx.xml" ) atmgeom_checkedCalc atmfields_checkedCalc( bad_partition_functions_ok = 1 ) cloudbox_checkedCalc sensor_checkedCalc propmat_clearsky_agenda_checkedCalc StringSet( info, "RT4Calc" ) Print( info, 0 ) RT4Calc #( pfct_aa_grid_size=28, nstreams=36 ) # sufficient to make Z-norm accurate # enough for monodisperse 1cm raindrop #WriteXMLIndexed( in=doit_i_field, file_index=ybatch_index, filename="RT4.doit_i_field" ) #WriteXML( in=scat_za_grid, filename="RT4.scat_za_grid.xml" ) Copy( iy_main_agenda, iy_main_agenda__Hybrid ) StringSet( info, "yCalc with iyHybrid" ) Print( info, 0 ) yCalc #WriteXMLIndexed( in=jacobian, file_index=ybatch_index, filename="RT4.jacobian" ) Copy( yref, y ) ##### # DOIT ##### # DOAngularGridsSet( N_aa_grid=19, N_za_grid=37 ) # # # 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="rel", dx=1e-3 ) # #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_density", unit="rel", dx=1e-3 ) # jacobianDoitClose # #WriteXMLIndexed( in=jacobian_quantities, file_index=ybatch_index ) # #WriteXMLIndexed( in=jacobian_indices, file_index=ybatch_index ) ####################################################################################### # # # Particle Number Density calculations # # keep using old microphys scheme here so far. as this is hard-wired inside # # jacobianDoit, i.e. we have to keep/set the scat_species_*_*_field(s) anyways # pnd_fieldCalcFromscat_speciesFields # WriteXML( in=pnd_field, filename="ScatSpec.pnd_field.xml" ) # #pnd_fieldZero # #Copy( dpnd_dummy, dpnd_field_dx ) #need this for getting dpnd_field in the correct size # #pnd_fieldCalcFromParticleBulkProps # #Copy( dpnd_field_dx, dpnd_dummy ) # # # Consistency checks # atmgeom_checkedCalc # atmfields_checkedCalc( bad_partition_functions_ok = 1 ) # cloudbox_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 # # # Initialize Doit variables # Copy( iy_main_agenda, iy_main_agenda__Emission ) # 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 # StringSet( info, "DoitCalc" ) # Print( info, 0 ) # DoitCalc # WriteXMLIndexed( in=doit_i_field, file_index=ybatch_index, filename="DOIT.doit_i_field" ) # WriteXML( in=scat_za_grid, filename="DOIT.scat_za_grid.xml" ) # # # # 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 ) # # StringSet( info, "jacobianDoit" ) # Print( info, 0 ) # jacobianDoit( ScatSpeciesMerge_do=1 ) # WriteXMLIndexed( in=jacobian_quantities, file_index=ybatch_index, filename="DOIT.jacobian_quantities" ) # WriteXMLIndexed( in=jacobian, file_index=ybatch_index, filename="DOIT.jacobian" ) # Compare( yref, y, 1. ) # # # jacobianOff # jacobianDoitAddSpecies( species="abs_species.H2O-PWR98", unit="abs", dx=1e-6 ) # jacobianDoitAddSpecies( species="scat_species.RR.mass_density", unit="abs", dx=1e-7 ) # jacobianDoitClose # Copy( scat_data, scat_data_orig ) # Copy( scat_meta, scat_meta_orig ) # Copy( scat_species, scat_species_orig ) # Copy( pnd_field, pnd_field_orig ) # ScatSpeciesMerge # cloudbox_checkedCalc # DoitInit # DoitGetIncoming(rigorous=0) # doit_i_fieldSetClearsky # Copy( doit_conv_test_agenda, doit_conv_test_agenda_refcase ) # StringSet( info, "DoitCalc2" ) # Print( info, 0 ) # DoitCalc # WriteXMLIndexed( in=doit_i_field, file_index=ybatch_index, filename="DOIT2.doit_i_field" ) # WriteXML( in=scat_za_grid, filename="DOIT2.scat_za_grid.xml" ) # Copy( doit_conv_test_agenda, doit_conv_test_agenda_perturb ) # Copy( scat_data, scat_data_orig ) # Copy( scat_meta, scat_meta_orig ) # Copy( scat_species, scat_species_orig ) # Copy( pnd_field, pnd_field_orig ) # StringSet( info, "jacobianDoit2" ) # Print( info, 0 ) # jacobianDoit( ScatSpeciesMerge_do=1 ) # WriteXMLIndexed( in=jacobian_quantities, file_index=ybatch_index, filename="DOIT2.jacobian_quantities" ) # WriteXMLIndexed( in=jacobian, file_index=ybatch_index, filename="DOIT2.jacobian" ) } # 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( "ascii", ybatch, "ybatchREFERENCE_iyHybrid.xml" ) #WriteXML( "ascii", ybatch_jacobians, "ybatch_jacobiansREFERENCE_iyHybrid.xml" ) # 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" ) ReadXML( ybatch_ref, "ybatchREFERENCE_iyHybrid.xml" ) CompareRelative( ybatch, ybatch_ref, 1e-9 ) ArrayOfMatrixCreate( ybatch_jacobians_ref ) #ReadXML( ybatch_jacobians_ref, "ybatch_jacobiansREFERENCE_DoitJacobians.xml" ) ReadXML( ybatch_jacobians_ref, "ybatch_jacobiansREFERENCE_iyHybrid.xml" ) CompareRelative( ybatch_jacobians, ybatch_jacobians_ref, 1e-6 ) }