Index: /trunk/LMDZ.TITAN/Tools/angmom.F90
===================================================================
--- /trunk/LMDZ.TITAN/Tools/angmom.F90	(revision 1453)
+++ /trunk/LMDZ.TITAN/Tools/angmom.F90	(revision 1454)
@@ -47,11 +47,11 @@
 integer,dimension(4) :: datashape4d ! shape of 4D datasets
 
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
+real :: miss_val ! special "missing value" to specify missing data
+real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
 real :: pi
 real,dimension(:),allocatable :: lon ! longitude
 integer lonlength ! # of grid points along longitude
 real,dimension(:),allocatable :: lat ! latitude
-real,dimension(:),allocatable :: latrad ! latitude in radian
+real,dimension(:),allocatable :: coslat ! cos of latitude
 integer latlength ! # of grid points along latitude
 real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
@@ -117,4 +117,5 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
@@ -160,9 +161,10 @@
 if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
 
-allocate(latrad(latlength))
-latrad = lat*pi/180.
+allocate(coslat(latlength))
+! Beware of rounding problems at poles...
+coslat(:) = max(0.,cos(lat(:)*pi/180.))
 
 ! Lat, lon pressure intervals
-deltalat = abs(latrad(2)-latrad(1))
+deltalat = abs(lat(2)-lat(1))*pi/180.
 deltalon = abs(lon(2)-lon(1))*pi/180.
 
@@ -281,5 +283,5 @@
 endif
 
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get U ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
@@ -295,5 +297,5 @@
 
 !text="zareoid"
-!call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
+!call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
 !if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
 !if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
@@ -309,5 +311,5 @@
   text="DUVDF"
 endif
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duvdf,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duvdf,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "Failed to get duvdf ID"
@@ -357,5 +359,5 @@
 
 text="dugwo"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwo,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwo,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "Failed to get dugwo ID"
@@ -367,5 +369,5 @@
 
 text="dugwno"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwno,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwno,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "Failed to get dugwno ID"
@@ -390,5 +392,5 @@
 
 text="duajs"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duajs,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duajs,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "Failed to get duajs ID"
@@ -414,5 +416,5 @@
   text="DUDYN"
 endif
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dudyn,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dudyn,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "Failed to get dudyn ID"
@@ -486,6 +488,6 @@
 enddo ! timelength
 
-call cellmass(infid,latlength,lonlength,altlength,timelength, &
-              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
+call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
+              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
               dmass )
 
@@ -504,8 +506,8 @@
       osam(ilon,ilat,ilev,itim) = &
          rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim) &
-       * cos(latrad(ilat))*cos(latrad(ilat)) &
+       * coslat(ilat)*coslat(ilat) &
        * omega
       rsam(ilon,ilat,ilev,itim) = vitu(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)
     else
       osam(ilon,ilat,ilev,itim) = miss_val
@@ -569,5 +571,5 @@
        tmpp = (ps(ilon+1,ilat,itim)  +ps(ilon,ilat,itim))/2.
       endif
-      tmou(itim) = tmou(itim) + a0*a0* deltalat*cos(latrad(ilat)) &
+      tmou(itim) = tmou(itim) + a0*a0* deltalat*coslat(ilat) &
        * signe*dz*tmpp &
        / hadley
@@ -582,5 +584,5 @@
     if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
       tdyn(itim) = tdyn(itim) + dudyn(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
        * dmass(ilon,ilat,ilev,itim) &
        / hadley
@@ -599,5 +601,5 @@
     if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
       tajs(itim) = tajs(itim) + duajs(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
        * dmass(ilon,ilat,ilev,itim) &
        / hadley
@@ -616,5 +618,5 @@
     if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
       tbls(itim) = tbls(itim) + duvdf(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
        * dmass(ilon,ilat,ilev,itim) &
        / hadley
@@ -633,5 +635,5 @@
     if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
       tgwo(itim) = tgwo(itim) + dugwo(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
        * dmass(ilon,ilat,ilev,itim) &
        / hadley
@@ -650,5 +652,5 @@
     if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
       tgwno(itim) = tgwno(itim) + dugwno(ilon,ilat,ilev,itim) &
-         * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))       &
+         * rayon(ilon,ilat,ilev,itim)*coslat(ilat)       &
          * dmass(ilon,ilat,ilev,itim)  &
        / hadley
Index: /trunk/LMDZ.TITAN/Tools/compile_pgf
===================================================================
--- /trunk/LMDZ.TITAN/Tools/compile_pgf	(revision 1453)
+++ /trunk/LMDZ.TITAN/Tools/compile_pgf	(revision 1454)
@@ -1,5 +1,5 @@
 # path for netcdf should be adapted to your configuration !
 
-\cp -f $2.h planet.h
+#\cp -f $2.h planet.h
  pgf95 -Bstatic cpdet.F90 moyzon.F moyzon2.F moytim.F dx_dp.F epflux.F90 \
     io.F90 dmass.F90 reverse.F90 $1.F90 \
Index: /trunk/LMDZ.TITAN/Tools/compilefft_pgf
===================================================================
--- /trunk/LMDZ.TITAN/Tools/compilefft_pgf	(revision 1453)
+++ /trunk/LMDZ.TITAN/Tools/compilefft_pgf	(revision 1454)
@@ -2,6 +2,6 @@
 # path for fftw3  should be adapted to your configuration !
 
-\rm planet.h
-ln -s $2.h planet.h
+#\rm planet.h
+#ln -s $2.h planet.h
  pgf95 -Bstatic cpdet.F90 moyzon.F moyzon2.F moytim.F dx_dp.F epflux.F90 \
     io.F90 dmass.F90 reverse.F90 $1.F90 \
Index: /trunk/LMDZ.TITAN/Tools/dmass.F90
===================================================================
--- /trunk/LMDZ.TITAN/Tools/dmass.F90	(revision 1453)
+++ /trunk/LMDZ.TITAN/Tools/dmass.F90	(revision 1454)
@@ -1,4 +1,4 @@
-subroutine cellmass(infid,latlength,lonlength,altlength,timelength, &
-                     lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
+subroutine cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
+                     miss_val,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
                      dmass )
 
@@ -14,4 +14,5 @@
 integer timelength ! # of points along time
 logical lmdflag ! true=LMD, false=CAM
+real :: miss_val ! missing value
 real :: deltalat,deltalon ! lat and lon intervals in radians
 real,dimension(latlength) :: latrad ! latitudes in radians
@@ -32,6 +33,4 @@
 real :: p0 ! GCM reference pressure
 integer tmpvarid
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
 
 include "planet.h"
@@ -46,5 +45,5 @@
 if(lmdflag) then  ! LMD vs CAM
   text="pres"
-  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,press,ierr1,ierr2)
+  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,press,miss_val,ierr1,ierr2)
   if (ierr1.ne.NF_NOERR) then
 ! assume we are using _P files 
Index: /trunk/LMDZ.TITAN/Tools/dx_dp.F
===================================================================
--- /trunk/LMDZ.TITAN/Tools/dx_dp.F	(revision 1453)
+++ /trunk/LMDZ.TITAN/Tools/dx_dp.F	(revision 1454)
@@ -30,5 +30,5 @@
 c-----------------------------------------------------------------------
       do j=1,jjp1
-         if ((x(j,1).lt.indefini).and.(x(j,2).lt.indefini)) then  
+         if ((x(j,1).ne.indefini).and.(x(j,2).ne.indefini)) then  
            dxdp(j,1) = (x(j,2)-x(j,1))/(pniv(2) - pniv(1))
          else
@@ -37,9 +37,9 @@
         
          do l=2,llm-1
-           if ((x(j,l-1).lt.indefini).and.(x(j,l+1).lt.indefini))then
+           if ((x(j,l-1).ne.indefini).and.(x(j,l+1).ne.indefini))then
              dxdp(j,l)= (x(j,l+1)-x(j,l-1))/(pniv(l+1) - pniv(l-1))
-         else if((x(j,l+1).lt.indefini).and.(x(j,l).lt.indefini))then
+         else if((x(j,l+1).ne.indefini).and.(x(j,l).ne.indefini))then
              dxdp(j,l)= (x(j,l+1)-x(j,l))  /(pniv(l+1) - pniv(l))
-         else if((x(j,l-1).lt.indefini).and.(x(j,l).lt.indefini))then
+         else if((x(j,l-1).ne.indefini).and.(x(j,l).ne.indefini))then
              dxdp(j,l)= (x(j,l)-x(j,l-1))  /(pniv(l)   - pniv(l-1))
            else 
@@ -47,5 +47,5 @@
            end if
          end do
-         if ((x(j,llm).lt.indefini).and.(x(j,llm-1).lt.indefini)) then
+         if ((x(j,llm).ne.indefini).and.(x(j,llm-1).ne.indefini)) then
          dxdp(j,llm)= (x(j,llm)-x(j,llm-1))/(pniv(llm)-pniv(llm-1))
          else
@@ -58,5 +58,5 @@
       do j=1,jjp1
          do l = 1, llm
