Changeset 476 for trunk

Timestamp:

12/18/11 23:24:49 (12 years ago)

Author:

lelod

Message:

enlever ligne perimee mais encore bug

File:

• trunk/src/cresamsu.pro (modified) (2 diffs)

trunk/src/cresamsu.pro

@@ -224 +224 @@
 testfilename=''
 ; lecture des fichiers
-   for yyyy=yyyyb,yyyye do begin
-      for mm=mmb,mme do begin
-         for dd=ddb,dde do begin
+for yyyy=yyyyb,yyyye do begin
+   for mm=mmb,mme do begin
+      for dd=ddb,dde do begin
 ; ouverture fichier journalier
-            result=file_amsu_t2_to_mem( yyyy, mm, dd, numch $
-                                 , lon_min, lon_max, lat_min, lat_max $
-                                 , testfilename)
-            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
+         result=file_amsu_t2_to_mem( yyyy, mm, dd, numch $
+                                     , lon_min, lon_max, lat_min, lat_max $
+                                     , testfilename)
+         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)
+         nn=n_elements(tb1)
 ; egalement aussi des bornes en lon / lat (a retirer des parametres d'entree)++
 
 ; boucle sur les points du fichier
 
-            for i=0,nn do begin
+         for i=0,nn 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-lo1-rhdeg)/pash)>0
-                     kfn=fix((lon-lo1+rhdeg)/pash)+1<kmax
-                     ldp=fix((lat-la1-rhdeg)/pash)>0
-                     lfn=fix((lat-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
-                     if (fov gt fov_int2 and fov lt fov_int3) then begin
-                        yy=lat*coef
-                        xx=lon*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*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*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
-               endif            ; fin test donnees Nan
-            endfor
+                     endfor
+                  endfor
+               endif
+            endif
+         endfor
+         
            ; fin lecture donnees fichier journalier
-         endfor
       endfor
    endfor
+endfor
+
 ; fin boucle sur le temps
 
@@ -329 +329 @@
 jpie = n_elements(xxe)
 jpje = n_elements(yye)
+maskfin=fltarr(nblon,nblat)
 maskfin=interpolate(bate,indgen(nblon),indgen(nblat),/grid)
 ind=where(maskfin gt 0.5,nland)
-if nland ne 0 then maskfin(ind)=1
+if nland ne 0 then begin
+   maskfin(ind)=1
+endif
 ind=where(maskfin le 0.5,nsea)
-if nsea ne 0 then  maskfin(ind)=0
-
+if nsea ne 0 then begin
+   maskfin(ind)=0
+endif
 ;toto=reform(tb1fin[*,*,0])
                                 ;plt,xxlon,yylat,toto