#DEFINITIONS: -*-sh-*- # -------------------------------------------------------------------- # Test batch calculations for HIRS. # Viju Oommen John and Stefan Buehler, August to October 2008. # -------------------------------------------------------------------- Arts2 { INCLUDE "general.arts" INCLUDE "continua.arts" StringCreate(satellite) ArrayOfIndexCreate(channels) ArrayOfIndexCreate(views) StringCreate(hitran_file) NumericCreate(f_grid_spacing) # 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]) #ArrayOfIndexSet(channels, [11]) # Select here which views you want # --- ArrayOfIndexSet(views, [0, 7, 14, 21, 27]) # You have to give here the location of the HITRAN catalogue # --- # FIXME: This has to be replaced by a little piece of the catalog # explicitly included here, or? # We have to discuss with Oliver (and perhaps Patrick) how to handle # this. # We can also include a pre-calculated absorption table for HIRS with # ARTS. StringSet(hitran_file,"/storage3/data/catalogue/hitran/hitran2004/HITRAN04.par") #StringSet(hitran_file,"/home/patrick/Data/HITRAN_2004/HITRAN04.par") # Set frequency grid spacing # (See comments in hirs_reference.arts concerning useful values here) # --- NumericSet(f_grid_spacing, 5e8) INCLUDE "hirs/hirs_fast.arts" # Set up absorption # ================= # Atmospheric profiles # --- # Atmospheric profiles are stored in an ArrayOfMatrix. # It contains one matrix for each atmospheric state. # Each matrix row corresponds to one pressure level. The # meaning of the columns is: # p[Pa] T[K] z[m] H2O[VMR] O3[VMR] # ArrayOfMatrixCreate( arrayofmatrix_1 ) #ReadXML( arrayofmatrix_1, # "/storage2/home/sbuehler/checkouts/arts-xml-data/atmosphere/chevallier_91L/chevallierl91_clear_q.xml" ) #ReadXML( arrayofmatrix_1, "../AMSU/chevallierl91_clear_q_extract.xml" ) ReadXML( arrayofmatrix_1, "garand_profiles.xml" ) # Storage in an array of matrix is handy for profile data, because # it is very compact. However, the more general internal # representation of the data is in batch_atm_fields_compact. # Convert to batch_atm_fields_compact # --- # The values taken for O2 and N2 are from Wallace&Hobbs, 2nd edition. batch_atm_fields_compactFromArrayOfMatrix( batch_atm_fields_compact, atmosphere_dim, arrayofmatrix_1, ["T", "z", "H2O", "O3", "CO2", "N2O", "CO", "CH4"], ["O2", "N2"], [0.2095, 0.7808] ) # Delete original data array to conserve memory: # --- Delete( arrayofmatrix_1 ) # Set parameters for lookup table # --- # Arguments omitted for better maintainability of this test file. abs_lookupSetupBatch # abs_lookupSetupBatch( abs_p, # abs_t, # abs_t_pert, # abs_vmrs, # abs_nls, # abs_nls_pert, # abs_species, # batch_atm_fields_compact, # abs_t_interp_order, # abs_nls_interp_order, # , # , # , # ) # 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_lookupCreate #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 abs_scalar_gas_agenda to use lookup table # --- AgendaSet( abs_scalar_gas_agenda ){ abs_scalar_gasExtractFromLookup } # Set up RT calculation # ===================== # If you do not want to use a blackbody surface, uncomment this # --- #VectorSetConstant( surface_scalar_reflectivity, 1, 0.01 ) #AgendaSet( surface_prop_agenda ){ # Ignore( rte_pos ) # Ignore( rte_gp_p ) # InterpSurfaceFieldToRteGps( surface_skin_t, atmosphere_dim, # rte_gp_lat, rte_gp_lon, t_surface ) # surfaceFlatReflectivity( surface_los, surface_rmatrix, surface_emission, # f_grid, stokes_dim, atmosphere_dim, rte_los, # surface_skin_t, surface_scalar_reflectivity ) #} # 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 two lines if you want to set it to something else. #MatrixSetConstant( sensor_pos, 1, 1, 850e3 ) #sensor_posAddRgeoid # 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: # # 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 # get some surface properties from corresponding atmospheric fields Extract( z_surface, z_field, 0 ) Extract( t_surface, t_field, 0 ) # Perform RT calculations # --- basics_checkedCalc cloudbox_checkedCalc yCalc # Optionally write out jacobian: # WriteXMLIndexed( output_file_format, ybatch_index, jacobian, "" ) # Convert the measurement from radiance units to Planck Tb: StringSet( y_unit, "PlanckBT" ) y_unitApply } # Set number of batch cases: nelemGet( ybatch_n, batch_atm_fields_compact ) #IndexSet(ybatch_n, 1) # Execute the batch calculations: # --- ybatchCalc # Store result matrix: # --- WriteXML( "ascii", ybatch ) }