Index: trunk/LMDZ.PLUTO/util/zrecast_plut.F90
===================================================================
--- trunk/LMDZ.PLUTO/util/zrecast_plut.F90	(revision 3376)
+++ trunk/LMDZ.PLUTO/util/zrecast_plut.F90	(revision 3438)
@@ -40,5 +40,5 @@
 !              'phisinit.nc')
 ! EM 02/2007 : Changed behavior for "altitude above surface" case
-!              (for MCD RMS computation). Number of levels is then set as 
+!              (for MCD RMS computation). Number of levels is then set as
 !              number of levels in initial file,
 !              and the new set of above surface levels follow a more elaborate
@@ -99,5 +99,5 @@
 real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
 real,dimension(:,:,:),allocatable :: ps ! GCM surface pressure
-real,dimension(:,:,:,:),allocatable :: press ! GCM atmospheric pressure 
+real,dimension(:,:,:,:),allocatable :: press ! GCM atmospheric pressure
 real,dimension(:,:,:,:),allocatable :: temp ! GCM atmospheric temperature
 real,dimension(:,:,:,:),allocatable :: rho ! GCM atmospheric density
@@ -195,5 +195,5 @@
     if (tmpndims.eq.4) then
       nbvar=nbvar+1
-      var(nbvar)=tmpvarname    
+      var(nbvar)=tmpvarname
     else
       write(*,*) 'Error: ',trim(tmpvarname),' is not a 4D variable'
@@ -490,5 +490,5 @@
 endif
 ! look for 'phisinit' in current file
-ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid) 
+ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
 if (ierr.ne.NF_NOERR) then
   write(*,*) "Warning: Failed to get phisinit ID from file ",trim(infile)
@@ -598,5 +598,5 @@
     if (nblev.eq.1) then ! in case only one level is asked for
       zareoid(nblev)=zamin
-    else 
+    else
       do i=1,nblev
         zareoid(i)=zamin+(real(i)-1.0)*((zamax-zamin)/(real(nblev)-1.0))
@@ -625,5 +625,5 @@
     if (nblev.eq.1) then ! in case only one level is asked for
       zradius(nblev)=zamin
-    else 
+    else
       do i=1,nblev
         zradius(i)=zamin+(real(i)-1.0)*((zamax-zamin)/(real(nblev)-1.0))
@@ -977,5 +977,5 @@
 ! 3.3.1 Define 1D variables
 
-! longitude 
+! longitude
 datashape(1)=lon_dimid
 ierr=NF_DEF_VAR(outfid,"longitude",NF_REAL,1,datashape(1),lon_varid)
@@ -1156,5 +1156,6 @@
   stop "Error: Problem writing long_name for Time"
 endif
-text='days since 0000-01-1 00:00:00'
+!text='days since 0000-01-1 00:00:00'
+text='days since 00000'
 ierr=NF_PUT_ATT_TEXT(outfid,time_varid,'units',len_trim(text),text)
 if (ierr.ne.NF_NOERR) then
@@ -1259,5 +1260,5 @@
     write(*,*) 'Error, could not define variable ',trim(var(i))
     write(*,*) NF_STRERROR(ierr)
-    stop 
+    stop
   endif
 
@@ -1320,5 +1321,5 @@
     endif
   endif
-  
+
   ! look for a "units" attribute
   text='   '
@@ -1334,5 +1335,5 @@
     endif
   endif
-  
+
   ! look for a "missing_value" attribute
   ierr=NF_GET_ATT_REAL(infid,tmpvarid,"missing_value",miss_val)
@@ -1343,5 +1344,5 @@
     miss_val=miss_val_def
   endif
-  
+
   ! write the missing_value attribute to output
   ierr=NF_PUT_ATT_REAL(outfid,var_id(i),'missing_value',NF_REAL,1,miss_val)
@@ -1356,5 +1357,5 @@
     write(*,*)'nbattr:',nbattr,' j:',j
   endif
-    
+
 enddo ! of do=1,nbvar
 
@@ -1499,9 +1500,9 @@
     stop
   endif
-  
+
   ! interpolate dataset onto new grid
   call p_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
                       miss_val,ps,press,indata,plevel,outdata)
-  
+
   ! write new dataset to output file
   ierr=NF_PUT_VAR_REAL(outfid,var_id(i),outdata)
@@ -1511,7 +1512,7 @@
     stop
   endif
