#DEFINITIONS: -*-sh-*- # -------------------------------------------------------------------- # Test batch calculations for HIRS. # Viju Oommen John and Stefan Buehler, August to October 2008. # -------------------------------------------------------------------- Arts2 { INCLUDE "general/general.arts" INCLUDE "general/continua.arts" INCLUDE "general/agendas.arts" INCLUDE "general/planet_earth.arts" # 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 ) # blackbody surface with skin temperature interpolated from t_surface field Copy( surface_rtprop_agenda, surface_rtprop_agenda__Blackbody_SurfTFromt_field ) jacobianOff StringCreate( satellite ) ArrayOfIndexCreate( channels ) ArrayOfIndexCreate( views ) # Select here which satellite you want # --- StringSet( satellite, "NOAA14" ) # Select here which channels you want # --- # # HIRS Channels 13-19 are shortwave channels. Simulating them with # ARTS for thermal radiation only is pointless. # # WATCH OUT! We use the usual zero-based ARTS indexing here, so index # 11 is actually channel 12, etc. #ArrayOfIndexSet( channels, [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18] ) ArrayOfIndexSet( channels, [0,1,2,3,4,5,6,7,8,9,10,11] ) # Select here which views you want # --- ArrayOfIndexSet( views, [0, 7, 14, 21, 27] ) # Basic settings # --- # This example assumes 1D AtmosphereSet1D # Scalar RT IndexSet( stokes_dim, 1 ) # HIRS spectral response functions assume radiance. StringSet( iy_unit, "1" ) # Set up absorption # --- abs_speciesSet( species=[ "H2O, H2O-SelfContCKDMT100, H2O-ForeignContCKDMT100", "O3", "CO2, CO2-CKDMT100", "N2O", "CO", "CH4", "O2, O2-CIAfunCKDMT100", "N2, N2-CIAfunCKDMT100, N2-CIArotCKDMT100" ] ) abs_lineshapeDefine( abs_lineshape, "Voigt_Kuntz6", "VVH", 750e9 ) # deriving abs_lines separate, because one might want to use a lookup table # # no use of HITRAN in test cases # (uncomment when you want to use HITRAN) #abs_linesReadFromHitran( abs_lines, # "HITRAN2012.par", # 2.0022234e+12, 7.9855733e+14 ) # instead we use a HITRAN-extract converted to ARTSCAT # (comment out when you want to use HITRAN) ReadXML( abs_lines, "testdata/abs_lines_IR.xml.gz" ) abs_lines_per_speciesCreateFromLines # Sensor setup # --- # Definition of sensor position and LOS ReadXML( sensor_los, "hirs.sensor_los.xml" ) # Select those views that are requested by the user Select( sensor_los, sensor_los, views ) nrowsGet( nrows, sensor_los ) ncolsGet( ncols, sensor_los ) MatrixSetConstant( sensor_pos, nrows, ncols, 850e3 ) # Start sensor response setup # Normalise the sensor response IndexSet( sensor_norm, 1 ) AntennaOff # Construct names of sensor description files # Nominal channel frequencies StringCreate( f_backend_file) StringJoin( f_backend_file, satellite, "_HIRS.f_backend.xml" ) # Fast frequency grid StringCreate( f_grid_file) StringJoin( f_grid_file, satellite, "_HIRS.f_grid_fast.xml" ) # Weights associated with each frequency StringCreate( weights_file) StringJoin( weights_file, satellite, "_HIRS.W_fast.xml" ) # Spectrometer ReadXML( f_backend, f_backend_file ) # Read in optimized frequency grid ReadXML( f_grid, f_grid_file ) # Read in and WMRF weights: ReadXML( wmrf_weights, weights_file ) # Select only the active channels: Copy( wmrf_channels, channels ) # The method acts on several variables: # f_grid, wmrf_weights, and f_backend WMRFSelectChannels # Initialize sensor variables. sensor_responseInit # Add WMRF weigths to sensor response: sensor_responseWMRF # End of sensor response setup # Compact line list, to kick out lines that are outside there # cutoff range for all frequencies. abs_lines_per_speciesCompact # Atmospheric profiles # --- # Atmospheric profiles are stored in an ArrayOfGriddedField4. # It contains one GriddedField4 for each atmospheric state. # ReadXML( batch_atm_fields_compact, "testdata/garand_profiles.xml.gz" ) # 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 ) # Set up lookup table # --- # Set parameters for lookup table # Arguments omitted for better maintainability of this test file. abs_lookupSetupBatch # Optional manual override of T and VMR perturbations # # If your input data contains extreme outliers, the lookup table will # get unreasonably large. It is suggested that you instead set them # manually here. The Chevallier 91L data (humidity set) contains # temperature variations from -70 to +55 (K) and humidity variations from # 0 to 6 (fractional units). This should be the reasonable range of # atmospheric variability. You will get failures from individual jobs in # the batch, that are outside this range. #VectorLinSpace( abs_t_pert, -70, 55, 5 ) #VectorLinSpace( abs_nls_pert, 0, 6, 0.5 ) # Create the lookup table abs_xsec_agenda_checkedCalc abs_lookupCalc #WriteXML( "binary", abs_lookup ) # Test (and print) lookup table accuracy # # This method is not necessary, it just tests and prints the lookup # table accuracy. Comment it out if you do not want this # information. The test should take approximately as much time as the # lookup table generation, perhaps even more. So, it is not cheap! #abs_lookupTestAccuracy # Set propmat_clearsky_agenda to use lookup table # --- Copy( propmat_clearsky_agenda, propmat_clearsky_agenda__LookUpTable ) # Set up RT calculation # ===================== # Definition of sensor position and LOS # --- # Optionally set sensor_pos # # The sensor altitude is predefined in amsu.arts to 850 km above the geoid. # Comment out the next line if you want to set it to something else. #MatrixSetConstant( sensor_pos, 1, 1, 850e3 ) # Set the agenda for batch calculations: # --- AgendaSet( ybatch_calc_agenda ){ # Extract the atmospheric profiles for this case: Extract( atm_fields_compact, batch_atm_fields_compact, ybatch_index ) # Split up *atm_fields_compact* to # generate p_grid, t_field, z_field, vmr_field: AtmFieldsFromCompact # Optionally set Jacobian parameters. # uncomment this for NO jacobian calculations jacobianOff # Uncomment this block if you want Jacobians. Attention, it slows down the # computation a lot. # Also, you can add other Jacobians here, for example for temperature. # # jacobianInit # jacobianAddAbsSpecies( jacobian_quantities, jacobian_agenda, # atmosphere_dim, # p_grid, lat_grid, lon_grid, # p_grid, lat_grid, lon_grid, # "H2O, H2O-SelfContCKDMT100, H2O-ForeignContCKDMT100", # "analytical", # "rel", # 0.01 ) # jacobianClose # No scattering cloudboxOff # get some surface properties from corresponding atmospheric fields Extract( z_surface, z_field, 0 ) Extract( t_surface, t_field, 0 ) # Perform RT calculations # --- atmfields_checkedCalc( bad_partition_functions_ok = 1 ) atmgeom_checkedCalc cloudbox_checkedCalc sensor_checkedCalc yCalc # Optionally write out jacobian: # WriteXMLIndexed( output_file_format, ybatch_index, jacobian, "" ) # Convert the measurement from radiance units to Planck Tb: StringSet( iy_unit, "PlanckBT" ) yApplyUnit } # Set number of batch cases: nelemGet( ybatch_n, batch_atm_fields_compact ) #IndexSet( ybatch_n, 1 ) # Execute the batch calculations: # --- propmat_clearsky_agenda_checkedCalc ybatchCalc # Store result matrix: # --- WriteXML( "ascii", ybatch ) # Compare with reference: ArrayOfVectorCreate( ybatch_ref ) ReadXML( ybatch_ref, "TestHIRS.NOAA14.ybatch.ref.xml" ) Compare( ybatch_ref, ybatch, 0.01 ) }