2025-09-05 09:42:21 +02:00

222 lines
7.6 KiB
Prolog
Executable File

;+---------------------------------------------------------------------------
PRO read_cdfobs, Files, NumObs=NumObs, $
Latitudes=Latitudes, Longitudes=Longitudes, Depths=Depths, $
Obs=Observations, ModelVals=ModelVals, qcs=QCs, $
Ob2=Observations2, ModelVal2=ModelVals2, qc2=QCs2, $
Ob3=Observations3, $
Dates=Dates, rmdi=rmdi, iobs=iobs, jobs=jobs, kobs=kobs, $
Salinity=Salinity, nodates=nodates, types=types, $
filetype=ObsType, quiet=quiet, MDT=MDT, $
ProfileNum=ProfileNum, error=error, $
notfussy=notfussy,VarName=VarName
;+--------------------------------------------------------------------------
; Read in observation and feedback files
; detects filetype and calls the appropriate reading routine
;
; Author: D. J. Lea Feb 2008
;+--------------------------------------------------------------------------
; Declare error label.
ON_IOERROR, IOERROR
error=0
title=''
; Set types to undefined
types=-1 & tempvar=size(temporary(types))
; Read in netcdf data file and observation operator
;2. Work out which type of data is in the files by looking at the first one.
ncid = ncdf_open(Files(0), /nowrite)
if ncid lt 0 then message, 'Error opening file '+File
result = NCDF_ATTINQ( ncid, 'title', /global)
if result.datatype ne 'UNKNOWN' then ncdf_attget, ncid, 'title', Title, /global
if string(Title) eq "NEMO observation operator output" then ObsType='feedback' $
else $
if string(Title) eq "Forecast class 4 file" then ObsType='ForecastClass4' $
else $
ObsType = 'none'
if ObsType eq 'none' then begin
varid = ncdf_varid(ncid, 'POTM_CORRECTED')
if varid ge 0 then ObsType = 'Prof'
varid = ncdf_varid(ncid, 'SLA')
if varid ge 0 then ObsType = 'SSH'
varid = ncdf_varid(ncid, 'SST')
varid2 = ncdf_varid(ncid, 'sea_surface_temperature')
varid=max([varid,varid2])
if varid ge 0 then ObsType = 'SST'
varid = ncdf_varid(ncid, 'SEAICE')
varid2 = ncdf_varid(ncid, 'sea_ice_concentration')
varid=max([varid,varid2])
if varid ge 0 then ObsType = 'SEAICE'
varid = ncdf_varid(ncid, 'LOGCHL')
varid2 = ncdf_varid(ncid, 'LogChl')
varid3 = ncdf_varid(ncid, 'CHL1_mean')
varid=max([varid,varid2,varid3])
if varid ge 0 then ObsType = 'LOGCHL'
varid = ncdf_varid(ncid, 'PCO2')
if varid ge 0 then ObsType = 'PCO2'
varid = ncdf_varid(ncid, 'UCRT')
varid2 = ncdf_varid(ncid, 'VCRT')
varid=max([varid,varid2])
if varid ge 0 then ObsType = 'CRT'
endif
ncdf_close, ncid
if not keyword_set(quiet) then print, 'Reading in data as type ',ObsType
if ObsType eq 'feedback' then begin
if n_elements(DepRange) eq 0 then DepRange=[0,5000]
read_feedback, Files, DepRange=DepRange, VarName=VarName, NumData=numobs, $
OutLats=Latitudes, OutLons=Longitudes, $
OutInstruments=Instruments, OutPlatform=Platform, $
OutDeps=Depths, $
OutObs1=Observations, OutObs2=Observations2, $
OutMod1=ModelVals, OutMod2=ModelVals2, $
OutQC1=QCs, OutQC2=QCs2, $
MDT=MDT, OutDates=Dates, $
rmdi=rmdi, quiet=quiet, $
OutProfileNum=ProfileNum
; info, instruments
; info, platform
types=Platform
if n_elements(instruments) gt 0 then types=Platform+' '+Instruments
print, 'Varname: ',VarName
endif else if ObsType eq 'Prof' then begin
if keyword_set(PlotModel) then begin
OutTMod = 1
OutSMod = 1
endif
if n_elements(DepRange) eq 0 then DepRange=[0,5000]
read_enact, Files, DepRange=DepRange, NumData=NumData, $
OutLats=Latitudes, OutLons=Longitudes, $
Instruments=Instruments, Platform=Platform, $
OutDeps=Depths, OutSObs=OutSObs, OutSQC=OutSQC, $
OutTObs=OutTObs, OutTQC=OutTQC, OutDates=Dates, $
OutTMod=OutTMod, OutSMod=OutSMod, rmdi=rmdi, quiet=quiet, $
iobs=iobs, jobs=jobs, kobs=kobs, $
ProfileNum=ProfileNum
types=Instruments+' '+Platform
; if salinity keyword is set then put salinity values first
if keyword_set(Salinity) then begin
Observations = OutSObs
ModelVals = OutSMod
QCs = OutSQC
Observations2 = OutTObs
ModelVals2 = OutTMod
QCs2 = OutTQC
endif else begin
Observations = OutTObs
ModelVals = OutTMod
QCs = OutTQC
Observations2 = OutSObs
ModelVals2 = OutSMod
QCs2 = OutSQC
endelse
numobs=numdata
numobs=n_elements(latitudes)
endif else if ObsType eq 'SSH' then begin
read_sla, Files, NumObs=NumObs, $
Latitudes=Latitudes, Longitudes=Longitudes, $
ObsSLA=Observations, ModSLA=ModelVals, SLAQC=QCs, $
MDT=MDT, Satellites=Satellites, types=types, Dates=Dates, rmdi=rmdi, $
iobs=iobs, jobs=jobs, mstp=mstp, quiet=quiet
Depths = fltarr(NumObs)
endif else if ObsType eq 'SST' then begin
read_sst, Files, NumObs=NumObs, $
Latitudes=Latitudes, Longitudes=Longitudes, $
ObsSST=Observations, ModSST=ModelVals, SSTQC=QCs, $
Dates=Dates, rmdi=rmdi, types=types, iobs=iobs, jobs=jobs, $
quiet=quiet
Depths = fltarr(NumObs)
endif else if ObsType eq 'SEAICE' then begin
read_seaice, Files, NumObs=NumObs, $
Latitudes=Latitudes, Longitudes=Longitudes, $
Obs=Observations, Modarr=ModelVals, QC=QCs, $
Dates=Dates, rmdi=rmdi, iobs=iobs, jobs=jobs, nodates=nodates, $
quiet=quiet
Depths = fltarr(NumObs)
endif else if ObsType eq 'LOGCHL' then begin
read_chl, Files, NumObs=NumObs, $
Latitudes=Latitudes, Longitudes=Longitudes, $
Obs=Observations, Modarr=ModelVals, QC=QCs, $
Dates=Dates, rmdi=rmdi, iobs=iobs, jobs=jobs, nodates=nodates, $
quiet=quiet
Depths = fltarr(NumObs)
endif else if ObsType eq 'PCO2' then begin
read_pco2, Files, NumObs=NumObs, $
Latitudes=Latitudes, Longitudes=Longitudes, $
Obs=Observations, Modarr=ModelVals, QC=QCs, $
Dates=Dates, rmdi=rmdi, iobs=iobs, jobs=jobs, nodates=nodates, $
quiet=quiet
Depths = fltarr(NumObs)
endif else if ObsType eq 'CRT' then begin
read_crt, Files, NumObs=NumObs, $
OutLats=Latitudes, OutLons=Longitudes, $
OutTypes=Types, OutDates=Dates, $
OutUObs=Observations, OutVObs=Observations2, $
OutUMod=ModelVals, OutVMod=ModelVals2, Quiet=Quiet, $
FloatNum=FloatNum, QC=QCs, OutSpeed=Observations3
rmdi=0.0
QCs2=QCs
Depths = fltarr(NumObs)
endif else if ObsType eq 'ForecastClass4' then begin
print, 'reading forecast class 4 files'
read_forc_class4, Files, NumObs=NumObs, $
Latitudes=Latitudes, Longitudes=Longitudes, $
Depths=Depths, $
Types=types, Dates=Dates, $
Obsarr=Observations, Obs2arr=Observations2, $
Modarr=ModelVals, Mod2arr=ModelVals2, $
QCs1=QCs, QCs2=QCs2, rmdi=rmdi, $
notfussy=notfussy
endif else message, 'Error: ObsType is not set correctly'
if (n_elements(types) eq 0) then begin
types=replicate(rmdi,NumObs)
endif
goto, NOERROR
IOERROR: error=1
print,'read_cdfobs: an error occurred trying read files: ',files
NOERROR:
END