-  
+
   enddo ! of do i=1,nbvar
-  
+
   ! interpolate zareoid onto new grid
   call p_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
@@ -1540,5 +1541,5 @@
     stop
   endif
-  
+
   ! interpolate dataset onto new grid
   ! check if variable is "rho" (to set flag for interpolation below)
@@ -1548,5 +1549,5 @@
     j=0
   endif
-  
+
   call z_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
                       miss_val,za_gcm,indata,j,zareoid,outdata)
@@ -1561,5 +1562,5 @@
 
   enddo ! of do i=1,nbvar
-  
+
   ! interpolate pressure onto new grid
   write(*,*) ""
@@ -1593,5 +1594,5 @@
       stop
     endif
-  
+
     ! interpolate dataset onto new grid
     ! check if variable is "rho" (to set flag for interpolation below)
@@ -1601,5 +1602,5 @@
       j=0
     endif
-  
+
     call zs_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
                         miss_val,zs_gcm,indata,phisinit,press,temp,rho,&
@@ -1614,5 +1615,5 @@
     endif
   enddo ! of do i=1,nbvar
-  
+
   ! interpolate pressure onto new grid
   write(*,*) ""
@@ -1646,5 +1647,5 @@
     stop
   endif
-  
+
   ! interpolate dataset onto new grid
   ! check if variable is "rho" (to set flag for interpolation below)
@@ -1654,5 +1655,5 @@
     j=0
   endif
-  
+
   call z_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
                       miss_val,ra_gcm,indata,j,zradius,outdata)
@@ -1667,5 +1668,5 @@
 
   enddo ! of do i=1,nbvar
-  
+
   ! interpolate pressure onto new grid
   write(*,*) ""
@@ -1734,9 +1735,9 @@
 
 ! Parameters needed to integrate hydrostatic equation:
-real,parameter :: g0=0.62
+real,parameter :: g0=0.6198
 !g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
-real,parameter :: a0=1187.E3
+real,parameter :: a0=1188.E3
 !a0: 'mean' Mars radius=3396.km
-real :: gz 
+real :: gz
 ! gz: gravitational acceleration at a given (zareoid) altitude
 
@@ -1749,5 +1750,5 @@
 
 !==============================================================================
-! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point 
+! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point
 !   --later used to integrate the hydrostatic equation--
 !==============================================================================
@@ -1779,5 +1780,5 @@
         ! compute sigma level of layer
         sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
-        
+
         ! compute "mean" temperature of the layer
         if(temp(iloop,jloop,kloop,tloop).eq.  &
@@ -1790,17 +1791,17 @@
                      temp(iloop,jloop,kloop-1,tloop))
         endif
-        
+
         ! compute mean value of R of the layer
         Rmean=0.5*(R(iloop,jloop,kloop,tloop)+R(iloop,jloop,kloop-1,tloop))
-        
+
         ! compute gravitational acceleration (at altitude zaeroid(kloop-1))
         ! NB: zareoid=zsurface+phis/g0
         gz=g0*(a0**2)/ &
            (a0+zs_gcm(iloop,jloop,kloop-1,tloop)+phis(iloop,jloop)/g0)**2
-        
+
         ! compute zs_gcm(iloop,jloop,kloop,tloop)
         zs_gcm(iloop,jloop,kloop,tloop)=zs_gcm(iloop,jloop,kloop-1,tloop)- &
                 log(sigma(kloop)/sigma(kloop-1))*Rmean*Tmean/gz
-        
+
         ! compute za_gcm(kloop)
 !        za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
@@ -1860,9 +1861,9 @@
 
 ! Parameters needed to integrate hydrostatic equation:
-real,parameter :: g0=0.62
+real,parameter :: g0=0.6198
 !g0: exact mean gravity
-real,parameter :: a0=1187.E3
+real,parameter :: a0=1188.E3
 !a0: 'mean' radius
-real :: gz 
+real :: gz
 ! gz: gravitational acceleration at a given (zareoid) altitude
 
@@ -1875,5 +1876,5 @@
 
 !==============================================================================
-! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point 
+! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point
 !   --later used to integrate the hydrostatic equation--
 !==============================================================================
@@ -1905,5 +1906,5 @@
         ! compute sigma level of layer
         sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
-        
+
         ! compute "mean" temperature of the layer
         if(temp(iloop,jloop,kloop,tloop).eq.  &
@@ -1916,15 +1917,15 @@
                      temp(iloop,jloop,kloop-1,tloop))
         endif
