source: trunk/src/correc_amsu.pro @ 550

; docformat = 'rst'
;+
;
; ============
; correc_amsu.pro
; ============
;
; .. function: correc_amsu(numch, yyyyb, mmb, ddb, yyyye, mme, dde, lonmin, lonmax, latmin, latmax)
; DESCRIPTION
; ===========
;
; programme d'interpolation avec algorithme de Cressmann'ajustement
;des donnees nadir (en complément de la correction de FK)
;
;; grille de sortie en long / lat
; recoit en entree les jours de l'annee jdeb et jfin
;correspondant au debut et la fin du mois
;

; EXAMPLES
; ========
;
; ::
;
;    numch='a5'
;    yyyyb=2006L
;    yyyye=2006L
;    mmb=8
;    mme=8
;    ddb=1
;    dde=3
;    lonmin=-25.
;    lonmax=25.
;    latmin=-5.
;    latmax=20.

;    result = cresamsu(numch, yyyyb, mmb, ddb $
;            , yyyye, mme, dde $
;            , lonmin, lonmax, latmin, latmax)
; 
;
; SEE ALSO
; ========
;
;  Previous step :ref:`extract_amsu.pro` (:ref:`traite_amsuab.sh`)
; :ref:`traite_amsuab.sh`
;
; Use :func:`file_amsu_t2_to_mem
; :func:`idl_amsu_netcdf`
;
; TODO
; ====
;
;
;; EVOLUTIONS
; ==========
;

function correc_amsu $
   , numch $
   , yyyyb $
   , mmb $
   , ddb $
   , yyyye $
   , mme $
   , dde $
   , lonmin $
   , lonmax $
   , latmin $
   , latmax

compile_opt idl2, strictarrsubs
@cm_project
;
s = size(SCOPE_TRACEBACK(/STRUCTURE),/DIMENSION)
routine = (SCOPE_TRACEBACK(/STRUCTURE))[s-1].Routine
;
key_performance=0
if key_performance EQ 1 THEN BEGIN
    time1 = SYSTIME(1)
ENDIF
;
;
; Return to caller if errors
ON_ERROR, 2
result = -1
;
usage = 'result = correc_amsu(numch, yyyyb, mmb, ddb, yyyye, mme, dde, lonmin, lonmax, latmin, latmax)'
;
nparam = N_PARAMS()
IF (nparam NE 11) THEN BEGIN
    ras = report(['Incorrect number of arguments.' $
          + '!C' $
          + 'Usage : ' + usage])
    return, result
ENDIF
;
tbmin=100
tbmax=350

prefix=numch
look = 'filename'
geomin = geolocation_to_string_idl(lonmin, latmin, look,1)
geomax = geolocation_to_string_idl(lonmax, latmax, look,1)

testfilename=''
nbfile = 0L
; lecture des fichiers
ddmb=ddb
ddme=dde
mmyb=mmb
mmye=mme
indic=0
for yyyy=yyyyb,yyyye do begin
   if (yyyye ne yyyyb) then begin
      mmyb=1
      mmye=12
      if yyyy eq yyyyb then mmyb=mmb
      if yyyy eq yyyye then mmye=mme
   endif
   print,"traitement annee ",yyyy,'mois ',mmyb,'a ',mmye
   for mm=mmyb,mmye do begin
      if (mme ne mmb) then begin
         ddmb=1
         ddme=31
         if mm eq mmb then ddmb=ddb
         if mm eq mme then ddme=dde
     endif
      print,"mois ",mm,"jours ",ddmb,"a ",ddme
      for dd=ddmb,ddme do begin
         print, 'iii : traitement du jour ',yyyy,mm,dd
         ; lecture du fichier journalier
         result=file_amsu_t2_to_mem( yyyy, mm, dd, numch $
                                     , lonmin, lonmax, latmin, latmax $
                                     , testfilename)
         if (size(result,/TYPE) EQ 3) THEN BEGIN
            goto, nofile
         ENDIF
         desc=result.data.desc
         hour=result.data.hour
         fov=result.data.fov
         lon=result.data.lon
         lat=result.data.lat
         mask=result.data.landseamask
         tb1=result.data.tb
            ; decodage du nb de points (nn)
         nn=n_elements(tb1)

         if indic eq 0 then begin
            nfov=max(fov)
            fovv=indgen(nfov)+1
            sumtb=fltarr(nfov)
            vartb=fltarr(nfov)
            kont=intarr(nfov)*0L
         endif
         nfovday=max(fov)
         if nfovday gt nfov then goto, nextday
         for ifov=0,nfov-1 do begin 
            ind=where(fov eq ifov+1 and mask eq 1 ,nifov) 
            if (nifov ne 0) then begin 
               sumtb[ifov]=sumtb[ifov]+tb1[ind] 
               vartb[ifov]=vartb[ifov]+tb1[ind]*tb1[ind] 
               kont[ifov]=kont[ifov]+1L
            endif 
         endfor 

nofile:
         print, 'www : no data to read for this day'
         goto, nextday
nextday:
         indic=indic+1
      endfor
   endfor
; fin boucle sur le temps
endfor

moytb=sumtb/kont
sigtb=sqrt((vartb-moytb*moytb)/kont)
window,0
plot,fovv,moytb,psym=7,xstyle=1,ystyle=1,yrange=[min(moytb)-max(sigtb),max(moytb)+max(sigtb)]
errplot,moytb-sigtb,moytb+sigtb


print,moytb

realend:
        
result = 0
return, result  
end