source: trunk/src/dmean_mld.pro @ 5

;+
;
; =============
; dmean_mld.pro
; =============
;
; .. function:: dmean_mld(data_prn)
;
;    Compute daily mean and standart deviation of MLD from structure ++details
;
;    :param data_prn: structure ++details
;    :type data_prn: structure
;    :raise data_prn: required
;
;    :returns: data structure +todo+ details or -1 if error
;    :rtype: structure +todo+ details
;
; :examples:
;
; Realistic example with POMME files::
;
;    IDL> data_prn=file_prn_to_mem()
;    IDL> data_dmean=dmean_mld(data_prn)
;
; impression de controle ::
;
;    IDL> help, data_dmean, /structure
;    IDL> print, data_dmean[0]
;    IDL> print, data_dmean.mean_zhom001
;    IDL> print, data_dmean.std_zhom001
;
; :uses:
;
; :see also:
;
; :func:`file_prn_to_mem`
; :func:`join_mean_asc`
;
; :restrictions:
;
; :todo:
;
; test arguments
;
; missing values
;
; init returned structure to NaN vs 0L or 0.0
;
; usage of report
;
; - fplod 20101129T092544Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * add 1st station in result structure
;
; - fplod 20101124T143429Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * creation
;
; :version:
;
; $Id$
;
;-
FUNCTION dmean_mld, data_prn
;
compile_opt idl2, strictarrsubs
;
; Return to caller if errors
ON_ERROR, 2
;
usage = 'result=dmean_mld(data_prn)'
;
nparam = N_PARAMS()
IF (nparam LT 1) THEN BEGIN
   ras = report(['Incorrect number of arguments.' $
         + '!C' $
         + 'Usage : ' + usage])
   return, -1
ENDIF
;
; initialize returned structure
resultstruct = { doy : 0L $
             , first_station : 0L $
             , mean_zhom001: 0.0 $
             , mean_zhom002 : 0.0 $
             , mean_zhom005 : 0.0 $
             , mean_zhom010: 0.0 $
             , mean_zhomelch : 0.0 $
             , std_zhom001: 0.0 $
             , std_zhom002 : 0.0 $
             , std_zhom005 : 0.0 $
             , std_zhom010: 0.0 $
             , std_zhomelch : 0.0 $
               }
;
; found days list
multiple_days=ceil(data_prn.doy)
days=multiple_days[UNIQ(multiple_days, SORT(multiple_days))]
ndays=n_elements(days)
;
; initialize dimension of the returned structure
result = replicate(resultstruct, ndays)
result.doy=days
;
; loop on each day
FOR iday=0L, ndays - 1 DO BEGIN
    ; find julian day in data_prn
    ;++print, result[iday].doy
    indice_prn = WHERE(ceil(data_prn.doy) EQ result[iday].doy,count)
    ;++print, count
    ;++print, 'indice du jour', result[iday].doy,' dans data_prn = ',indice_prn
    ; store the first (in numeric number - supposed to be chronological order)
    result[iday].first_station = data_prn[indice_prn[0]].station
    ; compute daily mean of MLDs
    result[iday].mean_zhom001 = total(data_prn[indice_prn].zhom001)/count
    result[iday].mean_zhom002 = total(data_prn[indice_prn].zhom002)/count
    result[iday].mean_zhom005 = total(data_prn[indice_prn].zhom005)/count
    result[iday].mean_zhom010 = total(data_prn[indice_prn].zhom010)/count
    result[iday].mean_zhomelch = total(data_prn[indice_prn].zhomelch)/count
    IF (count GT 1) THEN BEGIN
       ; compute standart deviation of MLDs
       result[iday].std_zhom001 = stddev(data_prn[indice_prn].zhom001)
       result[iday].std_zhom002 = stddev(data_prn[indice_prn].zhom002)
       result[iday].std_zhom005 = stddev(data_prn[indice_prn].zhom005)
       result[iday].std_zhom010 = stddev(data_prn[indice_prn].zhom010)
       result[iday].std_zhomelch = stddev(data_prn[indice_prn].zhomelch)
    ENDIF ELSE BEGIN
       result[iday].std_zhom001 = 0.0
       result[iday].std_zhom002 = 0.0
       result[iday].std_zhom005 = 0.0
       result[iday].std_zhom010 = 0.0
       result[iday].std_zhomelch = 0.0
    ENDELSE
    ;++print, 'mean mld 1', result[iday].mean_zhom001
ENDFOR
;
return, result
;
END