Changeset 574

Timestamp:

06/04/12 23:00:37 (12 years ago)

Author:

lelod

Message:

correction bug sur dates

File:

• trunk/src/cresamsu.pro (modified) (1 diff)

trunk/src/cresamsu.pro

@@ -355 +355 @@
 nbfile = 0L
 ; lecture des fichiers
-ddmb=ddb
-ddme=dde
-mmyb=mmb
-mmye=mme
-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
+;ddmb=ddb
+;ddme=dde
+;mmyb=mmb
+;mmye=mme
+;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
+dated=julday(mmb,ddb,yyyyb)
+datee=julday(mme,dde,yyyye)
+print,'dates debut et fin en jours juliens',dateb,datee
+for jd=dateb,datee do begin
+   caldat,jd, mm,dd,yyyy
+   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
+   tbinit=result.data.tb
+                                ; decodage du nb de points (nn)
+   nn=n_elements(tbinit)
+; egalement aussi des bornes en lon / lat (a retirer des parametres d'entree)++
+
+; boucle sur les points du fichier
+   tb1=tbinit
+   if correcbis eq "oui" then begin
+      print,"ATTENTION: correction supplementaire activee"
+; application de la correction bis
+      tb1=fltarr(nn)
+      print,'fauchee',min(fov),max(fov)
+      for ifov=min(fov)+1,fovmax-1 do begin
+         ind=where(fov eq ifov and mask eq 1,nbfov)
+         if nbfov ne 0 then begin
+            tb1[ind]=tbinit[ind] -correction_t[ifov]
+         endif
+      endfor
+      ind=where(fov eq 1,smallfov)
+      if smallfov ne 0 then begin
+         tb1[ind]=tbinit[ind] -correction_t[0]
+      endif
+      ind=where(fov ge fovmax,bigfov)
+      if bigfov ne 0 then begin
+         tb1[ind]=tbinit[ind] -correction_t[fovmax-1]
+      endif
+      for ifov=min(fov)+1,fovmax-1 do begin
+         ind=where(fov eq ifov and mask eq 0,nbfov)
+         if nbfov ne 0 then begin
+            tb1[ind]=tbinit[ind]-correction_m[ifov]
+         endif
+      endfor
+      ind=where(fov eq 1,smallfov)
+      if smallfov ne 0 then begin
+         tb1[ind]=tbinit[ind] -correction_m[0]
+      endif
+      ind=where(fov ge fovmax,bigfov)
+      if bigfov ne 0 then begin
+         tb1[ind]=tbinit[ind] -correction_m[fovmax-1]
+      endif
    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
-         tbinit=result.data.tb
-            ; decodage du nb de points (nn)
-         nn=n_elements(tbinit)
-; egalement aussi des bornes en lon / lat (a retirer des parametres d'entree)++
-
-; boucle sur les points du fichier
-         tb1=tbinit
-         if correcbis eq "oui" then begin
-             print,"ATTENTION: correction supplementaire activee"
-; application de la correction bis
-             tb1=fltarr(nn)
-             print,'fauchee',min(fov),max(fov)
-             for ifov=min(fov)+1,fovmax-1 do begin
-                 ind=where(fov eq ifov and mask eq 1,nbfov)
-                 if nbfov ne 0 then begin
-                     tb1[ind]=tbinit[ind] -correction_t[ifov]
-                 endif
-             endfor
-             ind=where(fov eq 1,smallfov)
-             if smallfov ne 0 then begin
-                 tb1[ind]=tbinit[ind] -correction_t[0]
-             endif
-             ind=where(fov ge fovmax,bigfov)
-             if bigfov ne 0 then begin
-                 tb1[ind]=tbinit[ind] -correction_t[fovmax-1]
-             endif
-             for ifov=min(fov)+1,fovmax-1 do begin
-                 ind=where(fov eq ifov and mask eq 0,nbfov)
-                 if nbfov ne 0 then begin
-                     tb1[ind]=tbinit[ind]-correction_m[ifov]
-                 endif
-             endfor
-             ind=where(fov eq 1,smallfov)
-             if smallfov ne 0 then begin
-                 tb1[ind]=tbinit[ind] -correction_m[0]
-             endif
-             ind=where(fov ge fovmax,bigfov)
-             if bigfov ne 0 then begin
-                 tb1[ind]=tbinit[ind] -correction_m[fovmax-1]
-             endif
-         endif
-
-         for i=0,nn-1 do begin
+   
+   for i=0,nn-1 do begin
 ; recup de tb1,hour,lon,lat
-            if (tb1[i] ge tbmin and tb1[i] le tbmax) then begin
+      if (tb1[i] ge tbmin and tb1[i] le tbmax) then begin
 ;if  ((tb1 lt tbmin) or (tb1 gt tbmax)) then goto,suite0 ;; on ne s interesse pas aux donnees aberrantes
 ; a remplacer par if tb1 ne Nan then begin
 ;if (cont lt 10) then print,lon,lat,fov,tb1,hour
 ; test sur le debut et la fin des donnes utiles tenant compte du rayon temporel
-               jour=julday(mm,dd,yyyy,0)*24.+hour[i] ; en heures decimales
+         jour=julday(mm,dd,yyyy,0)*24.+hour[i] ; en heures decimales
                                 ;if (cont lt 10) then print,mm,dd,yyyy,'jour julien lu (en heure)',jour
-               if ((jour ge hrdeb-rtemphor) and (jour le hrfin+rtemphor) and (hour[i] le hfin) and (hour[i] ge hdeb) ) then begin
+         if ((jour ge hrdeb-rtemphor) and (jour le hrfin+rtemphor) and (hour[i] le hfin) and (hour[i] ge hdeb) ) then begin
                                 ; print,'on cherche les donnees dans la fenetre temporelle'
                                 ; utilisation seulement des orbites du matin - projet heat low
                                 ; grilles temporelle et spatiale utile
-                  jdp=fix((jour-hrdeb-rtemphor)/pasthor)>0
-                  jfn=fix((jour-hrdeb+rtemphor)/pasthor)+1<jmax
+            jdp=fix((jour-hrdeb-rtemphor)/pasthor)>0
+            jfn=fix((jour-hrdeb+rtemphor)/pasthor)+1<jmax
                                 ;print,'grille utile temporelle',jdp,jfn
-                  kdp=fix((lon[i]-lo1-rhdeg)/pash)>0
-                  kfn=fix((lon[i]-lo1+rhdeg)/pash)+1<kmax
-                  ldp=fix((lat[i]-la1-rhdeg)/pash)>0
-                  lfn=fix((lat[i]-la1+rhdeg)/pash)+1<lmax
+            kdp=fix((lon[i]-lo1-rhdeg)/pash)>0
+            kfn=fix((lon[i]-lo1+rhdeg)/pash)+1<kmax
+            ldp=fix((lat[i]-la1-rhdeg)/pash)>0
+            lfn=fix((lat[i]-la1+rhdeg)/pash)+1<lmax
                                 ;print,'grille spatiale utile', kdp,kfn,ldp,lfn
                                 ; calcul de la fonction de Cressman
-                  yy=lat[i]*coef
-                  xx=lon[i]*coef
-                  coyy=cos(yy)
-                  coxx=cos(xx)
-                  sixx=sin(xx)
-                  siyy=sin(yy)
+            yy=lat[i]*coef
+            xx=lon[i]*coef
+            coyy=cos(yy)
+            coxx=cos(xx)
+            sixx=sin(xx)
+            siyy=sin(yy)
                                 ; boucle pour le calcul des contributions aux points de la grille
                                 ; (los, las, t)
-                  for nj=jdp,jfn do begin
-                     dij=(jour-t[nj])/rtemphor ; normalisee par le rayon de Cressmann
-                     disj=dij*dij
-                     for l=ldp,lfn do begin
-                        coco=coyy*cos(las[l])
-                        sisi=siyy*sin(las[l])
-                        for k=kdp,kfn do begin
-                           xxs=los[k]
-                           ang=coco*(coxx*cos(xxs)+sixx*sin(xxs))+sisi
-                           vv=ang
-                           ac=acos(vv)
-                           dist=rterre*rterre*ac*ac/r2
-                           if (dist lt 1. and disj lt 1.) then begin
-                              rnorm=2.
-                              poids=(rnorm-dist-disj)/(rnorm+dist+disj)
-                              stb1[k,l,nj]=stb1[k,l,nj]+tb1[i]*poids
-                              z1[k,l,nj]=z1[k,l,nj]+poids
-                              kont1[k,l,nj]=kont1[k,l,nj]+1
+            for nj=jdp,jfn do begin
+               dij=(jour-t[nj])/rtemphor ; normalisee par le rayon de Cressmann
+               disj=dij*dij
+               for l=ldp,lfn do begin
+                  coco=coyy*cos(las[l])
+                  sisi=siyy*sin(las[l])
+                  for k=kdp,kfn do begin
+                     xxs=los[k]
+                     ang=coco*(coxx*cos(xxs)+sixx*sin(xxs))+sisi
+                     vv=ang
+                     ac=acos(vv)
+                     dist=rterre*rterre*ac*ac/r2
+                     if (dist lt 1. and disj lt 1.) then begin
+                        rnorm=2.
+                        poids=(rnorm-dist-disj)/(rnorm+dist+disj)
+                        stb1[k,l,nj]=stb1[k,l,nj]+tb1[i]*poids
+                        z1[k,l,nj]=z1[k,l,nj]+poids
+                        kont1[k,l,nj]=kont1[k,l,nj]+1
                                 ; print,kont1[k,l,nj]
-                           endif
-                        endfor
-                     endfor
+                     endif
                   endfor
-               endif
-            endif
+               endfor
+            endfor
+         endif
+      endif
             ; fin test sur donnees valides (inutile normalement)
-         endfor
+   endfor
          ; fin lecture donnees fichier journalier
-         nbfile = nbfile + 1L
-         goto, nextday
+   nbfile = nbfile + 1L
+   goto, nextday
 nofile:
-    print, 'www : no data to read for this day'
-    goto, nextday
+   print, 'www : no data to read for this day'
+   goto, nextday
 
 nextday:
-       IF key_performance EQ 1 THEN BEGIN
-          msg = report(['ppp : ' + routine + ' : durée intermediaire au passage au fichier suivant' $
-                   + string(SYSTIME(1)-time1,format='(F12.6)')])
-
-       ENDIF
-      endfor
-   endfor
+   IF key_performance EQ 1 THEN BEGIN
+      msg = report(['ppp : ' + routine + ' : durée intermediaire au passage au fichier suivant' $
+                    + string(SYSTIME(1)-time1,format='(F12.6)')])
+
+   ENDIF
+endfor
+
 ; fin boucle sur le temps
-endfor
+
 
 if nbfile GT 0 THEN BEGIN