-           if ( (abs(dxdp(j,l)).gt.(indefini/100.)).and.
+           if ( (abs(dxdp(j,l)).gt.(abs(indefini)/100.)).and.
      .          (dxdp(j,l).ne.indefini)) then
               write(*,*) '----> j= ', j , ' l= ' ,  l
Index: /trunk/LMDZ.TITAN/Tools/energy.F90
===================================================================
--- /trunk/LMDZ.TITAN/Tools/energy.F90	(revision 1453)
+++ /trunk/LMDZ.TITAN/Tools/energy.F90	(revision 1454)
@@ -36,11 +36,11 @@
 integer,dimension(4) :: datashape4d ! shape of 4D datasets
 
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
+real :: miss_val ! special "missing value" to specify missing data
+real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
 real :: pi
 real,dimension(:),allocatable :: lon ! longitude
 integer lonlength ! # of grid points along longitude
 real,dimension(:),allocatable :: lat ! latitude
-real,dimension(:),allocatable :: latrad ! latitude in radian
+real,dimension(:),allocatable :: coslat ! cos of latitude
 integer latlength ! # of grid points along latitude
 real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
@@ -80,4 +80,5 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
@@ -123,9 +124,10 @@
 if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
 
-allocate(latrad(latlength))
-latrad = lat*pi/180.
+allocate(coslat(latlength))
+! Beware of rounding problems at poles...
+coslat(:) = max(0.,cos(lat(:)*pi/180.))
 
 ! Lat, lon pressure intervals
-deltalat = abs(latrad(2)-latrad(1))
+deltalat = abs(lat(2)-lat(1))*pi/180.
 deltalon = abs(lon(2)-lon(1))*pi/180.
 
@@ -193,9 +195,9 @@
   text="T"
 endif
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "  looking for t instead... "
   text="t"
-  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
   if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
 endif
@@ -216,5 +218,5 @@
   text="U"
 endif
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
@@ -228,5 +230,5 @@
   text="V"
 endif
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
@@ -242,5 +244,5 @@
 
 !text="zareoid"
-!call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
+!call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
 !if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
 !if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
@@ -276,6 +278,6 @@
 enddo ! timelength
 
-call cellmass(infid,latlength,lonlength,altlength,timelength, &
-              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
+call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
+              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
               dmass )
 
Index: /trunk/LMDZ.TITAN/Tools/epflux.F90
===================================================================
--- /trunk/LMDZ.TITAN/Tools/epflux.F90	(revision 1453)
+++ /trunk/LMDZ.TITAN/Tools/epflux.F90	(revision 1454)
@@ -1,3 +1,3 @@
-subroutine  epflux(iip1,jjp1,llm,indefini,rlatu,rbar &
+subroutine  epflux(iip1,jjp1,llm,indefini,latdeg,rbar &
                   ,teta,u3d,v3d,w3d,press &
                   ,epy2d,epz2d,divep2d,vtem2d,wtem2d,acc_meridien2d &
@@ -41,5 +41,5 @@
       
       integer :: iip1,jjp1,llm
-      real :: indefini,rlatu(jjp1),rbar(jjp1,llm)
+      real :: indefini,latdeg(jjp1),rbar(jjp1,llm)
       REAL :: teta(iip1,jjp1,llm)
       REAL :: u3d(iip1,jjp1,llm),v3d(iip1,jjp1,llm)
@@ -75,4 +75,5 @@
 !  Autres
 !  ------
+      real :: rlatu(jjp1) ! lat in radians (beware of rounding effects at poles)
       real :: pi
       REAL :: tetas(llm)
@@ -88,4 +89,7 @@
       jjm = jjp1-1
       pi = 2.*asin(1.)
+
+! To avoid rounding effects at poles
+      rlatu(:)=latdeg(:)*pi/180.00001
 
 !     Calcul de moyennes zonales
@@ -108,5 +112,5 @@
          n= 0 
          do j=2,jjm
-            if (tetabar(j,l).lt.indefini) then 
+            if (tetabar(j,l).ne.indefini) then 
                tetas(l) = tetas(l) + tetabar(j,l)
                n = n+1 
@@ -237,5 +241,5 @@
       do l=1,llm
 
-       if ((ubar(2,l).lt.indefini).and.(ubar(1,l).lt.indefini))then
+       if ((ubar(2,l).ne.indefini).and.(ubar(1,l).ne.indefini))then
       ducosdlat(1,l) = ( ubar(2,l)*cos(rlatu(2))-ubar(1,l)*cos(rlatu(1)) ) & 
                       /( rlatu(2) - rlatu(1)) 
@@ -244,5 +248,5 @@
        end if
         do j=2,jjm
-       if ((ubar(j+1,l).lt.indefini).and.(ubar(j-1,l).lt.indefini))then
+       if ((ubar(j+1,l).ne.indefini).and.(ubar(j-1,l).ne.indefini))then
       ducosdlat(j,l) = ( ubar(j+1,l)*cos(rlatu(j+1))-ubar(j-1,l)*cos(rlatu(j-1))) &
                       /( rlatu(j+1) - rlatu(j-1))
@@ -251,5 +255,5 @@
        end if
         enddo
-       if ((ubar(jjp1,l).lt.indefini).and.(ubar(jjm,l).lt.indefini))then
+       if ((ubar(jjp1,l).ne.indefini).and.(ubar(jjm,l).ne.indefini))then
       ducosdlat(jjp1,l) = ( ubar(jjp1,l)*cos(rlatu(jjp1))  &
                                -ubar(jjm,l)*cos(rlatu(jjm))  ) &
@@ -260,6 +264,6 @@
            
         do j=1,jjp1
-       if ((vptetapbar(j,l).lt.indefini).and. &
-           (dtetabardp(j,l).lt.indefini)) then
+       if ((vptetapbar(j,l).ne.indefini).and. &
+           (dtetabardp(j,l).ne.indefini)) then
            vz(j,l) = vptetapbar(j,l)/dtetabardp(j,l)
          wlat(j,l) = cos(rlatu(j))*vz(j,l)
@@ -284,8 +288,8 @@
 !     ------------
 
-            if (   (rbar(j,l).lt.indefini).and. &
-                (vpupbar(j,l).lt.indefini).and. &
-                (dubardp(j,l).lt.indefini).and. &
-                     (vz(j,l).lt.indefini)) then
+            if (   (rbar(j,l).ne.indefini).and. &
+                (vpupbar(j,l).ne.indefini).and. &
+                (dubardp(j,l).ne.indefini).and. &
+                     (vz(j,l).ne.indefini)) then
            
      epy2d(j,l) = rbar(j,l) * cos(rlatu(j))* ( vpupbar(j,l) &
@@ -303,8 +307,8 @@
 !     ------------
 
-            if (   (rbar(j,l).lt.indefini).and. &
-                (wpupbar(j,l).lt.indefini).and. &
-              (ducosdlat(j,l).lt.indefini).and. &
-                     (vz(j,l).lt.indefini)) then
+            if (   (rbar(j,l).ne.indefini).and. &
+                (wpupbar(j,l).ne.indefini).and. &
+              (ducosdlat(j,l).ne.indefini).and. &
+                     (vz(j,l).ne.indefini)) then
 
      epz2d(j,l) = rbar(j,l) * cos(rlatu(j))*                    &
@@ -330,6 +334,6 @@
       do l=1,llm
 
-          if ( (wlat(2,l).lt.indefini) .and. &
-               (wlat(1,l).lt.indefini) ) then
+          if ( (wlat(2,l).ne.indefini) .and. &
+               (wlat(1,l).ne.indefini) ) then
                dwlatdlat(1,l)=(wlat(2,l)-wlat(1,l)) &
                      /(rlatu(2)-rlatu(1))
@@ -338,6 +342,6 @@
           end if
         do j=2,jjm
-          if ( (wlat(j+1,l).lt.indefini) .and. &
-               (wlat(j-1,l).lt.indefini) ) then
+          if ( (wlat(j+1,l).ne.indefini) .and. &
+               (wlat(j-1,l).ne.indefini) ) then
                dwlatdlat(j,l)=(wlat(j+1,l)-wlat(j-1,l)) &
                      /(rlatu(j+1)-rlatu(j-1))
@@ -346,6 +350,6 @@
           end if
         enddo 
-          if ( (wlat(jjp1,l).lt.indefini) .and. & 
-               (wlat(jjm, l).lt.indefini) ) then
+          if ( (wlat(jjp1,l).ne.indefini) .and. & 
+               (wlat(jjm, l).ne.indefini) ) then
                dwlatdlat(jjp1,l)=(wlat(jjp1,l)-wlat(jjm,l)) &
                      /(rlatu(jjp1)-rlatu(jjm))
@@ -355,6 +359,6 @@
 
         do j=1,jjp1
-          if ( (dvzdp(j,l).lt.indefini) &
-           .and.(vbar(j,l).lt.indefini)) then
+          if ( (dvzdp(j,l).ne.indefini) &
+           .and.(vbar(j,l).ne.indefini)) then
             vtem2d(j,l) = vbar(j,l)-dvzdp(j,l)
           else
@@ -364,7 +368,7 @@
 
         do j=1,jjp1
-          if (   (rbar(j,l).lt.indefini).and. &
-            (dwlatdlat(j,l).lt.indefini).and. &
-                 (wbar(j,l).lt.indefini)) then
+          if (   (rbar(j,l).ne.indefini).and. &
+            (dwlatdlat(j,l).ne.indefini).and. &
+                 (wbar(j,l).ne.indefini)) then
            wtem2d(j,l) = wbar(j,l)+dwlatdlat(j,l)/(rbar(j,l)*cos(rlatu(j)))
           else
@@ -385,9 +389,9 @@
       do l=1,llm
         do j=1,jjp1
-          if (   (rbar(j,l).lt.indefini) &
-          .and.(vtem2d(j,l).lt.indefini) &
-          .and.(wtem2d(j,l).lt.indefini) &
-       .and.(ducosdlat(j,l).lt.indefini) &
-         .and.(dubardp(j,l).lt.indefini)) then
+          if (   (rbar(j,l).ne.indefini) &
+          .and.(vtem2d(j,l).ne.indefini) &
+          .and.(wtem2d(j,l).ne.indefini) &
+       .and.(ducosdlat(j,l).ne.indefini) &
+         .and.(dubardp(j,l).ne.indefini)) then
        acc_meridien2d(j,l) =  &
            vtem2d(j,l)*(ducosdlat(j,l)/(rbar(j,l)*cos(rlatu(j))) &
@@ -410,8 +414,8 @@
         divep2d(1,l) =0
         do j=2,jjm
-          if (    (rbar(j,l).lt.indefini) .and. &
-               (epy2d(j+1,l).lt.indefini) .and. &
-               (epy2d(j-1,l).lt.indefini) .and. &
-                (depzdp(j,l).lt.indefini)    ) then
+          if (    (rbar(j,l).ne.indefini) .and. &
+               (epy2d(j+1,l).ne.indefini) .and. &
+               (epy2d(j-1,l).ne.indefini) .and. &
+                (depzdp(j,l).ne.indefini)    ) then
 
            depycosdlat=(epy2d(j+1,l)*cos(rlatu(j+1))   &
Index: /trunk/LMDZ.TITAN/Tools/fft.F90
===================================================================
--- /trunk/LMDZ.TITAN/Tools/fft.F90	(revision 1453)
+++ /trunk/LMDZ.TITAN/Tools/fft.F90	(revision 1454)
@@ -63,11 +63,10 @@
 integer,dimension(4) :: datashape4d ! shape of 4D datasets
 
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
+real :: miss_val ! special "missing value" to specify missing data
+real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
 real :: pi
 real,dimension(:),allocatable :: lon ! longitude
 integer lonlength ! # of grid points along longitude
 real,dimension(:),allocatable :: lat ! latitude
-real,dimension(:),allocatable :: latrad ! latitude in radian
 integer latlength ! # of grid points along latitude
 real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
@@ -125,6 +124,11 @@
 logical :: lmdflag
 
+! Tuning parameters
+real :: fcoup1,fcoup2,width
+real :: fcoup1tmp,fcoup2tmp,widthtmp
+logical,dimension(4) :: ok_out
+character (len=1) :: ok_outtmp
+
 include "planet.h"
-include "filter.h"
 
 #include <fftw3.f> 
@@ -135,7 +139,58 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
 write(*,*) "You are working on the atmosphere of ",planet
+
+! initialisation
+!----------------
+
+! Par defaut
+
+! Define the filters
+! Low  cutting frequency, in Hz    : fcoup1
+fcoup1=2.5e-6
+! High cutting frequency, in Hz    : fcoup2
+fcoup2=6.5e-6
+! Half-width of the filters, in Hz : width
+width=4.e-7
+! Outputs (U,     V,      W,     T)
+ok_out=(/.true.,.true.,.false.,.true./)
+
+print*,"Low  cutting frequency, in Hz ?"
+print*,"between 1e-5 and 1e-7, 0 for default => 2.5e-6"
+read(*,*) fcoup1tmp
+if ((fcoup1tmp.lt.1e-5).and.(fcoup1tmp.gt.1e-7)) fcoup1=fcoup1tmp
+print*,"=",fcoup1
+print*,"High cutting frequency, in Hz ?"
+print*,"between 1e-5 and 1e-7, 0 for default => 6.5e-6"
+read(*,*) fcoup2tmp
+if ((fcoup2tmp.lt.1e-5).and.(fcoup2tmp.gt.1e-7)) fcoup2=fcoup2tmp
+print*,"=",fcoup2
+print*,"Half-width of the filters, in Hz ?"
+print*,"between 1e-6 and 1e-8, 0 for default => 4e-7)"
+read(*,*) widthtmp
+if ((widthtmp.lt.1e-6).and.(widthtmp.gt.1e-8)) width=widthtmp
+print*,"=",width
+!width = 3./time(timelength)
+
+! Outputs 
+print*,"Output of zonal wind ? (y or n, default is y)"
+read(*,'(a1)') ok_outtmp
+if (ok_outtmp.eq."n") ok_out(1)=.false.
+print*,"=",ok_out(1)
+print*,"Output of meridional wind ? (y or n, default is y)"
+read(*,'(a1)') ok_outtmp
+if (ok_outtmp.eq."n") ok_out(2)=.false.
+print*,"=",ok_out(2)
+print*,"Output of vertical wind ? (y or n, default is n)"
+read(*,'(a1)') ok_outtmp
+if (ok_outtmp.eq."y") ok_out(3)=.true.
+print*,"=",ok_out(3)
+print*,"Output of temperature ? (y or n, default is y)"
+read(*,'(a1)') ok_outtmp
+if (ok_outtmp.eq."n") ok_out(4)=.false.
+print*,"=",ok_out(4)
 
 !===============================================================================
@@ -172,7 +227,4 @@
 ierr=NF_GET_VAR_REAL(infid,lat_varid,lat)
 if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
-
-allocate(latrad(latlength))
-latrad = lat*pi/180.
 
 allocate(plev(altlength))
@@ -212,9 +264,9 @@
 
 text="temp"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "  looking for t instead... "
   text="t"
-  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
   if (ierr1.ne.NF_NOERR) then
     print*,"Error: Failed to get temperature ID"
@@ -236,5 +288,5 @@
 
 text="vitu"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   print*,"Error: Failed to get vitu ID"
@@ -252,5 +304,5 @@
 
 text="vitv"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   print*,"Error: Failed to get vitv ID"
@@ -268,5 +320,5 @@
 
 text="vitw"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   print*,"Error: Failed to get vitw ID"
@@ -288,4 +340,5 @@
 ! 2.2.1 FFT and filtering
 !===============================================================================
+
 ! allocations
 !-------------
@@ -405,19 +458,4 @@
 enddo
 
-! Define the filters
-
-! IN filter.h :
-! Low  cutting frequency, in Hz    : fcoup1
-! High cutting frequency, in Hz    : fcoup2
-! Half-width of the filters, in Hz : width
-
-!print*,"Low  cutting frequency, in Hz ?"
-!read(*,*) fcoup1
-!print*,"High cutting frequency, in Hz ?"
-!read(*,*) fcoup2
-!print*,"Half-width of the filters, in Hz ?"
-!read(*,*) width
-!width = 3./time(timelength)
-
 do itim=1,M_fft+1
   if (freq(itim).lt.(fcoup1-width)) then
Index: /trunk/LMDZ.TITAN/Tools/io.F90
===================================================================
--- /trunk/LMDZ.TITAN/Tools/io.F90	(revision 1453)
+++ /trunk/LMDZ.TITAN/Tools/io.F90	(revision 1454)
@@ -218,5 +218,5 @@
 !===========================================================================
 
-subroutine get_var4d(infid,dim1,dim2,dim3,dim4,text,var,ierr1,ierr2)
+subroutine get_var4d(infid,dim1,dim2,dim3,dim4,text,var,missing,ierr1,ierr2)
 
 implicit none
@@ -229,5 +229,6 @@
 character (len=64) :: text ! name of variable to read
 real,dimension(dim1,dim2,dim3,dim4) :: var ! variable to read
-integer :: ierr1,ierr2 ! NetCDF routines return code
+real :: missing ! missing value
+integer :: ierr1,ierr2,miss ! NetCDF routines return code
 
 ! local
@@ -240,4 +241,5 @@
   write(*,*) "ID ok for ",trim(text)
   ierr2=NF_GET_VAR_REAL(infid,tmpvarid,var)
+  miss=NF_GET_ATT_REAL(infid,tmpvarid,"missing_value",missing)
 endif
 
Index: /trunk/LMDZ.TITAN/Tools/moytim.F
===================================================================
--- /trunk/LMDZ.TITAN/Tools/moytim.F	(revision 1453)
+++ /trunk/LMDZ.TITAN/Tools/moytim.F	(revision 1454)
@@ -33,5 +33,5 @@
            n = 0
            do t=1,nt
-              if (x(i,j,l,t).lt.indefini) then 
+              if (x(i,j,l,t).ne.indefini) then 
                  xmean(i,j,l) = xmean(i,j,l) + x(i,j,l,t)
                  n = n+1
Index: /trunk/LMDZ.TITAN/Tools/moyzon.F
===================================================================
--- /trunk/LMDZ.TITAN/Tools/moyzon.F	(revision 1453)
+++ /trunk/LMDZ.TITAN/Tools/moyzon.F	(revision 1454)
@@ -32,5 +32,5 @@
            n = 0
            do i=1,iim
-              if (x(i,j,l).lt.indefini) then 
+              if (x(i,j,l).ne.indefini) then 
                  xbar(j,l) = xbar(j,l) + x(i,j,l)
                  n = n+1
Index: /trunk/LMDZ.TITAN/Tools/moyzon2.F
===================================================================
--- /trunk/LMDZ.TITAN/Tools/moyzon2.F	(revision 1453)
+++ /trunk/LMDZ.TITAN/Tools/moyzon2.F	(revision 1454)
@@ -34,5 +34,5 @@
            n = 0
            do i=1,iim
-              if (x(i,j,l).lt.indefini) then 
+              if (x(i,j,l).ne.indefini) then 
                  xbar(j,l) = xbar(j,l) + x(i,j,l)
                  n = n+1
Index: /trunk/LMDZ.TITAN/Tools/psi.F90
===================================================================
--- /trunk/LMDZ.TITAN/Tools/psi.F90	(revision 1453)
+++ /trunk/LMDZ.TITAN/Tools/psi.F90	(revision 1454)
@@ -31,11 +31,11 @@
 integer,dimension(4) :: datashape4d ! shape of 4D datasets
 
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
+real :: miss_val ! special "missing value" to specify missing data
+real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
 real :: pi
 real,dimension(:),allocatable :: lon ! longitude
 integer lonlength ! # of grid points along longitude
 real,dimension(:),allocatable :: lat ! latitude
-real,dimension(:),allocatable :: latrad ! latitude in radian
+real,dimension(:),allocatable :: coslat ! cos of latitude
 integer latlength ! # of grid points along latitude
 real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
@@ -66,4 +66,5 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
@@ -104,9 +105,10 @@
 if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
 
-allocate(latrad(latlength))
-latrad = lat*pi/180.
+allocate(coslat(latlength))
+! Beware of rounding problems at poles...
+coslat(:) = max(0.,cos(lat(:)*pi/180.))
 
 ! Lat, lon pressure intervals
-deltalat = abs(latrad(2)-latrad(1))
+deltalat = abs(lat(2)-lat(1))*pi/180.
 deltalon = abs(lon(2)-lon(1))*pi/180.
 
@@ -160,5 +162,5 @@
 ! meridional wind vitv (in m/s)
 text="vitv"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
@@ -173,9 +175,9 @@
 
 text="zareoid"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "  looking for geop instead... "
   text="geop"
-  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
+  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
   if (ierr1.ne.NF_NOERR) stop "Error: Failed to get geop ID"
 do itim=1,timelength
@@ -225,6 +227,6 @@
 enddo ! timelength
 
-call cellmass(infid,latlength,lonlength,altlength,timelength, &
-              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
+call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
+              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
               dmass )
 
Index: /trunk/LMDZ.TITAN/Tools/stability.F90
===================================================================
--- /trunk/LMDZ.TITAN/Tools/stability.F90	(revision 1453)
+++ /trunk/LMDZ.TITAN/Tools/stability.F90	(revision 1454)
@@ -31,11 +31,11 @@
 integer,dimension(4) :: datashape4d ! shape of 4D datasets
 
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
+real :: miss_val ! special "missing value" to specify missing data
+real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
 real :: pi
 real,dimension(:),allocatable :: lon ! longitude
 integer lonlength ! # of grid points along longitude
 real,dimension(:),allocatable :: lat ! latitude
-real,dimension(:),allocatable :: latrad ! latitude in radian
+real,dimension(:),allocatable :: coslat ! cos of latitude
 integer latlength ! # of grid points along latitude
 real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
@@ -91,4 +91,5 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
@@ -129,9 +130,10 @@
 if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
 
-allocate(latrad(latlength))
-latrad = lat*pi/180.
+allocate(coslat(latlength))
+! Beware of rounding problems at poles...
+coslat(:) = max(0.,cos(lat(:)*pi/180.))
 
 ! Lat, lon pressure intervals
-deltalat = abs(latrad(2)-latrad(1))
+deltalat = abs(lat(2)-lat(1))*pi/180.
 deltalon = abs(lon(2)-lon(1))*pi/180.
 
@@ -180,9 +182,9 @@
 
 text="temp"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "  looking for t instead... "
   text="t"
-  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
   if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
 endif
@@ -197,5 +199,5 @@
 ! zonal wind vitu (in m/s)
 text="vitu"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
@@ -203,5 +205,5 @@
 ! meridional wind vitv (in m/s)
 text="vitv"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
@@ -215,5 +217,5 @@
 
 ! text="zareoid"
-! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
+! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
 ! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
 ! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
@@ -276,6 +278,6 @@
  grav(lonlength+1,:,:,:) =  grav(1,:,:,:)
 
-call cellmass(infid,latlength,lonlength+1,altlength,timelength, &
-              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
+call cellmass(infid,latlength,lonlength+1,altlength,timelength,lmdflag, &
+              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
               dmass )
 
@@ -364,6 +366,6 @@
 
 do ilat=2,latlength-1
-   if (tan(latrad(ilat)).ne.0.) then
-     fac1 = R0/tan(latrad(ilat))
+   if (tan(lat(ilat)*pi/180.).ne.0.) then
+     fac1 = R0/tan(lat(ilat)*pi/180.)
    else
      fac1 = miss_val
Index: /trunk/LMDZ.TITAN/Tools/tem.F90
===================================================================
--- /trunk/LMDZ.TITAN/Tools/tem.F90	(revision 1453)
+++ /trunk/LMDZ.TITAN/Tools/tem.F90	(revision 1454)
@@ -37,11 +37,11 @@
 integer,dimension(3) :: datashape3d ! shape of 3D datasets
 
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
+real :: miss_val ! special "missing value" to specify missing data
+real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
 real :: pi
 real,dimension(:),allocatable :: lon ! longitude
 integer lonlength ! # of grid points along longitude
 real,dimension(:),allocatable :: lat ! latitude
-real,dimension(:),allocatable :: latrad ! latitude in radian
+real,dimension(:),allocatable :: coslat ! cos of latitude
 integer latlength ! # of grid points along latitude
 real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
@@ -107,4 +107,5 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
@@ -145,9 +146,10 @@
 if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
 
-allocate(latrad(latlength))
-latrad = lat*pi/180.
+allocate(coslat(latlength))
+! Beware of rounding problems at poles...
+coslat(:) = max(0.,cos(lat(:)*pi/180.))
 
 ! Lat, lon pressure intervals
-deltalat = abs(latrad(2)-latrad(1))
+deltalat = abs(lat(2)-lat(1))*pi/180.
 deltalon = abs(lon(2)-lon(1))*pi/180.
 
@@ -196,9 +198,9 @@
 
 text="temp"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "  looking for t instead... "
   text="t"
-  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
   if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
 endif
@@ -214,5 +216,5 @@
 ! zonal wind vitu (in m/s)
 text="vitu"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
@@ -220,5 +222,5 @@
 ! meridional wind vitv (in m/s)
 text="vitv"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
@@ -226,5 +228,5 @@
 ! vertical wind vitw (in Pa/s)
 text="vitw"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitw ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading vertical wind"
@@ -238,5 +240,5 @@
 
 ! text="zareoid"
-! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
+! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
 ! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
 ! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
@@ -312,6 +314,6 @@
  grav(lonlength+1,:,:,:) =  grav(1,:,:,:)
 
-call cellmass(infid,latlength,lonlength+1,altlength,timelength, &
-              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
+call cellmass(infid,latlength,lonlength+1,altlength,timelength,lmdflag, &
+              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
               dmass )
 
@@ -359,5 +361,5 @@
 call moyzon(lonlength,latlength,altlength,miss_val,u3d,ubar)
 
-call epflux(lonlength+1,latlength,altlength,miss_val,latrad,rbar &
+call epflux(lonlength+1,latlength,altlength,miss_val,lat,rbar &
             ,teta,u3d,v3d,w3d,plev &
             ,epy2d,epz2d,div2d,vtem2d,wtem2d,ammc2d &
@@ -376,5 +378,5 @@
 ! rsurg: r*dp/g = dm/(r cos(lat) dlon dlat) !!!
      rsurg(ilon,ilat,ilev) = dmass(ilon,ilat,ilev,itim) &
-          / (r3d(ilon,ilat,ilev)*cos(latrad(ilat))*deltalon*deltalat)
+          / (r3d(ilon,ilat,ilev)*coslat(ilat)*deltalon*deltalat)
     else
      rsurg(ilon,ilat,ilev) = miss_val
@@ -391,5 +393,5 @@
         (rsurgbar(ilat,ilev).ne.miss_val) ) then
       vm(ilat,ilev) = vtem2d(ilat,ilev) &
-            * 2.*pi*rsurgbar(ilat,ilev)*cos(latrad(ilat))
+            * 2.*pi*rsurgbar(ilat,ilev)*coslat(ilat)
     else
       vm(ilat,ilev) = miss_val
Index: /trunk/LMDZ.TITAN/Tools/tmc.F90
===================================================================
--- /trunk/LMDZ.TITAN/Tools/tmc.F90	(revision 1453)
+++ /trunk/LMDZ.TITAN/Tools/tmc.F90	(revision 1454)
@@ -49,11 +49,11 @@
 integer,dimension(4) :: datashape4d ! shape of 4D datasets
 
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
+real :: miss_val ! special "missing value" to specify missing data
+real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
 real :: pi
 real,dimension(:),allocatable :: lon ! longitude
 integer lonlength ! # of grid points along longitude
 real,dimension(:),allocatable :: lat ! latitude
-real,dimension(:),allocatable :: latrad ! latitude in radian
+real,dimension(:),allocatable :: coslat ! cos of latitude
 integer latlength ! # of grid points along latitude
 real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
@@ -141,4 +141,5 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
@@ -179,9 +180,10 @@
 if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
 
-allocate(latrad(latlength))
-latrad = lat*pi/180.
+allocate(coslat(latlength))
+! Beware of rounding problems at poles...
+coslat(:) = max(0.,cos(lat(:)*pi/180.))
 
 ! Lat, lon pressure intervals
-deltalat = abs(latrad(2)-latrad(1))
+deltalat = abs(lat(2)-lat(1))*pi/180.
 deltalon = abs(lon(2)-lon(1))*pi/180.
 
@@ -233,5 +235,5 @@
 ! zonal wind vitu (in m/s)
 text="vitu"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
@@ -239,5 +241,5 @@
 ! meridional wind vitv (in m/s)
 text="vitv"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
@@ -245,5 +247,5 @@
 ! vertical wind vitw (in Pa/s)
 text="vitw"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitw ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading vertical wind"
@@ -257,5 +259,5 @@
 
 ! text="zareoid"
-! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
+! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
 ! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
 ! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
@@ -293,6 +295,6 @@
 enddo ! timelength
 
-call cellmass(infid,latlength,lonlength,altlength,timelength, &
-              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
+call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
+              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
               dmass )
 
@@ -316,8 +318,8 @@
       osam(ilon,ilat,ilev,itim) = &
          rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim) &
-       * cos(latrad(ilat))*cos(latrad(ilat)) &
+       * coslat(ilat)*coslat(ilat) &
        * omega
       rsam(ilon,ilat,ilev,itim) = vitu(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)
       tsam(ilon,ilat,ilev,itim) = osam(ilon,ilat,ilev,itim)&
                                 + rsam(ilon,ilat,ilev,itim)
Index: /trunk/LMDZ.VENUS/Tools/angmom.F90
===================================================================
--- /trunk/LMDZ.VENUS/Tools/angmom.F90	(revision 1453)
+++ /trunk/LMDZ.VENUS/Tools/angmom.F90	(revision 1454)
@@ -47,11 +47,11 @@
 integer,dimension(4) :: datashape4d ! shape of 4D datasets
 
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
+real :: miss_val ! special "missing value" to specify missing data
+real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
 real :: pi
 real,dimension(:),allocatable :: lon ! longitude
 integer lonlength ! # of grid points along longitude
 real,dimension(:),allocatable :: lat ! latitude
-real,dimension(:),allocatable :: latrad ! latitude in radian
+real,dimension(:),allocatable :: coslat ! cos of latitude
 integer latlength ! # of grid points along latitude
 real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
@@ -117,4 +117,5 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
@@ -160,9 +161,10 @@
 if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
 
-allocate(latrad(latlength))
-latrad = lat*pi/180.
+allocate(coslat(latlength))
+! Beware of rounding problems at poles...
+coslat(:) = max(0.,cos(lat(:)*pi/180.))
 
 ! Lat, lon pressure intervals
-deltalat = abs(latrad(2)-latrad(1))
+deltalat = abs(lat(2)-lat(1))*pi/180.
 deltalon = abs(lon(2)-lon(1))*pi/180.
 
@@ -281,5 +283,5 @@
 endif
 
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get U ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
@@ -295,5 +297,5 @@
 
 !text="zareoid"
-!call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
+!call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
 !if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
 !if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
@@ -309,5 +311,5 @@
   text="DUVDF"
 endif
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duvdf,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duvdf,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "Failed to get duvdf ID"
@@ -357,5 +359,5 @@
 
 text="dugwo"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwo,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwo,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "Failed to get dugwo ID"
@@ -367,5 +369,5 @@
 
 text="dugwno"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwno,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dugwno,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "Failed to get dugwno ID"
@@ -390,5 +392,5 @@
 
 text="duajs"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duajs,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duajs,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "Failed to get duajs ID"
@@ -414,5 +416,5 @@
   text="DUDYN"
 endif
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dudyn,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dudyn,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "Failed to get dudyn ID"
@@ -486,6 +488,6 @@
 enddo ! timelength
 
-call cellmass(infid,latlength,lonlength,altlength,timelength, &
-              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
+call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
+              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
               dmass )
 
@@ -504,8 +506,8 @@
       osam(ilon,ilat,ilev,itim) = &
          rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim) &
-       * cos(latrad(ilat))*cos(latrad(ilat)) &
+       * coslat(ilat)*coslat(ilat) &
        * omega
       rsam(ilon,ilat,ilev,itim) = vitu(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)
     else
       osam(ilon,ilat,ilev,itim) = miss_val
@@ -569,5 +571,5 @@
        tmpp = (ps(ilon+1,ilat,itim)  +ps(ilon,ilat,itim))/2.
       endif
-      tmou(itim) = tmou(itim) + a0*a0* deltalat*cos(latrad(ilat)) &
+      tmou(itim) = tmou(itim) + a0*a0* deltalat*coslat(ilat) &
        * signe*dz*tmpp &
        / hadley
@@ -582,5 +584,5 @@
     if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
       tdyn(itim) = tdyn(itim) + dudyn(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
        * dmass(ilon,ilat,ilev,itim) &
        / hadley
@@ -599,5 +601,5 @@
     if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
       tajs(itim) = tajs(itim) + duajs(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
        * dmass(ilon,ilat,ilev,itim) &
        / hadley
@@ -616,5 +618,5 @@
     if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
       tbls(itim) = tbls(itim) + duvdf(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
        * dmass(ilon,ilat,ilev,itim) &
        / hadley
@@ -633,5 +635,5 @@
     if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
       tgwo(itim) = tgwo(itim) + dugwo(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))      &
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)      &
        * dmass(ilon,ilat,ilev,itim) &
        / hadley
@@ -650,5 +652,5 @@
     if (rayon(ilon,ilat,ilev,itim).ne.miss_val) then
       tgwno(itim) = tgwno(itim) + dugwno(ilon,ilat,ilev,itim) &
-         * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))       &
+         * rayon(ilon,ilat,ilev,itim)*coslat(ilat)       &
          * dmass(ilon,ilat,ilev,itim)  &
        / hadley
Index: /trunk/LMDZ.VENUS/Tools/compilefft_pgf
===================================================================
--- /trunk/LMDZ.VENUS/Tools/compilefft_pgf	(revision 1453)
+++ /trunk/LMDZ.VENUS/Tools/compilefft_pgf	(revision 1454)
@@ -2,6 +2,6 @@
 # path for fftw3  should be adapted to your configuration !
 
-\rm planet.h
-ln -s $2.h planet.h
+#\rm planet.h
+#ln -s $2.h planet.h
  pgf95 -Bstatic cpdet.F90 moyzon.F moyzon2.F moytim.F dx_dp.F epflux.F90 \
     io.F90 dmass.F90 reverse.F90 $1.F90 \
Index: /trunk/LMDZ.VENUS/Tools/dmass.F90
===================================================================
--- /trunk/LMDZ.VENUS/Tools/dmass.F90	(revision 1453)
+++ /trunk/LMDZ.VENUS/Tools/dmass.F90	(revision 1454)
@@ -1,4 +1,4 @@
-subroutine cellmass(infid,latlength,lonlength,altlength,timelength, &
-                     lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
+subroutine cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
+                     miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
                      dmass )
 
@@ -14,6 +14,7 @@
 integer timelength ! # of points along time
 logical lmdflag ! true=LMD, false=CAM
+real :: miss_val ! missing value
 real :: deltalat,deltalon ! lat and lon intervals in radians
-real,dimension(latlength) :: latrad ! latitudes in radians
+real,dimension(latlength) :: coslat ! cos of latitudes
 real,dimension(altlength) :: plev ! Pressure levels (Pa)
 real,dimension(lonlength,latlength,timelength),intent(in) :: ps ! surface pressure
@@ -32,6 +33,4 @@
 real :: p0 ! GCM reference pressure
 integer tmpvarid
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
 
 include "planet.h"
@@ -46,5 +45,5 @@
 if(lmdflag) then  ! LMD vs CAM
   text="pres"
-  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,press,ierr1,ierr2)
+  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,press,miss_val,ierr1,ierr2)
   if (ierr1.ne.NF_NOERR) then
 ! assume we are using _P files 
@@ -119,5 +118,12 @@
  do ilat=1,latlength
   tmpp=ps(ilon,ilat,itim)
-  pp=press(ilon,ilat,:,itim)
+! Beware when miss_val is negative...
+  do ilev=1,altlength
+    if (press(ilon,ilat,ilev,itim).eq.miss_val) then
+      pp(ilev)=1.e33
+    else
+      pp(ilev)=press(ilon,ilat,ilev,itim)
+    endif
+  enddo
 if (tmpp.ge.pp(1)) then
   deltap(ilon,ilat,1,itim)= tmpp - exp((log(pp(1))+log(pp(2)))/2.) ! initialization, rectified with ps below
@@ -136,5 +142,12 @@
  do ilat=1,latlength
   tmpp=ps(ilon,ilat,itim)
-  pp=press(ilon,ilat,:,itim)
+! Beware when miss_val is negative...
+  do ilev=1,altlength
+    if (press(ilon,ilat,ilev,itim).eq.miss_val) then
+      pp(ilev)=1.e33
+    else
+      pp(ilev)=press(ilon,ilat,ilev,itim)
+    endif
+  enddo
   do ilev=altlength,2,-1
     if ((pp(ilev).le.tmpp).and.(pp(ilev-1).gt.tmpp)) then
@@ -163,5 +176,5 @@
    .and.  (grav(ilon,ilat,ilev,itim).ne.miss_val)) then
      dmass(ilon,ilat,ilev,itim) = &
-               rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))*deltalon &
+               rayon(ilon,ilat,ilev,itim)*coslat(ilat)*deltalon &
              * rayon(ilon,ilat,ilev,itim)*deltalat &
              * deltap(ilon,ilat,ilev,itim)/grav(ilon,ilat,ilev,itim)
Index: /trunk/LMDZ.VENUS/Tools/dx_dp.F
===================================================================
--- /trunk/LMDZ.VENUS/Tools/dx_dp.F	(revision 1453)
+++ /trunk/LMDZ.VENUS/Tools/dx_dp.F	(revision 1454)
@@ -30,5 +30,5 @@
 c-----------------------------------------------------------------------
       do j=1,jjp1
-         if ((x(j,1).lt.indefini).and.(x(j,2).lt.indefini)) then  
+         if ((x(j,1).ne.indefini).and.(x(j,2).ne.indefini)) then  
            dxdp(j,1) = (x(j,2)-x(j,1))/(pniv(2) - pniv(1))
          else
@@ -37,9 +37,9 @@
         
          do l=2,llm-1
-           if ((x(j,l-1).lt.indefini).and.(x(j,l+1).lt.indefini))then
+           if ((x(j,l-1).ne.indefini).and.(x(j,l+1).ne.indefini))then
              dxdp(j,l)= (x(j,l+1)-x(j,l-1))/(pniv(l+1) - pniv(l-1))
-         else if((x(j,l+1).lt.indefini).and.(x(j,l).lt.indefini))then
+         else if((x(j,l+1).ne.indefini).and.(x(j,l).ne.indefini))then
              dxdp(j,l)= (x(j,l+1)-x(j,l))  /(pniv(l+1) - pniv(l))
-         else if((x(j,l-1).lt.indefini).and.(x(j,l).lt.indefini))then
+         else if((x(j,l-1).ne.indefini).and.(x(j,l).ne.indefini))then
              dxdp(j,l)= (x(j,l)-x(j,l-1))  /(pniv(l)   - pniv(l-1))
            else 
@@ -47,5 +47,5 @@
            end if
          end do
-         if ((x(j,llm).lt.indefini).and.(x(j,llm-1).lt.indefini)) then
+         if ((x(j,llm).ne.indefini).and.(x(j,llm-1).ne.indefini)) then
          dxdp(j,llm)= (x(j,llm)-x(j,llm-1))/(pniv(llm)-pniv(llm-1))
          else
@@ -58,5 +58,5 @@
       do j=1,jjp1
          do l = 1, llm
-           if ( (abs(dxdp(j,l)).gt.(indefini/100.)).and.
+           if ( (abs(dxdp(j,l)).gt.(abs(indefini)/100.)).and.
      .          (dxdp(j,l).ne.indefini)) then
               write(*,*) '----> j= ', j , ' l= ' ,  l
Index: /trunk/LMDZ.VENUS/Tools/energy.F90
===================================================================
--- /trunk/LMDZ.VENUS/Tools/energy.F90	(revision 1453)
+++ /trunk/LMDZ.VENUS/Tools/energy.F90	(revision 1454)
@@ -36,11 +36,11 @@
 integer,dimension(4) :: datashape4d ! shape of 4D datasets
 
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
+real :: miss_val ! special "missing value" to specify missing data
+real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
 real :: pi
 real,dimension(:),allocatable :: lon ! longitude
 integer lonlength ! # of grid points along longitude
 real,dimension(:),allocatable :: lat ! latitude
-real,dimension(:),allocatable :: latrad ! latitude in radian
+real,dimension(:),allocatable :: coslat ! cos of latitude
 integer latlength ! # of grid points along latitude
 real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
@@ -80,4 +80,5 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
@@ -123,9 +124,10 @@
 if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
 
-allocate(latrad(latlength))
-latrad = lat*pi/180.
+allocate(coslat(latlength))
+! Beware of rounding problems at poles...
+coslat(:) = max(0.,cos(lat(:)*pi/180.))
 
 ! Lat, lon pressure intervals
-deltalat = abs(latrad(2)-latrad(1))
+deltalat = abs(lat(2)-lat(1))*pi/180.
 deltalon = abs(lon(2)-lon(1))*pi/180.
 
@@ -193,9 +195,9 @@
   text="T"
 endif
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "  looking for t instead... "
   text="t"
-  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
   if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
 endif
@@ -216,5 +218,5 @@
   text="U"
 endif
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
@@ -228,5 +230,5 @@
   text="V"
 endif
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
@@ -242,5 +244,5 @@
 
 !text="zareoid"
-!call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
+!call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
 !if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
 !if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
@@ -276,6 +278,6 @@
 enddo ! timelength
 
-call cellmass(infid,latlength,lonlength,altlength,timelength, &
-              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
+call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
+              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
               dmass )
 
Index: /trunk/LMDZ.VENUS/Tools/epflux.F90
===================================================================
--- /trunk/LMDZ.VENUS/Tools/epflux.F90	(revision 1453)
+++ /trunk/LMDZ.VENUS/Tools/epflux.F90	(revision 1454)
@@ -1,3 +1,3 @@
-subroutine  epflux(iip1,jjp1,llm,indefini,rlatu,rbar &
+subroutine  epflux(iip1,jjp1,llm,indefini,latdeg,rbar &
                   ,teta,u3d,v3d,w3d,press &
                   ,epy2d,epz2d,divep2d,vtem2d,wtem2d,acc_meridien2d &
@@ -41,5 +41,5 @@
       
       integer :: iip1,jjp1,llm
-      real :: indefini,rlatu(jjp1),rbar(jjp1,llm)
+      real :: indefini,latdeg(jjp1),rbar(jjp1,llm)
       REAL :: teta(iip1,jjp1,llm)
       REAL :: u3d(iip1,jjp1,llm),v3d(iip1,jjp1,llm)
@@ -75,4 +75,5 @@
 !  Autres
 !  ------
+      real :: rlatu(jjp1) ! lat in radians (beware of rounding effects at poles)
       real :: pi
       REAL :: tetas(llm)
@@ -88,4 +89,7 @@
       jjm = jjp1-1
       pi = 2.*asin(1.)
+
+! To avoid rounding effects at poles
+      rlatu(:)=latdeg(:)*pi/180.00001
 
 !     Calcul de moyennes zonales
@@ -108,5 +112,5 @@
          n= 0 
          do j=2,jjm
-            if (tetabar(j,l).lt.indefini) then 
+            if (tetabar(j,l).ne.indefini) then 
                tetas(l) = tetas(l) + tetabar(j,l)
                n = n+1 
@@ -237,5 +241,5 @@
       do l=1,llm
 
-       if ((ubar(2,l).lt.indefini).and.(ubar(1,l).lt.indefini))then
+       if ((ubar(2,l).ne.indefini).and.(ubar(1,l).ne.indefini))then
       ducosdlat(1,l) = ( ubar(2,l)*cos(rlatu(2))-ubar(1,l)*cos(rlatu(1)) ) & 
                       /( rlatu(2) - rlatu(1)) 
@@ -244,5 +248,5 @@
        end if
         do j=2,jjm
-       if ((ubar(j+1,l).lt.indefini).and.(ubar(j-1,l).lt.indefini))then
+       if ((ubar(j+1,l).ne.indefini).and.(ubar(j-1,l).ne.indefini))then
       ducosdlat(j,l) = ( ubar(j+1,l)*cos(rlatu(j+1))-ubar(j-1,l)*cos(rlatu(j-1))) &
                       /( rlatu(j+1) - rlatu(j-1))
@@ -251,5 +255,5 @@
        end if
         enddo
-       if ((ubar(jjp1,l).lt.indefini).and.(ubar(jjm,l).lt.indefini))then
+       if ((ubar(jjp1,l).ne.indefini).and.(ubar(jjm,l).ne.indefini))then
       ducosdlat(jjp1,l) = ( ubar(jjp1,l)*cos(rlatu(jjp1))  &
                                -ubar(jjm,l)*cos(rlatu(jjm))  ) &
@@ -260,6 +264,6 @@
            
         do j=1,jjp1
-       if ((vptetapbar(j,l).lt.indefini).and. &
-           (dtetabardp(j,l).lt.indefini)) then
+       if ((vptetapbar(j,l).ne.indefini).and. &
+           (dtetabardp(j,l).ne.indefini)) then
            vz(j,l) = vptetapbar(j,l)/dtetabardp(j,l)
          wlat(j,l) = cos(rlatu(j))*vz(j,l)
@@ -284,8 +288,8 @@
 !     ------------
 
-            if (   (rbar(j,l).lt.indefini).and. &
-                (vpupbar(j,l).lt.indefini).and. &
-                (dubardp(j,l).lt.indefini).and. &
-                     (vz(j,l).lt.indefini)) then
+            if (   (rbar(j,l).ne.indefini).and. &
+                (vpupbar(j,l).ne.indefini).and. &
+                (dubardp(j,l).ne.indefini).and. &
+                     (vz(j,l).ne.indefini)) then
            
      epy2d(j,l) = rbar(j,l) * cos(rlatu(j))* ( vpupbar(j,l) &
@@ -303,8 +307,8 @@
 !     ------------
 
-            if (   (rbar(j,l).lt.indefini).and. &
-                (wpupbar(j,l).lt.indefini).and. &
-              (ducosdlat(j,l).lt.indefini).and. &
-                     (vz(j,l).lt.indefini)) then
+            if (   (rbar(j,l).ne.indefini).and. &
+                (wpupbar(j,l).ne.indefini).and. &
+              (ducosdlat(j,l).ne.indefini).and. &
+                     (vz(j,l).ne.indefini)) then
 
      epz2d(j,l) = rbar(j,l) * cos(rlatu(j))*                    &
@@ -330,6 +334,6 @@
       do l=1,llm
 
-          if ( (wlat(2,l).lt.indefini) .and. &
-               (wlat(1,l).lt.indefini) ) then
+          if ( (wlat(2,l).ne.indefini) .and. &
+               (wlat(1,l).ne.indefini) ) then
                dwlatdlat(1,l)=(wlat(2,l)-wlat(1,l)) &
                      /(rlatu(2)-rlatu(1))
@@ -338,6 +342,6 @@
           end if
         do j=2,jjm
-          if ( (wlat(j+1,l).lt.indefini) .and. &
-               (wlat(j-1,l).lt.indefini) ) then
+          if ( (wlat(j+1,l).ne.indefini) .and. &
+               (wlat(j-1,l).ne.indefini) ) then
                dwlatdlat(j,l)=(wlat(j+1,l)-wlat(j-1,l)) &
                      /(rlatu(j+1)-rlatu(j-1))
@@ -346,6 +350,6 @@
           end if
         enddo 
-          if ( (wlat(jjp1,l).lt.indefini) .and. & 
-               (wlat(jjm, l).lt.indefini) ) then
+          if ( (wlat(jjp1,l).ne.indefini) .and. & 
+               (wlat(jjm, l).ne.indefini) ) then
                dwlatdlat(jjp1,l)=(wlat(jjp1,l)-wlat(jjm,l)) &
                      /(rlatu(jjp1)-rlatu(jjm))
@@ -355,6 +359,6 @@
 
         do j=1,jjp1
-          if ( (dvzdp(j,l).lt.indefini) &
-           .and.(vbar(j,l).lt.indefini)) then
+          if ( (dvzdp(j,l).ne.indefini) &
+           .and.(vbar(j,l).ne.indefini)) then
             vtem2d(j,l) = vbar(j,l)-dvzdp(j,l)
           else
@@ -364,7 +368,7 @@
 
         do j=1,jjp1
-          if (   (rbar(j,l).lt.indefini).and. &
-            (dwlatdlat(j,l).lt.indefini).and. &
-                 (wbar(j,l).lt.indefini)) then
+          if (   (rbar(j,l).ne.indefini).and. &
+            (dwlatdlat(j,l).ne.indefini).and. &
+                 (wbar(j,l).ne.indefini)) then
            wtem2d(j,l) = wbar(j,l)+dwlatdlat(j,l)/(rbar(j,l)*cos(rlatu(j)))
           else
@@ -385,9 +389,9 @@
       do l=1,llm
         do j=1,jjp1
-          if (   (rbar(j,l).lt.indefini) &
-          .and.(vtem2d(j,l).lt.indefini) &
-          .and.(wtem2d(j,l).lt.indefini) &
-       .and.(ducosdlat(j,l).lt.indefini) &
-         .and.(dubardp(j,l).lt.indefini)) then
+          if (   (rbar(j,l).ne.indefini) &
+          .and.(vtem2d(j,l).ne.indefini) &
+          .and.(wtem2d(j,l).ne.indefini) &
+       .and.(ducosdlat(j,l).ne.indefini) &
+         .and.(dubardp(j,l).ne.indefini)) then
        acc_meridien2d(j,l) =  &
            vtem2d(j,l)*(ducosdlat(j,l)/(rbar(j,l)*cos(rlatu(j))) &
@@ -410,8 +414,8 @@
         divep2d(1,l) =0
         do j=2,jjm
-          if (    (rbar(j,l).lt.indefini) .and. &
-               (epy2d(j+1,l).lt.indefini) .and. &
-               (epy2d(j-1,l).lt.indefini) .and. &
-                (depzdp(j,l).lt.indefini)    ) then
+          if (    (rbar(j,l).ne.indefini) .and. &
+               (epy2d(j+1,l).ne.indefini) .and. &
+               (epy2d(j-1,l).ne.indefini) .and. &
+                (depzdp(j,l).ne.indefini)    ) then
 
            depycosdlat=(epy2d(j+1,l)*cos(rlatu(j+1))   &
Index: /trunk/LMDZ.VENUS/Tools/fft.F90
===================================================================
--- /trunk/LMDZ.VENUS/Tools/fft.F90	(revision 1453)
+++ /trunk/LMDZ.VENUS/Tools/fft.F90	(revision 1454)
@@ -63,11 +63,10 @@
 integer,dimension(4) :: datashape4d ! shape of 4D datasets
 
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
+real :: miss_val ! special "missing value" to specify missing data
+real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
 real :: pi
 real,dimension(:),allocatable :: lon ! longitude
 integer lonlength ! # of grid points along longitude
 real,dimension(:),allocatable :: lat ! latitude
-real,dimension(:),allocatable :: latrad ! latitude in radian
 integer latlength ! # of grid points along latitude
 real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
@@ -125,6 +124,11 @@
 logical :: lmdflag
 
+! Tuning parameters
+real :: fcoup1,fcoup2,width
+real :: fcoup1tmp,fcoup2tmp,widthtmp
+logical,dimension(4) :: ok_out
+character (len=1) :: ok_outtmp
+
 include "planet.h"
-include "filter.h"
 
 #include <fftw3.f> 
@@ -135,7 +139,58 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
 write(*,*) "You are working on the atmosphere of ",planet
+
+! initialisation
+!----------------
+
+! Par defaut
+
+! Define the filters
+! Low  cutting frequency, in Hz    : fcoup1
+fcoup1=6.e-7
+! High cutting frequency, in Hz    : fcoup2
+fcoup2=1.4e-6
+! Half-width of the filters, in Hz : width
+width=2.e-7
+! Outputs (U,     V,      W,     T)
+ok_out=(/.true.,.true.,.false.,.true./)
+
+print*,"Low  cutting frequency, in Hz ?"
+print*,"between 1e-5 and 1e-7, 0 for default => 6.e-7"
+read(*,*) fcoup1tmp
+if ((fcoup1tmp.lt.1e-5).and.(fcoup1tmp.gt.1e-7)) fcoup1=fcoup1tmp
+print*,"=",fcoup1
+print*,"High cutting frequency, in Hz ?"
+print*,"between 1e-5 and 1e-7, 0 for default => 1.4e-6"
+read(*,*) fcoup2tmp
+if ((fcoup2tmp.lt.1e-5).and.(fcoup2tmp.gt.1e-7)) fcoup2=fcoup2tmp
+print*,"=",fcoup2
+print*,"Half-width of the filters, in Hz ?"
+print*,"between 1e-6 and 1e-8, 0 for default => 2e-7)"
+read(*,*) widthtmp
+if ((widthtmp.lt.1e-6).and.(widthtmp.gt.1e-8)) width=widthtmp
+print*,"=",width
+!width = 3./time(timelength)
+
+! Outputs 
+print*,"Output of zonal wind ? (y or n, default is y)"
+read(*,'(a1)') ok_outtmp
+if (ok_outtmp.eq."n") ok_out(1)=.false.
+print*,"=",ok_out(1)
+print*,"Output of meridional wind ? (y or n, default is y)"
+read(*,'(a1)') ok_outtmp
+if (ok_outtmp.eq."n") ok_out(2)=.false.
+print*,"=",ok_out(2)
+print*,"Output of vertical wind ? (y or n, default is n)"
+read(*,'(a1)') ok_outtmp
+if (ok_outtmp.eq."y") ok_out(3)=.true.
+print*,"=",ok_out(3)
+print*,"Output of temperature ? (y or n, default is y)"
+read(*,'(a1)') ok_outtmp
+if (ok_outtmp.eq."n") ok_out(4)=.false.
+print*,"=",ok_out(4)
 
 !===============================================================================
@@ -172,7 +227,4 @@
 ierr=NF_GET_VAR_REAL(infid,lat_varid,lat)
 if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
-
-allocate(latrad(latlength))
-latrad = lat*pi/180.
 
 allocate(plev(altlength))
@@ -212,9 +264,9 @@
 
 text="temp"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "  looking for t instead... "
   text="t"
-  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
   if (ierr1.ne.NF_NOERR) then
     print*,"Error: Failed to get temperature ID"
@@ -236,5 +288,5 @@
 
 text="vitu"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   print*,"Error: Failed to get vitu ID"
@@ -252,5 +304,5 @@
 
 text="vitv"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   print*,"Error: Failed to get vitv ID"
@@ -268,5 +320,5 @@
 
 text="vitw"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   print*,"Error: Failed to get vitw ID"
@@ -288,4 +340,5 @@
 ! 2.2.1 FFT and filtering
 !===============================================================================
+
 ! allocations
 !-------------
@@ -405,19 +458,4 @@
 enddo
 
-! Define the filters
-
-! IN filter.h :
-! Low  cutting frequency, in Hz    : fcoup1
-! High cutting frequency, in Hz    : fcoup2
-! Half-width of the filters, in Hz : width
-
-!print*,"Low  cutting frequency, in Hz ?"
-!read(*,*) fcoup1
-!print*,"High cutting frequency, in Hz ?"
-!read(*,*) fcoup2
-!print*,"Half-width of the filters, in Hz ?"
-!read(*,*) width
-!width = 3./time(timelength)
-
 do itim=1,M_fft+1
   if (freq(itim).lt.(fcoup1-width)) then
Index: /trunk/LMDZ.VENUS/Tools/io.F90
===================================================================
--- /trunk/LMDZ.VENUS/Tools/io.F90	(revision 1453)
+++ /trunk/LMDZ.VENUS/Tools/io.F90	(revision 1454)
@@ -218,5 +218,5 @@
 !===========================================================================
 
-subroutine get_var4d(infid,dim1,dim2,dim3,dim4,text,var,ierr1,ierr2)
+subroutine get_var4d(infid,dim1,dim2,dim3,dim4,text,var,missing,ierr1,ierr2)
 
 implicit none
@@ -229,5 +229,6 @@
 character (len=64) :: text ! name of variable to read
 real,dimension(dim1,dim2,dim3,dim4) :: var ! variable to read
-integer :: ierr1,ierr2 ! NetCDF routines return code
+real :: missing ! missing value
+integer :: ierr1,ierr2,miss ! NetCDF routines return code
 
 ! local
@@ -240,4 +241,5 @@
   write(*,*) "ID ok for ",trim(text)
   ierr2=NF_GET_VAR_REAL(infid,tmpvarid,var)
+  miss=NF_GET_ATT_REAL(infid,tmpvarid,"missing_value",missing)
 endif
 
Index: /trunk/LMDZ.VENUS/Tools/moytim.F
===================================================================
--- /trunk/LMDZ.VENUS/Tools/moytim.F	(revision 1453)
+++ /trunk/LMDZ.VENUS/Tools/moytim.F	(revision 1454)
@@ -33,5 +33,5 @@
            n = 0
            do t=1,nt
-              if (x(i,j,l,t).lt.indefini) then 
+              if (x(i,j,l,t).ne.indefini) then 
                  xmean(i,j,l) = xmean(i,j,l) + x(i,j,l,t)
                  n = n+1
Index: /trunk/LMDZ.VENUS/Tools/moyzon.F
===================================================================
--- /trunk/LMDZ.VENUS/Tools/moyzon.F	(revision 1453)
+++ /trunk/LMDZ.VENUS/Tools/moyzon.F	(revision 1454)
@@ -32,5 +32,5 @@
            n = 0
            do i=1,iim
-              if (x(i,j,l).lt.indefini) then 
+              if (x(i,j,l).ne.indefini) then 
                  xbar(j,l) = xbar(j,l) + x(i,j,l)
                  n = n+1
Index: /trunk/LMDZ.VENUS/Tools/moyzon2.F
===================================================================
--- /trunk/LMDZ.VENUS/Tools/moyzon2.F	(revision 1453)
+++ /trunk/LMDZ.VENUS/Tools/moyzon2.F	(revision 1454)
@@ -34,5 +34,5 @@
            n = 0
            do i=1,iim
-              if (x(i,j,l).lt.indefini) then 
+              if (x(i,j,l).ne.indefini) then 
                  xbar(j,l) = xbar(j,l) + x(i,j,l)
                  n = n+1
Index: /trunk/LMDZ.VENUS/Tools/psi.F90
===================================================================
--- /trunk/LMDZ.VENUS/Tools/psi.F90	(revision 1453)
+++ /trunk/LMDZ.VENUS/Tools/psi.F90	(revision 1454)
@@ -31,11 +31,11 @@
 integer,dimension(4) :: datashape4d ! shape of 4D datasets
 
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
+real :: miss_val ! special "missing value" to specify missing data
+real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
 real :: pi
 real,dimension(:),allocatable :: lon ! longitude
 integer lonlength ! # of grid points along longitude
 real,dimension(:),allocatable :: lat ! latitude
-real,dimension(:),allocatable :: latrad ! latitude in radian
+real,dimension(:),allocatable :: coslat ! cos of latitude
 integer latlength ! # of grid points along latitude
 real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
@@ -66,4 +66,5 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
@@ -104,9 +105,10 @@
 if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
 
-allocate(latrad(latlength))
-latrad = lat*pi/180.
+allocate(coslat(latlength))
+! Beware of rounding problems at poles...
+coslat(:) = max(0.,cos(lat(:)*pi/180.))
 
 ! Lat, lon pressure intervals
-deltalat = abs(latrad(2)-latrad(1))
+deltalat = abs(lat(2)-lat(1))*pi/180.
 deltalon = abs(lon(2)-lon(1))*pi/180.
 
@@ -160,5 +162,5 @@
 ! meridional wind vitv (in m/s)
 text="vitv"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
@@ -173,9 +175,9 @@
 
 text="zareoid"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "  looking for geop instead... "
   text="geop"
-  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
+  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
   if (ierr1.ne.NF_NOERR) stop "Error: Failed to get geop ID"
 do itim=1,timelength
@@ -226,5 +228,5 @@
 
 call cellmass(infid,latlength,lonlength,altlength,timelength, &
-              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
+              lmdflag,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
               dmass )
 
Index: /trunk/LMDZ.VENUS/Tools/stability.F90
===================================================================
--- /trunk/LMDZ.VENUS/Tools/stability.F90	(revision 1453)
+++ /trunk/LMDZ.VENUS/Tools/stability.F90	(revision 1454)
@@ -31,11 +31,11 @@
 integer,dimension(4) :: datashape4d ! shape of 4D datasets
 
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
+real :: miss_val ! special "missing value" to specify missing data
+real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
 real :: pi
 real,dimension(:),allocatable :: lon ! longitude
 integer lonlength ! # of grid points along longitude
 real,dimension(:),allocatable :: lat ! latitude
-real,dimension(:),allocatable :: latrad ! latitude in radian
+real,dimension(:),allocatable :: coslat ! cos of latitude
 integer latlength ! # of grid points along latitude
 real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
@@ -91,4 +91,5 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
@@ -129,9 +130,10 @@
 if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
 
-allocate(latrad(latlength))
-latrad = lat*pi/180.
+allocate(coslat(latlength))
+! Beware of rounding problems at poles...
+coslat(:) = max(0.,cos(lat(:)*pi/180.))
 
 ! Lat, lon pressure intervals
-deltalat = abs(latrad(2)-latrad(1))
+deltalat = abs(lat(2)-lat(1))*pi/180.
 deltalon = abs(lon(2)-lon(1))*pi/180.
 
@@ -180,9 +182,9 @@
 
 text="temp"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "  looking for t instead... "
   text="t"
-  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
   if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
 endif
@@ -197,5 +199,5 @@
 ! zonal wind vitu (in m/s)
 text="vitu"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
@@ -203,5 +205,5 @@
 ! meridional wind vitv (in m/s)
 text="vitv"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
@@ -215,5 +217,5 @@
 
 ! text="zareoid"
-! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
+! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
 ! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
 ! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
@@ -276,6 +278,6 @@
  grav(lonlength+1,:,:,:) =  grav(1,:,:,:)
 
-call cellmass(infid,latlength,lonlength+1,altlength,timelength, &
-              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
+call cellmass(infid,latlength,lonlength+1,altlength,timelength,lmdflag, &
+              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
               dmass )
 
@@ -364,6 +366,6 @@
 
 do ilat=2,latlength-1
-   if (tan(latrad(ilat)).ne.0.) then
-     fac1 = R0/tan(latrad(ilat))
+   if (tan(lat(ilat)*pi/180.).ne.0.) then
+     fac1 = R0/tan(lat(ilat)*pi/180.)
    else
      fac1 = miss_val
Index: /trunk/LMDZ.VENUS/Tools/tem.F90
===================================================================
--- /trunk/LMDZ.VENUS/Tools/tem.F90	(revision 1453)
+++ /trunk/LMDZ.VENUS/Tools/tem.F90	(revision 1454)
@@ -37,11 +37,11 @@
 integer,dimension(3) :: datashape3d ! shape of 3D datasets
 
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
+real :: miss_val ! special "missing value" to specify missing data
+real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
 real :: pi
 real,dimension(:),allocatable :: lon ! longitude
 integer lonlength ! # of grid points along longitude
 real,dimension(:),allocatable :: lat ! latitude
-real,dimension(:),allocatable :: latrad ! latitude in radian
+real,dimension(:),allocatable :: coslat ! cos of latitude
 integer latlength ! # of grid points along latitude
 real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
@@ -107,4 +107,5 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
@@ -145,9 +146,10 @@
 if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
 
-allocate(latrad(latlength))
-latrad = lat*pi/180.
+allocate(coslat(latlength))
+! Beware of rounding problems at poles...
+coslat(:) = max(0.,cos(lat(:)*pi/180.))
 
 ! Lat, lon pressure intervals
-deltalat = abs(latrad(2)-latrad(1))
+deltalat = abs(lat(2)-lat(1))*pi/180.
 deltalon = abs(lon(2)-lon(1))*pi/180.
 
@@ -196,9 +198,9 @@
 
 text="temp"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) then
   write(*,*) "  looking for t instead... "
   text="t"
-  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,ierr1,ierr2)
+  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,temp,miss_val,ierr1,ierr2)
   if (ierr1.ne.NF_NOERR) stop "Error: Failed to get temperature ID"
 endif
@@ -214,5 +216,5 @@
 ! zonal wind vitu (in m/s)
 text="vitu"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
@@ -220,5 +222,5 @@
 ! meridional wind vitv (in m/s)
 text="vitv"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
@@ -226,5 +228,5 @@
 ! vertical wind vitw (in Pa/s)
 text="vitw"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitw ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading vertical wind"
@@ -238,5 +240,5 @@
 
 ! text="zareoid"
-! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
+! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
 ! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
 ! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
@@ -312,6 +314,6 @@
  grav(lonlength+1,:,:,:) =  grav(1,:,:,:)
 
-call cellmass(infid,latlength,lonlength+1,altlength,timelength, &
-              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
+call cellmass(infid,latlength,lonlength+1,altlength,timelength,lmdflag, &
+              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
               dmass )
 
@@ -359,5 +361,5 @@
 call moyzon(lonlength,latlength,altlength,miss_val,u3d,ubar)
 
-call epflux(lonlength+1,latlength,altlength,miss_val,latrad,rbar &
+call epflux(lonlength+1,latlength,altlength,miss_val,lat,rbar &
             ,teta,u3d,v3d,w3d,plev &
             ,epy2d,epz2d,div2d,vtem2d,wtem2d,ammc2d &
@@ -376,5 +378,5 @@
 ! rsurg: r*dp/g = dm/(r cos(lat) dlon dlat) !!!
      rsurg(ilon,ilat,ilev) = dmass(ilon,ilat,ilev,itim) &
-          / (r3d(ilon,ilat,ilev)*cos(latrad(ilat))*deltalon*deltalat)
+          / (r3d(ilon,ilat,ilev)*coslat(ilat)*deltalon*deltalat)
     else
      rsurg(ilon,ilat,ilev) = miss_val
@@ -391,5 +393,5 @@
         (rsurgbar(ilat,ilev).ne.miss_val) ) then
       vm(ilat,ilev) = vtem2d(ilat,ilev) &
-            * 2.*pi*rsurgbar(ilat,ilev)*cos(latrad(ilat))
+            * 2.*pi*rsurgbar(ilat,ilev)*coslat(ilat)
     else
       vm(ilat,ilev) = miss_val
Index: /trunk/LMDZ.VENUS/Tools/tmc.F90
===================================================================
--- /trunk/LMDZ.VENUS/Tools/tmc.F90	(revision 1453)
+++ /trunk/LMDZ.VENUS/Tools/tmc.F90	(revision 1454)
@@ -49,11 +49,11 @@
 integer,dimension(4) :: datashape4d ! shape of 4D datasets
 
-real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
-real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
+real :: miss_val ! special "missing value" to specify missing data
+real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
 real :: pi
 real,dimension(:),allocatable :: lon ! longitude
 integer lonlength ! # of grid points along longitude
 real,dimension(:),allocatable :: lat ! latitude
-real,dimension(:),allocatable :: latrad ! latitude in radian
+real,dimension(:),allocatable :: coslat ! cos of latitude
 integer latlength ! # of grid points along latitude
 real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
@@ -141,4 +141,5 @@
 
 pi = 2.*asin(1.)
+miss_val = miss_val_def
 
 write(*,*) ""
@@ -179,9 +180,10 @@
 if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
 
-allocate(latrad(latlength))
-latrad = lat*pi/180.
+allocate(coslat(latlength))
+! Beware of rounding problems at poles...
+coslat(:) = max(0.,cos(lat(:)*pi/180.))
 
 ! Lat, lon pressure intervals
-deltalat = abs(latrad(2)-latrad(1))
+deltalat = abs(lat(2)-lat(1))*pi/180.
 deltalon = abs(lon(2)-lon(1))*pi/180.
 
@@ -233,5 +235,5 @@
 ! zonal wind vitu (in m/s)
 text="vitu"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitu ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
@@ -239,5 +241,5 @@
 ! meridional wind vitv (in m/s)
 text="vitv"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitv,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitv ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading meridional wind"
@@ -245,5 +247,5 @@
 ! vertical wind vitw (in Pa/s)
 text="vitw"
-call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,ierr1,ierr2)
+call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitw,miss_val,ierr1,ierr2)
 if (ierr1.ne.NF_NOERR) stop "Error: Failed to get vitw ID"
 if (ierr2.ne.NF_NOERR) stop "Error: Failed reading vertical wind"
@@ -257,5 +259,5 @@
 
 ! text="zareoid"
-! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,ierr1,ierr2)
+! call get_var4d(infid,lonlength,latlength,altlength,timelength,text,za,miss_val,ierr1,ierr2)
 ! if (ierr1.ne.NF_NOERR) stop "Error: Failed to get za ID"
 ! if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zareoid"
@@ -293,6 +295,6 @@
 enddo ! timelength
 
-call cellmass(infid,latlength,lonlength,altlength,timelength, &
-              lmdflag,deltalon,deltalat,latrad,plev,ps,grav,rayon, &
+call cellmass(infid,latlength,lonlength,altlength,timelength,lmdflag, &
+              miss_val,deltalon,deltalat,coslat,plev,ps,grav,rayon, &
               dmass )
 
@@ -316,8 +318,8 @@
       osam(ilon,ilat,ilev,itim) = &
          rayon(ilon,ilat,ilev,itim)*rayon(ilon,ilat,ilev,itim) &
-       * cos(latrad(ilat))*cos(latrad(ilat)) &
+       * coslat(ilat)*coslat(ilat) &
        * omega
       rsam(ilon,ilat,ilev,itim) = vitu(ilon,ilat,ilev,itim) &
-       * rayon(ilon,ilat,ilev,itim)*cos(latrad(ilat))
+       * rayon(ilon,ilat,ilev,itim)*coslat(ilat)
       tsam(ilon,ilat,ilev,itim) = osam(ilon,ilat,ilev,itim)&
                                 + rsam(ilon,ilat,ilev,itim)
