# # Testing functionality (meeting format requirements, etc.) of supplemental # atmospheric scenario data. # # General test setup: reading in raw data (including a basic atmosphere), # expanding (where necessary), regridding (incl. extracting), executing # standard pre-RT calc internal test method atmfields_checkedCalc. # # # This case is for Earth and specifically tests # # - electron densities (given: 4seasons x 8dayttimes x 2solaractivity global 3D # cases) # - all 64 cases in global 3D with p_grid from altitude grid taken from # standard FASCOD, i.e. extending up to 95km only (CASEs A) and p_grid # from altitude grid taken from higher altitude extensions (up to 2000km) # to FASCOD (CASEs B) # - magnetic field (given: 4decadal 3D cases, each with with 3 components) # - global 3D case with p_grid from standard and expanded FASCOD altitude # grids (CASEs A and B) # # NOTE: beside applying them here as "background" atmosphere, the FASCOD # atmospheric field data is not tested. This dataset has been part of ARTS # xml-data for long time, applied in calculations, hence further testing is not # necessary. # # Jana Mendrok 2013-02-26 Arts2 { INCLUDE "general/general.arts" INCLUDE "general/agendas.arts" INCLUDE "general/planet_earth.arts" Tensor3Create( edensity_field ) GriddedField3Create( gf3tmp ) StringCreate( caseext ) StringCreate( casefull ) StringCreate( basename ) StringCreate( atmcase ) IndexCreate( ncases ) IndexCreate( interp_order ) IndexSet( interp_order, 1 ) # Array with base case names ArrayOfStringCreate( basecasearray ) ArrayOfStringSet( basecasearray, ["planets/Earth/Fascod/", "planets/Earth/IRI/IRI", "planets/Earth/IGRF/IGRF11"] ) # "/storage4/home/mendrok/projects/MicrowavePropagationToolbox/WorkPackageData/WP2000/Earth/MagField/ToolBoxdata/IGRF11"] ) #StringSet( atmcase, "midlatitude-summer" ) StringSet( atmcase, "tropical" ) # Arrays with Ne (sub)case names ArrayOfStringCreate( seasoncasearray ) ArrayOfStringSet( seasoncasearray, ["_spring", "_summer", "_fall", "_winter"] ) ArrayOfStringCreate( timecasearray ) ArrayOfStringSet( timecasearray, ["_00UTC.Ne", "_03UTC.Ne", "_06UTC.Ne", "_09UTC.Ne", "_12UTC.Ne", "_15UTC.Ne", "_18UTC.Ne", "_21UTC.Ne"] ) ArrayOfStringCreate( solarcasearray ) ArrayOfStringSet( solarcasearray, ["_solmax", "_solmin"] ) # Arrays with B (sub)case names ArrayOfStringCreate( yearcasearray ) ArrayOfStringSet( yearcasearray, ["_1980", "_1990", "_2000", "_2010"] ) #, "_2012"] ) ArrayOfStringCreate( resolcasearray ) ArrayOfStringSet( resolcasearray, ["_200km-5deg-5deg"] ) # the 3D geo grid to test VectorCreate( lat_grid3D ) VectorCreate( lon_grid3D ) VectorLinSpace( lat_grid3D, -90, 90, 18 ) VectorLinSpace( lon_grid3D, -20, 350, 18 ) VectorLinSpace( lon_grid3D, 0, 360, 18 ) ##### # first electron densities, 3D ##### ### #this is the real stuff, partI: electron density # We have to repeat that, as we have multiple Ne fields given, but each atm case # can have only one (in contrast to vmr profiles, where we could test all # profiles in a single run). We use a bunch of forloops for this, which we # define here and execute later on. AgendaCreate( forloop_agenda_time ) AgendaSet( forloop_agenda_time ){ # construct case name Extract( casefull, timecasearray, forloop_index ) Append( atmcase, casefull ) Print( atmcase, 0 ) # readin in raw data ReadXML( gf3tmp, atmcase ) #this is 3D data. we need to regrid the raw data to the calculation grid(s). # for supplemental atm data (Ne, B, winds) this requires manual regridding (in # contrast to basic atm data, where this is handled by AtmFieldsCalc(Expand1D). GriddedFieldPRegrid( gf3tmp, p_grid, gf3tmp, interp_order, 1 ) GriddedFieldLatLonRegrid( gf3tmp, lat_grid, lon_grid, gf3tmp, interp_order ) FieldFromGriddedField( edensity_field, p_grid, lat_grid, lon_grid, gf3tmp ) atmfields_checkedCalc atmgeom_checkedCalc #WriteXML( "ascii", p_grid ) #WriteXML( "ascii", z_field ) #WriteXML( "ascii", t_field ) #WriteXMLIndexed( "ascii", forloop_index, edensity_field ) } AgendaCreate( forloop_agenda_season ) AgendaSet( forloop_agenda_season ){ # construct case name Extract( casefull, seasoncasearray, forloop_index ) Append( atmcase, casefull ) #Print( atmcase, 0 ) Copy( forloop_agenda, forloop_agenda_time ) nelemGet( ncases, timecasearray ) IndexStepDown( ncases, ncases ) ForLoop( forloop_agenda, 0, ncases, 1 ) } AgendaCreate( forloop_agenda_solar ) AgendaSet( forloop_agenda_solar ){ # construct case name Copy( atmcase, basename ) Extract( casefull, solarcasearray, forloop_index ) Append( atmcase, casefull ) #Print( atmcase, 0 ) Copy( forloop_agenda, forloop_agenda_season ) nelemGet( ncases, seasoncasearray ) IndexStepDown( ncases, ncases ) ForLoop( forloop_agenda, 0, ncases, 1 ) } # now we actually execute things... # 3-dimensional atmosphere AtmosphereSet3D Copy( lat_grid, lat_grid3D ) Copy( lon_grid, lon_grid3D ) # set atmospheric scenario Extract( basename, basecasearray, 0 ) Append( basename, atmcase ) StringSet( caseext, "/" ) Append( basename, caseext ) Append( basename, atmcase ) # we manually select a minumim set of basic atm data (main atm constituents) abs_speciesSet( species=["O2", "N2"] ) AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw, abs_species, basename ) p_gridFromZRaw( p_grid, z_field_raw, 0 ) AtmFieldsCalcExpand1D( vmr_zeropadding=1 ) Extract( z_surface, z_field, 0 ) # now doing the Ne cases (reading, 1D regridding, checking). we use several # forloops to loop over the different aspects of the data Extract( basename, basecasearray, 1 ) Copy( forloop_agenda, forloop_agenda_solar ) nelemGet( ncases, solarcasearray ) IndexStepDown( ncases, ncases ) ForLoop( forloop_agenda, 0, ncases, 1 ) # now we want to use a larger p_grid (for using the full vertical extend of the # Ne data). that means, we need to get larger extending z_field and t_field. we # can derive a z_field from p_grid and t_field, but we need a larger t_field. as # manipulating fields within ARTS controlfiles is a hassle, we made extended raw # t/z fields that as far as possible agree with the assumptions made when # creating the Ne and B fields (IRI and IGRF data). we replace the standard t/z # raw field by these. # going back to clearsky atmosphere and derive expanded t/z Extract( basename, basecasearray, 0 ) Append( basename, atmcase ) StringSet( caseext, "/" ) Append( basename, caseext ) Append( basename, atmcase ) Copy( casefull, basename ) StringSet( caseext, ".expanded.t" ) Append( casefull, caseext ) Print( casefull, 0 ) ReadXML( t_field_raw, casefull ) Copy( casefull, basename ) StringSet( caseext, ".expanded.z" ) Append( casefull, caseext ) ReadXML( z_field_raw, casefull ) p_gridFromZRaw( p_grid, z_field_raw, 0 ) # make a higher resolution p_grid spread over input z_field (i.e.,with same # upper and lower boundaries, but more grid points) IndexCreate( itmp ) NumericCreate( ntmp1 ) NumericCreate( ntmp2 ) Extract( ntmp1, p_grid, 0 ) nelemGet( itmp, p_grid ) IndexStepDown( itmp, itmp ) Extract( ntmp2, p_grid, itmp ) VectorNLogSpace( p_grid, 201, ntmp1, ntmp2) #WriteXML( "ascii", p_grid ) AtmFieldsCalcExpand1D( vmr_zeropadding=1 ) Extract( z_surface, z_field, 0 ) # now doing the Ne cases (reading, 1D regridding, checking). we use several # forloops to loop over the different aspects of the data Extract( basename, basecasearray, 1 ) Copy( forloop_agenda, forloop_agenda_solar ) nelemGet( ncases, solarcasearray ) IndexStepDown( ncases, ncases ) ForLoop( forloop_agenda, 0, ncases, 1 ) ##### # second: magnetic field, 3D ##### AgendaCreate( forloop_agenda_Bres ) AgendaSet( forloop_agenda_Bres ){ # construct case name Extract( casefull, resolcasearray, forloop_index ) Append( atmcase, casefull ) Print( atmcase, 0 ) #this is 3D data. we need to regrid the raw data to the calculation grid(s). # for supplemental atm data (Ne, B, winds) this requires manual reading in and # regridding (in contrast to basic atm data, where this is handled by # AtmFieldsCalc(Expand1D). #now reading and regridding of each of the 3 B-field components separately Copy( casefull, atmcase ) StringSet( caseext, ".B_u" ) Append( casefull, caseext ) ReadXML( gf3tmp, casefull ) GriddedFieldPRegrid( gf3tmp, p_grid, gf3tmp, interp_order, 1 ) GriddedFieldLatLonRegrid( gf3tmp, lat_grid, lon_grid, gf3tmp, interp_order ) FieldFromGriddedField( mag_u_field, p_grid, lat_grid, lon_grid, gf3tmp ) #now reading and regridding of each of the 3 B-field components separately Copy( casefull, atmcase ) StringSet( caseext, ".B_v" ) Append( casefull, caseext ) ReadXML( gf3tmp, casefull ) GriddedFieldPRegrid( gf3tmp, p_grid, gf3tmp, interp_order, 1 ) GriddedFieldLatLonRegrid( gf3tmp, lat_grid, lon_grid, gf3tmp, interp_order ) FieldFromGriddedField( mag_v_field, p_grid, lat_grid, lon_grid, gf3tmp ) #now reading and regridding of each of the 3 B-field components separately Copy( casefull, atmcase ) StringSet( caseext, ".B_w" ) Append( casefull, caseext ) ReadXML( gf3tmp, casefull ) GriddedFieldPRegrid( gf3tmp, p_grid, gf3tmp, interp_order, 1 ) GriddedFieldLatLonRegrid( gf3tmp, lat_grid, lon_grid, gf3tmp, interp_order ) FieldFromGriddedField( mag_w_field, p_grid, lat_grid, lon_grid, gf3tmp ) atmfields_checkedCalc atmgeom_checkedCalc #WriteXML( "ascii", p_grid ) #WriteXML( "ascii", z_field ) #WriteXML( "ascii", t_field ) #WriteXMLIndexed( "ascii", forloop_index, mag_u_field ) #WriteXMLIndexed( "ascii", forloop_index, mag_v_field ) #WriteXMLIndexed( "ascii", forloop_index, mag_w_field ) } AgendaCreate( forloop_agenda_Byears ) AgendaSet( forloop_agenda_Byears ){ # construct case name Copy( atmcase, basename ) Extract( casefull, yearcasearray, forloop_index ) Append( atmcase, casefull ) Print( atmcase, 0 ) Copy( forloop_agenda, forloop_agenda_Bres ) nelemGet( ncases, resolcasearray ) IndexStepDown( ncases, ncases ) ForLoop( forloop_agenda, 0, ncases, 1 ) } # now we actually execute things... # 3-dimensional atmosphere AtmosphereSet3D Copy( lat_grid, lat_grid3D ) Copy( lon_grid, lon_grid3D ) # set atmospheric scenario Extract( basename, basecasearray, 0 ) Append( basename, atmcase ) StringSet( caseext, "/" ) Append( basename, caseext ) Append( basename, atmcase ) # we manually select a minumim set of basic atm data (main atm constituents) Delete( vmr_field_raw ) abs_speciesSet( species=["O2", "N2"] ) AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw, abs_species, basename ) Copy( casefull, basename ) StringSet( caseext, ".expanded.t" ) Append( casefull, caseext ) ReadXML( t_field_raw, casefull ) Copy( casefull, basename ) StringSet( caseext, ".expanded.z" ) Append( casefull, caseext ) ReadXML( z_field_raw, casefull ) p_gridFromZRaw( p_grid, z_field_raw, 0 ) # make a higher resolution p_grid spread over input z_field (i.e.,with same # upper and lower boundaries, but more grid points) Extract( ntmp1, p_grid, 0 ) nelemGet( itmp, p_grid ) IndexStepDown( itmp, itmp ) Extract( ntmp2, p_grid, itmp ) VectorNLogSpace( p_grid, 201, ntmp1, ntmp2) #WriteXML( "ascii", p_grid ) AtmFieldsCalcExpand1D( vmr_zeropadding=1 ) Extract( z_surface, z_field, 0 ) atmfields_checkedCalc atmgeom_checkedCalc # now doing the B cases (reading, 1D regridding, checking). we use several # forloops to loop over the different aspects of the data Extract( basename, basecasearray, 2 ) Copy( forloop_agenda, forloop_agenda_Byears ) nelemGet( ncases, yearcasearray ) IndexStepDown( ncases, ncases ) ForLoop( forloop_agenda, 0, ncases, 1 ) }