source: branches/UKMO/r6232_tracer_advection/NEMOGCM/TOOLS/OBSTOOLS/dataplot/read_cdfobs.pro @ 9295

;+---------------------------------------------------------------------------
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