# # Testing functionality (meeting format requirements, etc.) of surface related # data. # # General test setup: reading in raw data (including a basic atmosphere), # extracting/regridding, executing standard pre-RT calc internal test method # atmfields_checkedCalc, and performing some RT simulations in order to apply data # that has no dedicated check method (e.g., when only used through agendas like # surface reflectivity). # # # This case is for Earth and specifically tests # # - surface altitude: reading in raw 3D data, preprocess into z_surface variable # - within a 3D (global) case (CASE A) # - as 1D (extract one point) case (CASE B) # - land-sea-mask: used as basic data to derive land/sea dependent # reflectivities (maps mask value to refl=0.05 for land and 0.4 for sea) # - reading 3D raw data # - applying surface_rt_prop agenda on-the run extraction of mask value and # derivation of reflectivity at surface reflection point (CASE A-2; # CASE A-1 using fixed surface reflectivity for comparison) # - extracting at single point & used in 1D case (CASE B; sensor viewing # and lat/lon_true adjusted such that similar region observed as for # CASEs A) # # - CASEs A-1, A-2, and B RT calculation results compared for consistency # (setup such that A-1 and A-2 have the exactly same sensor setup, A-2 # looking to ocean and A-1 reflectivity set to the one of ocean; B has same # viewing geometry looking at similar point (but one point only!) as # cases A. Hence, results of A-1 and A-2 shall be identical, while results # from B shall be "close".) # # Jana Mendrok 2013-02-26 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) # parameters for converting land/sea mask to surface reflectivities # mask has land=1, sea=0 # later on we will convert those values to refl_land and refl_sea by # refl=a*maskval+b, where a=refl_land-refl_sea and b=refl_sea) NumericCreate( refl_land ) NumericCreate( refl_sea ) NumericSet( refl_land, 0.05 ) #emiss=0.95 NumericSet( refl_sea, 0.4 ) #emiss=0.60 IndexSet( lambertian_nza, 4 ) # the 3D geo grid to test VectorCreate( lat_grid3D ) VectorCreate( lon_grid3D ) VectorLinSpace( lat_grid3D, -90, 90, 18 ) VectorLinSpace( lon_grid3D, 0, 360, 18 ) #VectorLinSpace( lon_grid3D, 20, 90, 2 ) # alternatively, the geocoordinates for the 1D case (before switching on, check # to not overwrite them further down) #VectorSet( lat_true, [32.] ) #VectorSet( lon_true, [40.] ) GriddedField2Create( gf2tmp ) MatrixCreate( mtmp ) VectorCreate( vtmp ) NumericCreate( ntmp ) IndexCreate( itmp ) StringCreate( atmcase ) StringSet( atmcase, "planets/Earth/Fascod/tropical/tropical" ) #StringSet( atmcase, "planets/Earth/Fascod/midlatitude-summer/midlatitude-summer" ) StringCreate( surfpath ) StringSet( surfpath, "planets/Earth/ERA40/" ) StringCreate( zsurfname ) StringSet( zsurfname, "SurfaceAltitude_ERA40_1.0Degree" ) StringCreate( surfmaskname ) StringSet( surfmaskname, "LandSeaMask_ERA40_1.0Degree" ) #StringSet( surfmaskname, "LandSeaMask_ERA40_0.25Degree" ) StringCreate( casefull ) StringCreate( caseext ) # some stuff to get a basic atmosphere ##### # 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, atmcase ) p_gridFromZRaw( p_grid, z_field_raw, 0 ) # and some further settings in order to be able to do an RT calc ##### cloudboxOff jacobianOff IndexSet( stokes_dim, 1 ) VectorSet( f_grid, [300e9] ) sensorOff StringSet( iy_unit, "PlanckBT" ) abs_linesReadFromSplitArtscat(abs_lines, abs_species, "spectroscopy/Perrin/", 0, 1e12) abs_lines_per_speciesCreateFromLines # and agenda settings needed for RT calc ##### Copy( iy_main_agenda, iy_main_agenda__Emission ) Copy( ppath_agenda, ppath_agenda__FollowSensorLosPath ) Copy( blackbody_radiation_agenda, blackbody_radiation_agenda__Planck ) Copy( propmat_clearsky_agenda, propmat_clearsky_agenda__OnTheFly ) Copy( iy_space_agenda, iy_space_agenda__CosmicBackground ) Copy( iy_surface_agenda, iy_surface_agenda__UseSurfaceRtprop ) Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath ) # LOS zenith angle MatrixSet( sensor_los, [180;130;115;113.8] ) # LOS azimuth angle #MatrixSet( mtmp, []) nrowsGet( itmp, sensor_los ) MatrixSetConstant( mtmp, itmp, 1, 90. ) Append( sensor_los, mtmp, "trailing" ) # sensor altitude MatrixSetConstant( sensor_pos, itmp, 1, 600e3 ) # sensor latitude MatrixSetConstant( mtmp, itmp, 1, 0. ) Append( sensor_pos, mtmp, "trailing" ) # sensor longitutde MatrixSetConstant( mtmp, itmp, 1, 170. ) Append( sensor_pos, mtmp, "trailing" ) ##### # CASEs A: we start out with 3D (both fields are 3D) ##### AtmosphereSet3D Copy( lat_grid, lat_grid3D ) Copy( lon_grid, lon_grid3D ) AtmFieldsCalcExpand1D( vmr_zeropadding=1 ) # reading the surface altitude field Copy( casefull, surfpath ) Append( casefull, zsurfname ) ReadXML( gf2tmp, casefull ) GriddedFieldLatLonRegrid( gf2tmp, lat_grid, lon_grid, gf2tmp ) FieldFromGriddedField( z_surface, p_grid, lat_grid, lon_grid, gf2tmp ) atmfields_checkedCalc atmgeom_checkedCalc # reading surface mask; no regridding here. we need to do that as 1D inside the # surface_rtprop_agenda, as we only know then which exact point(s) we need. for # better traceability (and since here this isn't just a temporary field), we # create a dedicated workspace variable for the mask data. Copy( casefull, surfpath ) Append( casefull, surfmaskname ) GriddedField2Create( lsmask ) ReadXML( lsmask, casefull ) abs_xsec_agenda_checkedCalc propmat_clearsky_agenda_checkedCalc # now we need to do some RT calc in order to APPLY the reflectivity data cloudbox_checkedCalc sensor_checkedCalc # CASE A-1 ##### VectorSet( surface_scalar_reflectivity, [0.4] ) Copy( surface_rtprop_agenda, surface_rtprop_agenda__Specular_NoPol_ReflFix_SurfTFromt_field ) yCalc #WriteXML( "ascii", y, "TestSurf_Earth.y.3D.fixRef.xml" ) VectorCreate( y_3D_fixRef ) Copy( y_3D_fixRef, y ) #Print( y ) # CASE A-2 ##### AgendaSet( surface_rtprop_agenda ){ specular_losCalc InterpAtmFieldToPosition( out=surface_skin_t, field=t_field ) Select( lat_true, rtp_pos, [1] ) Select( lon_true, rtp_pos, [2] ) GriddedFieldLatLonRegrid( lsmask, lat_true, lon_true, lsmask ) FieldFromGriddedField( mtmp, p_grid, lat_true, lon_true, lsmask ) VectorExtractFromMatrix( surface_scalar_reflectivity, mtmp, 0, "row" ) Copy( ntmp, refl_sea ) NumericScale( ntmp, ntmp, -1. ) NumericAdd( ntmp, ntmp, refl_land ) VectorScale( surface_scalar_reflectivity, surface_scalar_reflectivity, ntmp ) VectorAddScalar( surface_scalar_reflectivity, surface_scalar_reflectivity, refl_sea ) #Print( rte_pos, 0 ) #Print( surface_scalar_reflectivity, 0 ) surfaceFlatScalarReflectivity } yCalc #WriteXML( "ascii", y, "TestSurf_Earth.y.3D.extRef.xml" ) VectorCreate( y_3D_extRef ) Copy( y_3D_extRef, y ) #Print( y ) Compare( y_3D_fixRef, y_3D_extRef, 1e-6 ) ##### # CASE B: now for the fun of it, we do a 1D case, too ##### # basic atmospheric fields AtmosphereSet1D AtmFieldsCalc( vmr_zeropadding=1 ) # sensor in 1D (zenith angles like for 3D. no azimtuh of course) VectorExtractFromMatrix( lat_true, sensor_pos, 1, "column") Select( lat_true, lat_true, [0] ) VectorExtractFromMatrix( lon_true, sensor_pos, 2, "column") Select( lon_true, lon_true, [0] ) VectorExtractFromMatrix( vtmp, sensor_los, 0, "column" ) Matrix1ColFromVector( sensor_los, vtmp ) nrowsGet( itmp, sensor_los ) MatrixSetConstant( sensor_pos, itmp, atmosphere_dim, 600e3 ) # reading the surface altitude field (need to do that again as we didn't keep # the "raw" gridded field) and extracting a single point for 1D Copy( casefull, surfpath ) Append( casefull, zsurfname ) ReadXML( gf2tmp, casefull ) GriddedFieldLatLonRegrid( gf2tmp, lat_true, lon_true, gf2tmp ) FieldFromGriddedField( z_surface, p_grid, lat_grid, lon_grid, gf2tmp ) #Print( z_surface ) atmfields_checkedCalc atmgeom_checkedCalc # reading surface mask (need to do that again as we didn't keep the "raw" # gridded field) and extracting a single point for 1D Copy( casefull, surfpath ) Append( casefull, surfmaskname ) ReadXML( gf2tmp, casefull ) GriddedFieldLatLonRegrid( gf2tmp, lat_true, lon_true, gf2tmp ) FieldFromGriddedField( mtmp, p_grid, lat_grid, lon_grid, gf2tmp ) VectorExtractFromMatrix( surface_scalar_reflectivity, mtmp, 0, "row" ) # converting the land/sea mask value to surface reflectivity Copy( ntmp, refl_sea ) NumericScale( ntmp, ntmp, -1. ) NumericAdd( ntmp, ntmp, refl_land ) VectorScale( surface_scalar_reflectivity, surface_scalar_reflectivity, ntmp ) VectorAddScalar( surface_scalar_reflectivity, surface_scalar_reflectivity, refl_sea ) #Print( surface_scalar_reflectivity ) # now we need to do some RT calc in order to APPLY the reflectivity data cloudbox_checkedCalc sensor_checkedCalc #Copy( surface_rtprop_agenda, surface_rtprop_agenda__Blackbody_SurfTFromt_field ) #Copy( surface_rtprop_agenda, surface_rtprop_agenda__lambertian_ReflFix_SurfTFromt_field ) Copy( surface_rtprop_agenda, surface_rtprop_agenda__Specular_NoPol_ReflFix_SurfTFromt_field ) yCalc #WriteXML( "ascii", y, "TestSurf_Earth.y.1D.xml" ) #Print( y ) Compare( y, y_3D_extRef, 0.1 ) }