-        
+
         ! compute mean value of R of the layer
         Rmean=0.5*(R(iloop,jloop,kloop,tloop)+R(iloop,jloop,kloop-1,tloop))
-        
+
         ! compute gravitational acceleration (at altitude zaeroid(kloop-1))
         gz=g0*(a0**2)/(a0+za_gcm(iloop,jloop,kloop-1,tloop))**2
-        
+
         ! compute zlocal(kloop)
         zlocal(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop-1,tloop)- &
                 log(sigma(kloop)/sigma(kloop-1))*Rmean*Tmean/gz
-        
+
         ! compute za_gcm(kloop)
         za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
@@ -1979,5 +1980,5 @@
 
 real,parameter :: P_ref=1 ! reference surface pressure used to build zsurface
-integer i 
+integer i
 
 ! 1. Build scale heights
@@ -2041,9 +2042,9 @@
 
 ! Parameters needed to integrate hydrostatic equation:
-real,parameter :: g0=0.62
+real,parameter :: g0=0.6198
 !g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
-real,parameter :: a0=1187.E3
+real,parameter :: a0=1188.E3
 !a0: 'mean' Mars radius=3396.km
-real :: gz 
+real :: gz
 ! gz: gravitational acceleration at a given (zareoid) altitude
 
@@ -2069,5 +2070,5 @@
         ! compute sigma level of layer
         sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
-        
+
         ! compute "mean" temperature of the layer
         if(temp(iloop,jloop,kloop,tloop).eq.  &
@@ -2080,14 +2081,14 @@
                      temp(iloop,jloop,kloop-1,tloop))
         endif
-                
+
         ! compute gravitational acceleration (at altitude zaeroid(kloop-1))
         ! NB: zareoid=zsurface+phis/g0
         gz=g0*(a0**2)/ &
            (a0+zs_gcm(iloop,jloop,kloop-1,tloop)+phis(iloop,jloop)/g0)**2
-        
+
         ! compute zs_gcm(iloop,jloop,kloop,tloop)
         zs_gcm(iloop,jloop,kloop,tloop)=zs_gcm(iloop,jloop,kloop-1,tloop)- &
                 log(sigma(kloop)/sigma(kloop-1))*R*Tmean/gz
-        
+
         ! compute za_gcm(kloop)
 !        za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
@@ -2143,9 +2144,9 @@
 
 ! Parameters needed to integrate hydrostatic equation:
-real,parameter :: g0=0.62
+real,parameter :: g0=0.6198
 !g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
-real,parameter :: a0=1187.E3
+real,parameter :: a0=1188.E3
 !a0: 'mean' Mars radius=3396.km
-real :: gz 
+real :: gz
 ! gz: gravitational acceleration at a given (zareoid) altitude
 
@@ -2171,5 +2172,5 @@
         ! compute sigma level of layer
         sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
-        
+
         ! compute "mean" temperature of the layer
         if(temp(iloop,jloop,kloop,tloop).eq.  &
@@ -2182,12 +2183,12 @@
                      temp(iloop,jloop,kloop-1,tloop))
         endif
-                
+
         ! compute gravitational acceleration (at altitude zaeroid(kloop-1))
         gz=g0*(a0**2)/(a0+za_gcm(iloop,jloop,kloop-1,tloop))**2
-        
+
         ! compute zlocal(kloop)
         zlocal(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop-1,tloop)- &
                 log(sigma(kloop)/sigma(kloop-1))*R*Tmean/gz
-        
+
         ! compute za_gcm(kloop)
         za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
@@ -2254,6 +2255,6 @@
         q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
       enddo
-      
-      ! Interpolation  
+
+      ! Interpolation
       do kloop=1,newaltlen
         ! compute, by interpolation, value at pressure level plevels(kloop)
@@ -2268,5 +2269,5 @@
         endif
       enddo
-      
+
     enddo !iloop
   enddo !jloop
@@ -2339,5 +2340,5 @@
         q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
       enddo !kloop
-      
+
       ! Interpolation
       do kloop=1,newaltlen
@@ -2351,9 +2352,9 @@
       else ! z_new(kloop) is out of range
         newdata(iloop,jloop,kloop,tloop)=missing
-      endif 
+      endif
       enddo !kloop
     enddo !iloop
   enddo !jloop
- enddo !tloop 
+ enddo !tloop
 else ! case when flag=1 (i.e.: rho)
  do tloop=1,tlen
@@ -2361,10 +2362,10 @@
     do iloop=1,lonlen
       do kloop=1,altlen
-        ! extract the vertical coordinates 
+        ! extract the vertical coordinates
         za(kloop)=z_gcm(iloop,jloop,kloop,tloop)
         ! store log values along altitude
         logq(kloop)=log(gcmdata(iloop,jloop,kloop,tloop))
       enddo !kloop
-      
+
       ! Interpolation
       do kloop=1,newaltlen
@@ -2382,5 +2383,5 @@
     enddo !iloop
   enddo !jloop
- enddo !tloop 
+ enddo !tloop
 endif
 
@@ -2402,5 +2403,5 @@
 ! grid values of altitude are known 'z_gcm', onto a new altitude grid 'z_new'.
 ! "Altitudes" must be above local surface altitudes.
-! Notes: 
+! Notes:
 !       'z_gcm' and 'znew' altitudes must be monotone increasing sequences.
 !       If altitudes in 'znew(i)' fall below those in 'z_gcm(:,:,1,:)', then
@@ -2408,5 +2409,5 @@
 !       flag=0, and extrapolated (exponentially) if flag=1.
 !       If altitudes in 'znew(i)' are above those in 'z_gcm(:,:,altlen,:)', then
-!       'newdata(:,:,i,:)' is set to constant value 'gcmdata(:,:,altlen,:)' 
+!       'newdata(:,:,i,:)' is set to constant value 'gcmdata(:,:,altlen,:)'
 !       if flag=0, and extrapolated (exponentially) if flag=1.
 !==============================================================================
@@ -2452,10 +2453,10 @@
 real :: x,y ! temporary variables
 integer :: iloop,jloop,kloop,tloop
-real,parameter :: g0=0.62
+real,parameter :: g0=0.6198
 !g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
-real,parameter :: a0=1187.E3
+real,parameter :: a0=1188.E3
 !a0: 'mean' Mars radius=3396.km
 real,parameter :: Rmean=296.9 ! molecular gas constant
-real :: gz 
+real :: gz
 ! gz: gravitational acceleration at a given (above areoid) altitude
 
@@ -2499,5 +2500,5 @@
         q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
       enddo !kloop
-      
+
       ! Interpolation
       do kloop=1,newaltlen
@@ -2516,5 +2517,5 @@
     enddo !iloop
   enddo !jloop
- enddo !tloop 
+ enddo !tloop
 else ! case when flag=1 (i.e.: rho or pressure)
  do tloop=1,tlen
@@ -2523,10 +2524,10 @@
       ! preliminary stuff
       do kloop=1,altlen
-        ! extract the vertical coordinates 
+        ! extract the vertical coordinates
         z(kloop)=z_gcm(iloop,jloop,kloop,tloop)
         ! store log values along altitude
         logq(kloop)=log(gcmdata(iloop,jloop,kloop,tloop))
       enddo !kloop
-      
+
       ! Interpolation
       do kloop=1,newaltlen
@@ -2553,5 +2554,5 @@
     enddo !iloop
   enddo !jloop
- enddo !tloop 
+ enddo !tloop
 endif
 
@@ -2714,5 +2715,5 @@
         double precision x,xi,sum
 
-	  ae= 1187000.d0
+	  ae= 1188000.d0
 	  gm =869.6 !stern et al 2015 gravitational parameter mu (km3/s-2) = G*M
 	  omega=1.138594423d-5    ! Pluto = 2*pi / P   rad/s
@@ -9296,5 +9297,5 @@
 	enddo
 	sum = sum+1.
-! centrifugal potential	
+! centrifugal potential
 	v0=(gm*sum/r + 0.5*omega**2 * r**2)
 !	write(out,*)'v0=', v0
@@ -9317,5 +9318,5 @@
         double precision dlon,dlat
         double precision rg
-        
+
         double precision pi,d2r
 	parameter (pi=3.141592653589792D0, d2r=pi/180.d0)
@@ -9358,5 +9359,5 @@
 	  enddo
 	  sum=sum+1.
-! centrifugal potential	
+! centrifugal potential
 	  rg=(gm*sum+0.5*(omega**2)*(r**3)*(cslt**2))/v0
 !          write(out,*) i,r,rg
