Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/caldyn0.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/caldyn0.F	(revision 1644)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/caldyn0.F	(revision 1644)
@@ -0,0 +1,1 @@
+link ../../dyn3d/caldyn0.F
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/callphysiq_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/callphysiq_mod.F90	(revision 1644)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/callphysiq_mod.F90	(revision 1644)
@@ -0,0 +1,95 @@
+!
+! $Id: $
+!
+MODULE callphysiq_mod
+
+IMPLICIT NONE
+
+CONTAINS
+
+SUBROUTINE call_physiq(klon,llm,nqtot,tname,                              &
+                       debut_split,lafin_split,                           &
+                       jD_cur,jH_cur_split,zdt_split,                     &
+                       zplev_omp,zplay_omp,                               &
+                       zpk_omp,zphi_omp,zphis_omp,                        &
+                       presnivs_omp,                                      &
+                       zufi_omp,zvfi_omp,zrfi_omp,ztfi_omp,zqfi_omp,      &
+                       flxwfi_omp,pducov,                                 &
+                       zdufi_omp,zdvfi_omp,zdtfi_omp,zdqfi_omp,zdpsrf_omp)
+
+  USE control_mod, ONLY: planet_type
+  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
+  USE physiq_mod, ONLY: physiq
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: klon ! (local) number of atmospheric columns
+  INTEGER,INTENT(IN) :: llm  ! number of atmospheric layers
+  INTEGER,INTENT(IN) :: nqtot ! number of tracers
+  CHARACTER(len=*),INTENT(IN) :: tname(nqtot) ! tracer names
+  LOGICAL,INTENT(IN) :: debut_split ! .true. if very first call to physics
+  LOGICAL,INTENT(IN) :: lafin_split ! .true. if last call to physics
+  REAL,INTENT(IN) :: JD_cur ! Julian day
+  REAL,INTENT(IN) :: JH_cur_split ! Julian hour (fraction of day)
+  REAL,INTENT(IN) :: zdt_split ! time step over which the physics are evaluated
+  REAL,INTENT(IN) :: zplev_omp(klon,llm+1) ! interlayer pressure (Pa)
+  REAL,INTENT(IN) :: zplay_omp(klon,llm) ! mid-layer pressure (Pa)
+  REAL,INTENT(IN) :: zpk_omp(klon,llm)
+  REAL,INTENT(IN) :: zphi_omp(klon,llm) ! geopotential at midlayer
+  REAL,INTENT(IN) :: zphis_omp(klon) ! surface geopotential
+  REAL,INTENT(IN) :: presnivs_omp(llm) ! approximate pressure of atm. layers
+  REAL,INTENT(IN) :: zufi_omp(klon,llm) ! zonal wind (m/s)
+  REAL,INTENT(IN) :: zvfi_omp(klon,llm) ! meridional wind (m/s)
+  REAL,INTENT(IN) :: zrfi_omp(klon,llm) ! relative wind vorticity, in s-1
+  REAL,INTENT(IN) :: ztfi_omp(klon,llm) ! temperature (K)
+  REAL,INTENT(IN) :: zqfi_omp(klon,llm,nqtot) ! tracers (*/kg of air)
+  REAL,INTENT(IN) :: flxwfi_omp(klon,llm) ! Vertical mass flux on lower mesh interfaces (kg/s) 
+  REAL,INTENT(IN) :: pducov(nbp_lon+1,nbp_lat,llm) ! dynamical tendency on ucov
+  ! tendencies (in */s) from the physics:
+  REAL,INTENT(OUT) :: zdufi_omp(klon,llm) ! tendency on zonal winds
+  REAL,INTENT(OUT) :: zdvfi_omp(klon,llm) ! tendency on meridional winds
+  REAL,INTENT(OUT) :: zdtfi_omp(klon,llm) ! tendency on temperature
+  REAL,INTENT(OUT) :: zdqfi_omp(klon,llm,nqtot) ! tendency on tracers
+  REAL,INTENT(OUT) :: zdpsrf_omp(klon) ! tendency on surface pressure
+
+  ! Local variables
+  CHARACTER(len=11) :: modname="call_physiq"
+  LOGICAL,SAVE :: firstcall=.true.
+!$OMP THREADPRIVATE(firstcall)
+
+! Sanity check on physics package type
+  IF (firstcall) THEN
+    IF (planet_type.ne."titan") THEN
+      CALL abort_gcm(modname,"wrong planet_type for this physics package",1)
+    ENDIF
+    firstcall=.false.
+  ENDIF
+
+
+! Call physics package with required inputs/outputs
+  CALL physiq(klon,           & ! ngrid
+              llm,            & ! nlayer
+              nqtot,          & ! nq
+              tname,          & ! nametrac
+              debut_split,    & ! firstcall
+              lafin_split,    & ! lastcall
+              jD_cur,         & ! pday
+              jH_cur_split,   & ! ptime
+              zdt_split,      & ! ptimestep
+              zplev_omp,      & ! pplev
+              zplay_omp,      & ! pplay
+              zphi_omp,       & ! pphi
+              zufi_omp,       & ! pu
+              zvfi_omp,       & ! pv
+              ztfi_omp,       & ! pt
+              zqfi_omp,       & ! pq
+              flxwfi_omp,     & ! flxw
+              zdufi_omp,      & ! pdu
+              zdvfi_omp,      & ! pdv
+              zdtfi_omp,      & ! pdt
+              zdqfi_omp,      & ! pdq
+              zdpsrf_omp)       ! pdpsrf
+
+
+END SUBROUTINE call_physiq
+
+END MODULE callphysiq_mod
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/datareadnc.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/datareadnc.F	(revision 1644)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/datareadnc.F	(revision 1644)
@@ -0,0 +1,292 @@
+c=======================================================================
+      SUBROUTINE datareadnc(relief,filename,phisinit,alb,ith,
+     .                    zmea,zstd,zsig,zgam,zthe)
+c=======================================================================
+c
+c
+c   Author: F. Hourdin      01/1997
+c   -------
+c
+c   Object: To read data from Martian surface to use in a GCM
+c   ------                from NetCDF file "surface.nc"
+c
+c
+c   Arguments:
+c   ----------
+c
+c     Inputs:
+c     ------
+c
+c     Outputs:
+c     --------
+c
+c=======================================================================
+c   donnees ALBEDO, INERTIE THERMIQUE, RELIEF:
+c
+c       Ces donnees sont au format NetCDF dans le fichier "surface.nc"
+c
+c   360 valeurs en longitude (de -179.5 a 179.5)
+c   180 valeurs en latitudes (de 89.5 a -89.5)
+c
+c   Pour les passer au format de la grille, on utilise "interp_horiz.F"
+c
+c   Il faut donc que ces donnees soient au format grille scalaire
+c               (imold+1 jmold+1)
+c       avec passage des coordonnees de la "boite" (rlonu, rlatv)
+c
+c   On prend imd (d pour donnees!) 
+c           imd = 360 avec copie de la 1ere valeur sur la imd+1 
+c                   (rlonud de -179 a -181)
+c           jmd = 179 
+c                   (rlatvd de 89 a -89)
+c=======================================================================
+
+      use datafile_mod, only: datadir, surfdir
+! to use  'getin'
+      USE ioipsl_getincom 
+      USE comconst_mod, ONLY: g,pi
+      implicit none
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comgeom.h"
+#include "netcdf.inc"
+
+c=======================================================================
+c   Declarations:
+C=======================================================================
+
+      INTEGER    imd,jmd,imdp1,jmdp1
+      parameter    (imd=360,jmd=179,imdp1=361,jmdp1=180)
+
+      INTEGER    iimp1
+      parameter    (iimp1=iim+1-1/iim)
+
+      character(len=3),intent(inout) :: relief*3
+      character(len=*),intent(in) :: filename ! surface.nc file
+      real,intent(out) :: phisinit(iimp1*jjp1) ! surface geopotential
+      real,intent(out) :: alb(iimp1*jjp1) ! albedo
+      real,intent(out) :: ith(iimp1*jjp1) ! thermal inertia
+      real,intent(out) :: zmea(imdp1*jmdp1)
+      real,intent(out) :: zstd(imdp1*jmdp1)
+      real,intent(out) :: zsig(imdp1*jmdp1)
+      real,intent(out) :: zgam(imdp1*jmdp1)
+      real,intent(out) :: zthe(imdp1*jmdp1)
+      
+      REAL        zdata(imd*jmdp1)
+      REAL        zdataS(imdp1*jmdp1)
+      REAL        pfield(iimp1*jjp1)
+
+      INTEGER   ierr
+
+      INTEGER   unit,nvarid
+
+      INTEGER    i,j,k
+
+      INTEGER klatdat,ngridmixgdat
+      PARAMETER (klatdat=180,ngridmixgdat=360)
+
+c    on passe une grille en rlonu rlatv et im+1 jm a interp_horiz)
+
+      REAL longitude(imd),latitude(jmdp1) ! Pour lecture des donnees
+      REAL rlonud(imdp1),rlatvd(jmd)
+
+      CHARACTER*20 string
+      DIMENSION string(4)
+!#include "fxyprim.h"
+
+      pi=2.*ASIN(1.)
+
+c=======================================================================
+c    rlonud, rlatvd
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c    Lecture NetCDF des donnees latitude et longitude
+c-----------------------------------------------------------------------
+      ierr = NF_OPEN (trim(datadir)//'/'//trim(surfdir)//'/'//
+     &                trim(adjustl(filename)),
+     &                NF_NOWRITE,unit)
+      IF (ierr.NE.NF_NOERR) THEN
+        ! In ye old days this file was stored in datadir;
+        ! let's be retro-compatible
+        ierr = NF_OPEN (trim(datadir)//'/'//
+     &                trim(adjustl(filename)),
+     &                NF_NOWRITE,unit)
+        
+      ENDIF
+      IF (ierr.NE.NF_NOERR) THEN
+        write(*,*)'Error : cannot open file '//trim(filename)
+        write(*,*)'(in phystd/datareadnc.F)'
+        write(*,*)'It should be in :',trim(datadir),'/',trim(surfdir)
+        write(*,*)'Check that your path to datagcm:',trim(datadir)
+        write(*,*)' is correct. You can change it in callphys.def with:'
+        write(*,*)' datadir = /absolute/path/to/datagcm'
+        write(*,*)'If necessary surface.nc (and other datafiles)'
+        write(*,*)' can be obtained online on:'
+        write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/'//
+     &             'LMDZ.GENERIC/datagcm/'
+        STOP
+      ENDIF
+
+c
+c Lecture des latitudes (coordonnees):
+c
+      ierr = NF_INQ_VARID (unit, "latitude", nvarid)
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(unit, nvarid, latitude)
+#else
+      ierr = NF_GET_VAR_REAL(unit, nvarid, latitude)
+#endif
+c
+c Lecture des longitudes (coordonnees):
+c
+      ierr = NF_INQ_VARID (unit, "longitude", nvarid)
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(unit, nvarid, longitude)
+#else
+      ierr = NF_GET_VAR_REAL(unit, nvarid, longitude)
+#endif
+
+c-----------------------------------------------------------------------
+c    Passage au format boites scalaires
+c-----------------------------------------------------------------------
+
+c-----------------------------------------------------------------------
+c       longitude(imd)        -->      rlonud(imdp1) 
+c-----------------------------------------------------------------------
+
+c Passage en coordonnees boites scalaires et en radian
+      do i=1,imd 
+          rlonud(i)=(longitude(i)+.5)*pi/180.
+      enddo
+
+c Repetition de la valeur im+1
+      rlonud(imdp1)=rlonud(1) + 2*pi
+
+c-----------------------------------------------------------------------
+c        latitude(jmdp1)         -->        rlonvd(jmd)
+c-----------------------------------------------------------------------
+
+c Passage en coordonnees boites scalaires et en radian
+      do j=1,jmd 
+          rlatvd(j)=(latitude(j)-.5)*pi/180.
+      enddo
+
+c=======================================================================
+c   lecture NetCDF de albedo, thermal, relief, zdtm (pour francois Lott)
+c=======================================================================
+
+      string(1) = 'albedo'
+      string(2) = 'thermal'
+      if (relief.ne.'pla') then
+        write(*,*) ' La topographie est celle de MOLA'
+        relief = 'MOL'
+          string(3) = 'z'//relief
+      else
+          string(3) = 'zMOL'  ! pour qu''il lise qqchose sur le fichier
+                            ! remise a 0 derriere
+      endif
+      string(4) = 'zMOL'    ! lecture pour calcul topog. sous-maille
+ 
+
+      DO k=1,4
+          write(*,*) 'string',k,string(k)
+          
+c-----------------------------------------------------------------------
+c    initialisation
+c-----------------------------------------------------------------------
+      pfield(1:iimp1*jjp1)=0
+      zdata(1:imd*jmdp1)=0
+      zdataS(1:iimp1*jjp1)=0
+
+c-----------------------------------------------------------------------
+c    Lecture NetCDF  
+c-----------------------------------------------------------------------
+
+      ierr = NF_INQ_VARID (unit, string(k), nvarid)
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(unit, nvarid, zdata)
+#else
+      ierr = NF_GET_VAR_REAL(unit, nvarid, zdata)
+#endif
+
+c-----------------------------------------------------------------------
+c        Cas particulier "Francois Lott" ( k=4 ) (relief sous-maille)
+c-----------------------------------------------------------------------
+      if (k.eq.4) then
+
+          zdata(:)=1000.*zdata(:)
+          longitude(:)=(pi/180.)*longitude(:)
+          latitude(:)=(pi/180.)*latitude(:)
+
+          call grid_noro1(360, 180, longitude, latitude, zdata,
+     .         iim, jjp1, rlonv, rlatu, zmea,zstd,zsig,zgam,zthe)
+
+          !CALL dump2d(iip1,jjp1,zmea,'zmea')
+          !CALL dump2d(iip1,jjp1,zstd,'zstd')
+          !CALL dump2d(iip1,jjp1,zsig,'zsig')
+          !CALL dump2d(iip1,jjp1,zgam,'zgam')
+          !CALL dump2d(iip1,jjp1,zthe,'zthe')
+
+      endif
+
+c-----------------------------------------------------------------------
+c   Passage de zdata en grille (imdp1 jmdp1)
+c-----------------------------------------------------------------------
+      do j=1,jmdp1
+          do i=1,imd
+              zdataS(i+imdp1*(j-1)) = zdata(i+ngridmixgdat*(j-1))
+          enddo
+          zdataS(imdp1+imdp1*(j-1)) = zdata(1+ngridmixgdat*(j-1))
+      enddo
+
+c-----------------------------------------------------------------------
+c    Interpolation
+c-----------------------------------------------------------------------
+      call interp_horiz(zdataS,pfield,imd,jmd,
+     .    iim, jjm,1,rlonud,rlatvd,rlonu,rlatv) 
+
+c-----------------------------------------------------------------------
+c    Periodicite    
+c-----------------------------------------------------------------------
+
+      do j=1,jjp1
+         pfield(iimp1*j) =  pfield(1+iimp1*(j-1))
+      enddo 
+ 
+c-----------------------------------------------------------------------
+c    Sauvegarde des champs    
+c-----------------------------------------------------------------------
+
+      if (k.eq.1) then                    ! albedo
+         do i=1,iimp1*jjp1
+              alb(i) = pfield(i)
+          enddo
+      elseif (k.eq.2) then                ! thermal
+         do i=1,iimp1*jjp1
+              ith(i) = pfield(i)
+          enddo
+      elseif (k.eq.3) then                ! relief
+        if (relief.eq.'pla') then
+             phisinit(1:iimp1*jjp1)=0
+        else
+             phisinit(1:iimp1*jjp1)=pfield(1:iimp1*jjp1)
+        endif
+      endif
+
+      ENDDO
+
+c-----------------------------------------------------------------------
+c    Traitement Phisinit
+c-----------------------------------------------------------------------
+
+      phisinit(1:iimp1*jjp1)=1000.*phisinit(1:iimp1*jjp1)
+      !CALL dump2d(iimp1,jjp1,phisinit,'Altitude en m')
+      phisinit(:)=g*phisinit(:)
+
+c-----------------------------------------------------------------------
+c    FIN
+c-----------------------------------------------------------------------
+
+      END
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/defrun_new.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/defrun_new.F	(revision 1644)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/defrun_new.F	(revision 1644)
@@ -0,0 +1,1 @@
+link ../../dyn3d/defrun_new.F
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/exner_hyb.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/exner_hyb.F	(revision 1644)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/exner_hyb.F	(revision 1644)
@@ -0,0 +1,1 @@
+link ../../dyn3d/exner_hyb.F
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/grid_noro1.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/grid_noro1.F	(revision 1644)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/grid_noro1.F	(revision 1644)
@@ -0,0 +1,425 @@
+      SUBROUTINE grid_noro1(imdep, jmdep, xdata, ydata, entree,
+     .                 imar, jmar, x, y, zmea,zstd,zsig,zgam,zthe)
+c=======================================================================
+c (F. Lott) (voir aussi z.x. Li, A. Harzallah et L. Fairhead)
+c
+c      Calcul des parametres de l'orographie sous-maille necessaires
+c      au nouveau shema de representation des montagnes meso-echelles
+c      dans le modele.  Les points sont mis sur une grille rectangulaire
+c      pseudo-physique.  Typiquement, il y a iim+1 latitudes incluant
+c      le pole nord et le pole sud.  Il y a jjm+1 longitudes, y compris
+c      aux poles.  Aux poles les champs peuvent ont une valeurs repetee
+c      jjm+1 fois.....  La valeur du champs en jjm+1 (jmar) est celle
+c      en j=1.  
+c      Les parametres a,b,c,d representent les limites de la region
+c      de point de grille correspondant a un point decrit precedemment.
+c      Les moyennes sur ces regions des valeurs calculees a partir de
+c      l'USN, sont ponderees par un poids, fonction de la surface
+c      occuppe par ces donnees a l'interieure de la grille du modele.
+c      Dans la plupart des cas ce poid est le rapport entre la surface
+c      de la region de point de grille USN et la surface de la region
+c      de point de grille du modele.
+c       
+c
+c           (c)
+c        ----d-----
+c        | . . . .|
+c        |        |
+c     (b)a . * . .b(a)
+c        |        |
+c        | . . . .|
+c        ----c-----
+c           (d)
+C=======================================================================
+c INPUT:
+c        imdep, jmdep: dimensions X et Y pour depart
+c        xdata, ydata: coordonnees X et Y pour depart
+c        entree: champ d'entree a transformer
+c        dans ce programme, on assume que les donnees sont les altitudes
+c        de l'USNavy: imdep=iusn=2160, jmdep=jusn=1080.
+c OUTPUT:
+c        imar, jmar: dimensions X et Y d'arrivee
+c        x, y: coordonnees X et Y d'arrivee
+c        les champs de sorties sont sur une grille physique:
+c             zmea:  orographie moyenne
+c             zstd:  deviation standard de l'orographie sous-maille
+c             zsig:  pente de l'orographie sous-maille 
+c             zgam:  anisotropy de l'orographie sous maille
+c             zthe:  orientation de l'axe oriente dans la direction
+c                    de plus grande pente de l'orographie sous maille
+C=======================================================================
+c     IMPLICIT INTEGER (I,J)
+c     IMPLICIT REAL(X,Z) 
+
+       USE comconst_mod, ONLY: rad
+
+       implicit none
+       integer iusn,jusn,iext
+       parameter(iusn=360,jusn=180,iext=40)
+c!-*-      include 'param1'
+c!-*-      include 'comcstfi.h'
+#include "dimensions.h"
+c!-*-
+c!-*-      parameter(iim=cols,jjm=rows)
+      REAL xusn(iusn+2*iext),yusn(jusn+2)	
+      REAL zusn(iusn+2*iext,jusn+2),zusnfi(iusn+2*iext,jusn+2)
+
+c   modif declarations pour implicit none
+      real zmeanor,zmeasud,zstdnor,zstdsud,zsignor
+      real zsigsud,zweinor,zweisud
+      real xk,xl,xm,xw,xp,xq
+      real zmaxmea,zmaxstd,zmaxsig,zmaxgam,zmaxthe,zminthe
+      real zbordnor,zbordsud,zbordest,zbordoue,xpi
+      real zdeltax,zdeltay,zlenx,zleny,weighx,weighy,xincr
+      integer i,j,ii,jj,ideltax,ihalph
+
+      INTEGER imdep, jmdep
+      REAL xdata(imdep),ydata(jmdep) 
+      REAL entree(imdep,jmdep)
+c
+      INTEGER imar, jmar
+  
+      REAL ztz(iim+1,jjm+1),zxtzx(iim+1,jjm+1)
+      REAL zytzy(iim+1,jjm+1),zxtzy(iim+1,jjm+1)
+      REAL zxtzxusn(iusn+2*iext,jusn+2),zytzyusn(iusn+2*iext,jusn+2)
+      REAL zxtzyusn(iusn+2*iext,jusn+2)
+      REAL weight(iim+1,jjm+1)
+      REAL x(imar+1),y(jmar)
+      REAL zmea(imar+1,jmar),zstd(imar+1,jmar)
+      REAL zsig(imar+1,jmar),zgam(imar+1,jmar),zthe(imar+1,jmar)
+c
+      REAL a(2200),b(2200),c(1100),d(1100)
+c
+c  quelques constantes:
+c
+      print *,' parametres de l orographie a l echelle sous maille' 
+      print*,'rad =',rad
+      print*,'Long et lat entree'
+      print*,(x(i),i=1,imar+1)
+      print*,(y(j),j=1,jmar)
+       print*,'Long et lat donnees'
+      print*,(xdata(i),i=1,imdep)
+      print*,(ydata(j),j=1,jmdep)
+
+      xpi=acos(-1.)
+      zdeltay=2.*xpi/float(jusn)*rad
+c
+c  quelques tests de dimensions:
+c    
+      IF (imar.GT.2200 .OR. jmar.GT.1100) THEN
+         PRINT*, 'imar ou jmar trop grand', imar, jmar
+         CALL ABORT
+      ENDIF
+
+      IF(imdep.ne.iusn.or.jmdep.ne.jusn)then
+         print *,' imdep ou jmdep mal dimensionnes:',imdep,jmdep
+         call abort
+      ENDIF
+
+      IF(imar+1.gt.iim+1.or.jmar.gt.jjm+1)THEN
+        print *,' imar ou jmar mal dimensionnes:',imar,jmar
+        call abort
+      ENDIF
+c
+C  Extension de la base de donnee de l'USN pour faciliter
+C  les calculs ulterieurs:
+c
+      DO j=1,jusn
+        yusn(j+1)=ydata(j)
+      DO i=1,iusn
+        zusn(i+iext,j+1)=entree(i,j)
+        xusn(i+iext)=xdata(i)
+      ENDDO
+      DO i=1,iext
+        zusn(i,j+1)=entree(iusn-iext+i,j)
+        xusn(i)=xdata(iusn-iext+i)-2.*xpi
+        zusn(iusn+iext+i,j+1)=entree(i,j)
+        xusn(iusn+iext+i)=xdata(i)+2.*xpi
+      ENDDO
+      ENDDO
+
+        yusn(1)=ydata(1)+(ydata(1)-ydata(2))
+        yusn(jusn+2)=ydata(jusn)+(ydata(jusn)-ydata(jusn-1))
+       DO i=1,iusn/2+iext
+        zusn(i,1)=zusn(i+iusn/2,2)
+        zusn(i+iusn/2+iext,1)=zusn(i,2)
+        zusn(i,jusn+2)=zusn(i+iusn/2,jusn+1)
+        zusn(i+iusn/2+iext,jusn+2)=zusn(i,jusn+1)
+       ENDDO
+c
+c  Calcul d'une orographie filtree aux hautes latitudes
+c  pour permettre des calculs plus isotropiques sur la pente
+c  des montagnes
+c
+       DO i=1,IUSN+2*iext
+       DO J=1,JUSN+2
+          zusnfi(i,j)=0.0
+       ENDDO
+       ENDDO
+
+      DO j=1,jusn
+            ideltax=1./cos(yusn(j+1))
+            ideltax=min(iusn/2-1,ideltax)
+            IF(MOD(IDELTAX,2).EQ.0)THEN
+              IDELTAX=IDELTAX+1
+            ENDIF
+            IHALPH=(IDELTAX-1)/2 
+c           print *,' ideltax=',ideltax
+         IF(ideltax.eq.1)THEN
+            DO i=1,iusn
+               zusnfi(i+iext,j+1)=entree(i,j)
+            ENDDO   
+         ELSE
+            DO i=1,ihalph
+               DO ii=1,i+ihalph
+               zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)+entree(ii,j)
+               ENDDO
+               DO ii=ihalph-i,0,-1
+               zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)+entree(iusn-ii,j)
+               ENDDO  
+               zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)/float(ideltax)
+            ENDDO   
+            DO i=iusn-ihalph+1,iusn
+               DO ii = i-ihalph,iusn
+               zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)+entree(ii,j)
+               ENDDO 
+               DO ii = 1,ihalph+i-iusn
+               zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)+entree(ii,j)
+               ENDDO
+               zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)/float(ideltax)
+            ENDDO
+            DO i=ihalph+1,iusn-ihalph
+               DO ii=-ihalph,ihalph
+               zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)+entree(i+ii,j)
+               ENDDO
+               zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)/float(ideltax)
+            ENDDO
+         ENDIF
+            DO i=1,iext
+               zusnfi(i,j+1)=zusnfi(iusn-iext+i,j+1)
+               zusnfi(i+iusn+iext,j+1)=zusnfi(i,j+1)
+            ENDDO
+      ENDDO
+c  
+c Calculer les limites des zones des nouveaux points
+c
+      a(1) = x(1) - (x(2)-x(1))/2.0
+      b(1) = (x(1)+x(2))/2.0
+      DO i = 2, imar-1
+         a(i) = b(i-1)
+         b(i) = (x(i)+x(i+1))/2.0
+      ENDDO
+      a(imar) = b(imar-1)
+      b(imar) = x(imar) + (x(imar)-x(imar-1))/2.0
+
+      c(1) = y(1) - (y(2)-y(1))/2.0
+      d(1) = (y(1)+y(2))/2.0
+      DO j = 2, jmar-1
+         c(j) = d(j-1)
+         d(j) = (y(j)+y(j+1))/2.0
+      ENDDO
+      c(jmar) = d(jmar-1)
+      d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0
+c
+c      quelques initialisations:
+      print*,'OKM1'
+c
+      DO i = 1, imar
+      DO j = 1, jmar
+         weight(i,j) = 0.0
+         zxtzx(i,j) = 0.0
+         zytzy(i,j) = 0.0
+         zxtzy(i,j) = 0.0
+         ztz(i,j) = 0.0
+         zmea(i,j) = 0.0
+         zstd(i,j)=0.0
+      ENDDO
+      ENDDO
+c
+c  calculs des correlations de pentes sur la grille de l'USN.
+c
+         DO j = 2,jusn+1 
+         DO i = 1, iusn+2*iext
+            zytzyusn(i,j)=0.0
+            zxtzxusn(i,j)=0.0
+            zxtzyusn(i,j)=0.0
+         ENDDO
+         ENDDO
+
+
+         DO j = 2,jusn+1 
+            zdeltax=zdeltay*cos(yusn(j))
+         DO i = 2, iusn+2*iext-1
+            zytzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))**2/zdeltay**2
+            zxtzxusn(i,j)=(zusnfi(i+1,j)-zusnfi(i-1,j))**2/zdeltax**2
+            zxtzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))/zdeltay
+     *                   *(zusnfi(i+1,j)-zusnfi(i-1,j))/zdeltax
+         ENDDO
+
+         ENDDO
+
+ 
+
+      print*,'OK0'
+c
+c  sommations des differentes quantites definies precedemment
+c  sur une grille du modele.
+c 
+      zleny=xpi/float(jusn)*rad
+      xincr=xpi/2./float(jusn)
+       DO ii = 1, imar
+       DO jj = 1, jmar
+c        PRINT *,' iteration ii jj:',ii,jj
+         DO j = 2,jusn+1 
+c         DO j = 3,jusn 
+            zlenx=zleny*cos(yusn(j))
+            zdeltax=zdeltay*cos(yusn(j))
+            zbordnor=(c(jj)-yusn(j)+xincr)*rad
+            zbordsud=(yusn(j)-d(jj)+xincr)*rad
+            weighy=amax1(0.,
+     *             amin1(zbordnor,zbordsud,zleny))
+         IF(weighy.ne.0)THEN
+         DO i = 2, iusn+2*iext-1
+            zbordest=(xusn(i)-a(ii)+xincr)*rad*cos(yusn(j))
+            zbordoue=(b(ii)+xincr-xusn(i))*rad*cos(yusn(j))
+            weighx=amax1(0.,
+     *             amin1(zbordest,zbordoue,zlenx))
+            IF(weighx.ne.0)THEN
+            weight(ii,jj)=weight(ii,jj)+weighx*weighy
+            zxtzx(ii,jj)=zxtzx(ii,jj)+zxtzxusn(i,j)*weighx*weighy
+            zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)*weighx*weighy
+            zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)*weighx*weighy
+            ztz(ii,jj)  =ztz(ii,jj)  +zusn(i,j)*zusn(i,j)*weighx*weighy
+            zmea(ii,jj) =zmea(ii,jj)+zusn(i,j)*weighx*weighy
+            ENDIF
+         ENDDO
+         ENDIF
+         ENDDO
+       ENDDO
+       ENDDO
+c
+c  calculs des differents parametres necessaires au programme
+c  de parametrisation de l'orographie a l'echelle moyenne:
+c
+      zmaxmea=0.
+      zmaxstd=0.
+      zmaxsig=0.
+      zmaxgam=0.
+      zmaxthe=0.
+      zminthe=0.
+c     print 100,' '
+c100  format(1X,A1,'II JJ',4X,'H',8X,'SD',8X,'SI',3X,'GA',3X,'TH') 
+       print*,'OK1'
+       DO ii = 1, imar
+       DO jj = 1, jmar
+c       print*,'ok0'
+         IF (weight(ii,jj) .NE. 0.0) THEN
+c  Orography moyenne:
+c         print*,'ok1'
+           zmea (ii,jj)=zmea (ii,jj)/weight(ii,jj)
+           zxtzx(ii,jj)=zxtzx(ii,jj)/weight(ii,jj)
+           zytzy(ii,jj)=zytzy(ii,jj)/weight(ii,jj)
+           zxtzy(ii,jj)=zxtzy(ii,jj)/weight(ii,jj)
+           ztz(ii,jj)  =ztz(ii,jj)/weight(ii,jj)
+c         print*,'ok2'
+c  Deviation standard:
+           zstd(ii,jj)=sqrt(amax1(0.,ztz(ii,jj)-zmea(ii,jj)**2))
+c  Coefficients K, L et M:
+           xk=(zxtzx(ii,jj)+zytzy(ii,jj))/2.
+           xl=(zxtzx(ii,jj)-zytzy(ii,jj))/2.
+           xm=zxtzy(ii,jj)
+           xp=xk-sqrt(xl**2+xm**2)
+           xq=xk+sqrt(xl**2+xm**2)
+           xw=1.e-8
+           if(xp.le.xw) xp=0.
+           if(xq.le.xw) xq=xw
+           if(abs(xm).le.xw) xm=xw*sign(1.,xm)
+c          print*,'ok3'
+c pente: 
+           zsig(ii,jj)=sqrt(xq)
+c           zsig(ii,jj)=sqrt(2.*xk)
+c isotropy:
+           zgam(ii,jj)=xp/xq
+c angle theta:
+           zthe(ii,jj)=57.29577951*atan2(xm,xl)/2.
+
+c          print 101,ii,jj,
+c    *           zmea(ii,jj),zstd(ii,jj),zsig(ii,jj),zgam(ii,jj),
+c    *           zthe(ii,jj)
+c101  format(1x,2(1x,i2),2(1x,f7.1),1x,f7.4,2x,f4.2,1x,f5.1)     
+c          print*,'ok4'
+         ELSE
+c           PRINT*, 'probleme,ii,jj=', ii,jj
+c          print*,'ok1b'
+         ENDIF
+      zmaxmea=amax1(zmea(ii,jj),zmaxmea)
+c         print*,'oka'
+      zmaxstd=amax1(zstd(ii,jj),zmaxstd)
+c         print*,'okb'
+      zmaxsig=amax1(zsig(ii,jj),zmaxsig)
+c         print*,'okc'
+      zmaxgam=amax1(zgam(ii,jj),zmaxgam)
+c         print*,'okd'
+      zmaxthe=amax1(zthe(ii,jj),zmaxthe)
+c         print*,'oke'
+      zminthe=amin1(zthe(ii,jj),zminthe)
+c      print*,'ok5'
+       ENDDO
+       ENDDO
+
+      print *,'  MEAN ORO:',zmaxmea
+	  print *,'  ST. DEV.:',zmaxstd
+      print *,'  PENTE:',zmaxsig
+      print *,' ANISOTROP:',zmaxgam
+      print *,'  ANGLE:',zminthe,zmaxthe	
+      
+C
+c  On passe ce donnees sur la grille dite physique....(?)
+c  On met gamma et theta a 1. et 0. aux poles ou ces quantites
+c  n'ont pas vraiment de sens
+c
+      DO jj=1,jmar
+      zmea(imar+1,jj)=zmea(1,jj)
+      zstd(imar+1,jj)=zstd(1,jj)
+      zsig(imar+1,jj)=zsig(1,jj)
+      zgam(imar+1,jj)=zgam(1,jj)
+      zthe(imar+1,jj)=zthe(1,jj)
+      ENDDO
+
+
+      zmeanor=0.0
+      zmeasud=0.0
+      zstdnor=0.0
+      zstdsud=0.0
+      zsignor=0.0
+      zsigsud=0.0
+      zweinor=0.0
+      zweisud=0.0
+
+      DO ii=1,imar
+      zweinor=zweinor+              weight(ii,   1)
+      zweisud=zweisud+              weight(ii,jmar)
+      zmeanor=zmeanor+zmea(ii,   1)*weight(ii,   1)
+      zmeasud=zmeasud+zmea(ii,jmar)*weight(ii,jmar)
+      zstdnor=zstdnor+zstd(ii,   1)*weight(ii,   1)
+      zstdsud=zstdsud+zstd(ii,jmar)*weight(ii,jmar)
+      zsignor=zsignor+zsig(ii,   1)*weight(ii,   1)
+      zsigsud=zsigsud+zsig(ii,jmar)*weight(ii,jmar)
+      ENDDO
+
+      DO ii=1,imar+1
+      zmea(ii,   1)=zmeanor/zweinor
+      zmea(ii,jmar)=zmeasud/zweisud
+      zstd(ii,   1)=zstdnor/zweinor
+      zstd(ii,jmar)=zstdsud/zweisud
+      zsig(ii,   1)=zsignor/zweinor
+      zsig(ii,jmar)=zsigsud/zweisud
+      zgam(ii,   1)=1.
+      zgam(ii,jmar)=1.
+      zthe(ii,   1)=0.
+      zthe(ii,jmar)=0.
+      ENDDO
+
+
+      RETURN
+      END
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/ini_archive.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/ini_archive.F	(revision 1644)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/ini_archive.F	(revision 1644)
@@ -0,0 +1,529 @@
+c=======================================================================
+      subroutine ini_archive(nid,idayref,phis,ith,tab_cntrl_fi,
+     &                       tab_cntrl_dyn)
+c=======================================================================
+c
+c
+c   Date:    01/1997
+c   ----
+c
+c   Objet:  ecriture de l'entete du fichier "start_archive"
+c   -----
+c
+c	 Proche de iniwrite.F
+c
+c	 On ajoute dans le tableau "tab_cntrl" (dynamique), a partir de 51, 
+c	 les valeurs de tab_cntrl_fi (les 38 parametres de controle physiques
+c	 du RUN + ptotal et cotoicetotal)
+c
+c			tab_cntrl(50+l)=tab_cntrl_fi(l)
+c
+c   Arguments:
+c   ---------
+c
+c	Inputs:
+c   ------
+c
+c       nid            unite logique du fichier "start_archive"
+c       idayref        Valeur du jour initial a mettre dans
+c                      l'entete du fichier "start_archive"
+c       phis           geopotentiel au sol
+c       ith            soil thermal inertia
+c       tab_cntrl_fi   tableau des param physiques
+c
+
+c=======================================================================
+
+      USE comsoil_h
+      USE slab_ice_h, only: noceanmx
+!      use control_mod
+      USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt
+      USE comconst_mod, ONLY: daysec,dtvr,rad,omeg,g,kappa,pi
+      implicit none
+
+#include "dimensions.h"
+!#include "dimphys.h"
+#include "paramet.h"
+#include "comgeom.h"
+!#include "control.h"
+
+#include "netcdf.inc"
+
+c-----------------------------------------------------------------------
+c   Declarations
+c-----------------------------------------------------------------------
+
+c   Local:
+c   ------
+      INTEGER	length,l
+      parameter (length = 100)
+      REAL		tab_cntrl(length) ! tableau des parametres du run
+      INTEGER	loop
+      INTEGER	ierr, setvdim, putvdim, putdat, setname,cluvdb
+      INTEGER	setdim
+      INTEGER	ind1,indlast
+
+c   Arguments:
+c   ----------
+      INTEGER*4	idayref
+      REAL		phis(ip1jmp1)
+      real ith(ip1jmp1,nsoilmx)
+      REAL		tab_cntrl_fi(length)
+      REAL tab_cntrl_dyn(length)
+
+!Mars --------Ajouts-----------
+c   Variables locales pour NetCDF:
+c
+      INTEGER dims2(2), dims3(3) !, dims4(4)
+      INTEGER idim_index
+      INTEGER idim_rlonu, idim_rlonv, idim_rlatu, idim_rlatv
+      INTEGER idim_llmp1,idim_llm
+      INTEGER idim_tim
+      INTEGER idim_nsoilmx ! "subsurface_layers" dimension ID #
+      INTEGER idim_noceanmx ! "ocean_layers" dimension ID #
+      INTEGER nid,nvarid
+      real sig_s(llm),s(llm)
+
+      pi  = 2. * ASIN(1.)
+
+
+c-----------------------------------------------------------------------
+c   Remplissage du tableau des parametres de controle du RUN  (dynamique)
+c-----------------------------------------------------------------------
+
+      DO l=1,length
+         tab_cntrl(l)=0.
+      ENDDO
+      
+      tab_cntrl(1:50)=tab_cntrl_dyn(1:50)
+
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!      tab_cntrl(1)  = FLOAT(iim) ! nombre de points en longitude
+!      tab_cntrl(2)  = FLOAT(jjm) ! nombre de points en latitude
+!      tab_cntrl(3)  = FLOAT(llm) ! nombre de couches
+!      tab_cntrl(4)  = FLOAT(idayref) ! jour 0
+!      tab_cntrl(5)  = rad ! rayon de mars(m) ~3397200
+!      tab_cntrl(6)  = omeg ! vitesse de rotation (rad.s-1)
+!      tab_cntrl(7)  = g   ! gravite (m.s-2) ~3.72
+!      tab_cntrl(8)  = cpp 
+!      tab_cntrl(8)  = 43.49 !mars temporaire Masse molaire de l''atm (g.mol-1) ~43.49
+!      tab_cntrl(9)  = kappa ! = r/cp  ~0.256793 (=rcp dans physique)
+!      tab_cntrl(10) = daysec ! duree du sol (s)  ~88775
+!      tab_cntrl(11) = dtvr ! pas de temps de la dynamique (s)
+!      tab_cntrl(12) = etot0 ! energie totale    !
+!      tab_cntrl(13) = ptot0 ! pression totalei   !    variables
+!      tab_cntrl(14) = ztot0 ! enstrophie totale   !  de controle
+!      tab_cntrl(15) = stot0 ! enthalpie totale   !    globales
+!      tab_cntrl(16) = ang0 ! moment cinetique  !
+!      tab_cntrl(17) = pa
+!      tab_cntrl(18) = preff
+
+c    .....    parametres  pour le zoom      ......   
+
+!      tab_cntrl(19)  = clon ! longitude en degres du centre du zoom
+!      tab_cntrl(20)  = clat ! latitude en degres du centre du zoom
+!      tab_cntrl(21)  = grossismx ! facteur de grossissement du zoom,selon longitude
+!      tab_cntrl(22)  = grossismy ! facteur de grossissement du zoom ,selon latitude
+
+!      IF ( fxyhypb )   THEN
+!       tab_cntrl(23) = 1.
+!       tab_cntrl(24) = dzoomx ! extension en longitude  de la zone du zoom
+!       tab_cntrl(25) = dzoomy ! extension en latitude  de la zone du zoom
+!      ELSE
+!       tab_cntrl(23) = 0.
+!       tab_cntrl(24) = dzoomx ! extension en longitude  de la zone du zoom
+!       tab_cntrl(25) = dzoomy ! extension en latitude  de la zone du zoom
+!       tab_cntrl(26) = 0.
+!       IF ( ysinus)  tab_cntrl(26) = 1.
+!      ENDIF
+
+c-----------------------------------------------------------------------
+c   Copie du tableau des parametres de controle du RUN  (physique)
+c		dans le tableau tab_cntrl
+c-----------------------------------------------------------------------
+
+      DO l=1,50
+         tab_cntrl(50+l)=tab_cntrl_fi(l)
+      ENDDO
+
+c=======================================================================
+c	Ecriture NetCDF de l''entete du fichier "start_archive"
+c=======================================================================
+
+c
+c Preciser quelques attributs globaux:
+c
+      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 21,
+     &                       "Fichier start_archive")
+c
+c Definir les dimensions du fichiers:
+c
+c     CHAMPS AJOUTES POUR LA VISUALISATION T,ps, etc... avec Grads ou ferret:
+      ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_rlatu)
+      ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_rlonv)
+      ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
+      ierr = NF_DEF_DIM (nid,"subsurface_layers",nsoilmx,idim_nsoilmx)
+      ierr = NF_DEF_DIM (nid,"ocean_layers",noceanmx,idim_noceanmx)
+
+      ierr = NF_DEF_DIM (nid,"index", length, idim_index)
+      ierr = NF_DEF_DIM (nid,"rlonu", iip1, idim_rlonu)
+      ierr = NF_DEF_DIM (nid,"rlatv", jjm, idim_rlatv)
+      ierr = NF_DEF_DIM (nid,"interlayer", llmp1, idim_llmp1)
+      ierr = NF_DEF_DIM (nid,"Time", NF_UNLIMITED, idim_tim)
+
+c
+      ierr = NF_ENDDEF(nid) ! sortir du mode de definition
+
+c-----------------------------------------------------------------------
+c  Ecriture du tableau des parametres du run
+c-----------------------------------------------------------------------
+
+      call def_var(nid,"Time","Time","days since 00:00:00",1,
+     .            idim_tim,nvarid,ierr)
+
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"controle",NF_DOUBLE,1,idim_index,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"controle",NF_FLOAT,1,idim_index,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22,
+     .                       "Parametres de controle")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
+#endif
+
+c-----------------------------------------------------------------------
+c  Ecriture des longitudes et latitudes
+c-----------------------------------------------------------------------
+
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"rlonu",NF_DOUBLE,1,idim_rlonu,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"rlonu",NF_FLOAT,1,idim_rlonu,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
+     .                       "Longitudes des points U")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonu)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonu)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"rlatu",NF_DOUBLE,1,idim_rlatu,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"rlatu",NF_FLOAT,1,idim_rlatu,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22,
+     .                       "Latitudes des points U")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"rlonv",NF_DOUBLE,1,idim_rlonv,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"rlonv",NF_FLOAT,1,idim_rlonv,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
+     .                       "Longitudes des points V")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"rlatv",NF_DOUBLE,1,idim_rlatv,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"rlatv",NF_FLOAT,1,idim_rlatv,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22,
+     .                       "Latitudes des points V")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatv)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatv)
+#endif
+
+c-----------------------------------------------------------------------
+c  Ecriture des niveaux verticaux
+c-----------------------------------------------------------------------
+
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"ap",NF_DOUBLE,1,idim_llmp1,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"ap",NF_FLOAT,1,idim_llmp1,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 32,
+     .                       "Coef A: niveaux pression hybride")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ap)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,ap)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"bp",NF_DOUBLE,1,idim_llmp1,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"bp",NF_FLOAT,1,idim_llmp1,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 35,
+     .                       "Coefficient B niveaux sigma hybride")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bp)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,bp)
+#endif
+c
+c ----------------------
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"aps",NF_DOUBLE,1,idim_llm,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"aps",NF_FLOAT,1,idim_llm,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 36,
+     .      "Coef AS: hybrid pressure in midlayers")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aps)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,aps)
+#endif
+c
+c ----------------------
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"bps",NF_DOUBLE,1,idim_llm,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"bps",NF_FLOAT,1,idim_llm,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 30,
+     .      "Coef BS: hybrid sigma midlayers")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bps)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,bps)
+#endif
+c
+c ----------------------
+
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"presnivs",NF_DOUBLE,1,idim_llm,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"presnivs",NF_FLOAT,1,idim_llm,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,presnivs)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,presnivs)
+#endif
+c ------------------------------------------------------------------
+c  Variable uniquement pour visualisation avec Grads ou Ferret
+c ------------------------------------------------------------------
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"latitude",NF_DOUBLE,1,idim_rlatu,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"latitude",NF_FLOAT,1,idim_rlatu,nvarid)
+#endif
+      ierr =NF_PUT_ATT_TEXT(nid,nvarid,'units',13,"degrees_north")
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 14,
+     .      "North latitude")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
+#endif
+c----------------------
+       ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr =NF_DEF_VAR(nid,"longitude", NF_DOUBLE, 1, idim_rlonv,nvarid)
+#else
+      ierr = NF_DEF_VAR(nid,"longitude", NF_FLOAT, 1, idim_rlonv,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 14,
+     .      "East longitude")
+      ierr = NF_PUT_ATT_TEXT(nid,nvarid,'units',12,"degrees_east")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180)
+#endif
+c--------------------------
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "altitude", NF_DOUBLE, 1,
+     .       idim_llm,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "altitude", NF_FLOAT, 1,
+     .       idim_llm,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",10,"pseudo-alt")
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,'units',2,"km")
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,'positive',2,"up")
+
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pseudoalt)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,pseudoalt)
+#endif
+
+!-------------------------------
+! (soil) depth variable mlayer() (known from comsoil.h)
+!-------------------------------
+      ierr=NF_REDEF (nid) ! Enter NetCDF (re-)define mode
+      ! define variable
+#ifdef NC_DOUBLE
+      ierr=NF_DEF_VAR(nid,"soildepth",NF_DOUBLE,1,idim_nsoilmx,nvarid)
+#else
+      ierr=NF_DEF_VAR(nid,"soildepth",NF_FLOAT,1,idim_nsoilmx,nvarid)
+#endif
+      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 20,
+     .                        "Soil mid-layer depth")
+      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",1,"m")
+      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"positive",4,"down")
+      ierr=NF_ENDDEF(nid) ! Leave NetCDF define mode
+      ! write variable
+#ifdef NC_DOUBLE
+      ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,mlayer)
+#else
+      ierr=NF_PUT_VAR_REAL (nid,nvarid,mlayer)
+#endif
+
+!---------------------
+! soil thermal inertia
+!---------------------
+      ierr=NF_REDEF (nid) ! Enter NetCDF (re-)define mode
+      dims3(1)=idim_rlonv
+      dims3(2)=idim_rlatu
+      dims3(3)=idim_nsoilmx
+      ! define variable
+#ifdef NC_DOUBLE
+      ierr=NF_DEF_VAR(nid,"inertiedat",NF_DOUBLE,3,dims3,nvarid)
+#else
+      ierr=NF_DEF_VAR(nid,"inertiedat",NF_FLOAT,3,dims3,nvarid)
+#endif
+      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 20,
+     &                        "Soil thermal inertia")
+      ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",15,
+     &                        "J.s-1/2.m-2.K-1")
+      ierr=NF_ENDDEF(nid) ! Leave NetCDF define mode
+      ! write variable
+#ifdef NC_DOUBLE
+      ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,ith)
+#else
+      ierr=NF_PUT_VAR_REAL (nid,nvarid,ith)
+#endif
+
+c-----------------------------------------------------------------------
+c  Ecriture aire et coefficients de passage cov. <-> contra. <--> naturel
+c-----------------------------------------------------------------------
+
+      ierr = NF_REDEF (nid)
+      dims2(1) = idim_rlonu
+      dims2(2) = idim_rlatu
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"cu",NF_DOUBLE,2,dims2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"cu",NF_FLOAT,2,dims2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 29,
+     .                       "Coefficient de passage pour U")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cu)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,cu)
+#endif
+c
+      ierr = NF_REDEF (nid)
+      dims2(1) = idim_rlonv
+      dims2(2) = idim_rlatv
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"cv",NF_DOUBLE,2,dims2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"cv",NF_FLOAT,2,dims2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 29,
+     .                       "Coefficient de passage pour V")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cv)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,cv)
+#endif
+c
+c Aire de chaque maille:
+c
+      ierr = NF_REDEF (nid)
+      dims2(1) = idim_rlonv
+      dims2(2) = idim_rlatu
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"aire",NF_DOUBLE,2,dims2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"aire",NF_FLOAT,2,dims2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22,
+     .                       "Aires de chaque maille")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aire)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,aire)
+#endif
+
+c-----------------------------------------------------------------------
+c  Ecriture du geopentiel au sol
+c-----------------------------------------------------------------------
+
+      ierr = NF_REDEF (nid)
+      dims2(1) = idim_rlonv
+      dims2(2) = idim_rlatu
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"phisinit",NF_DOUBLE,2,dims2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"phisinit",NF_FLOAT,2,dims2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19,
+     .                       "Geopotentiel au sol")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,phis)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,phis)
+#endif
+
+      PRINT*,'ini_archive: iim,jjm,llm,idayref',iim,jjm,llm,idayref
+      PRINT*,'ini_archive: rad,omeg,g,mugaz,kappa',
+     s rad,omeg,g,tab_cntrl_fi(8),kappa
+      PRINT*,'ini_archive: daysec,dtvr',daysec,dtvr
+
+      RETURN
+      END
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/iniphysiq_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/iniphysiq_mod.F90	(revision 1644)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/iniphysiq_mod.F90	(revision 1644)
@@ -0,0 +1,99 @@
+MODULE iniphysiq_mod
+
+CONTAINS
+
+subroutine iniphysiq(ii,jj,nlayer, &
+                     nbp, communicator, &
+                     punjours, pdayref,ptimestep, &
+                     rlatudyn,rlatvdyn,rlonudyn,rlonvdyn, &
+                     airedyn,cudyn,cvdyn, &
+                     prad,pg,pr,pcpp,iflag_phys)
+
+use control_mod, only: nday 
+use surf_heat_transp_mod, only: ini_surf_heat_transp
+use infotrac, only : nqtot ! number of advected tracers
+use planete_mod, only: ini_planete_mod
+USE comvert_mod, ONLY: ap,bp,preff
+use inifis_mod, only: inifis
+use ioipsl_getin_p_mod, only: getin_p
+
+use inigeomphy_mod, only: inigeomphy
+use geometry_mod, only: cell_area, & ! physics grid area (m2)
+                        longitude, & ! longitudes (rad)
+                        latitude ! latitudes (rad)
+! necessary to get klon_omp
+USE mod_phys_lmdz_para, ONLY: klon_omp ! number of columns (on local omp grid)
+USE dimphy, ONLY: init_dimphy
+
+implicit none
+include "dimensions.h"
+include "paramet.h"
+include "comgeom.h"
+include "iniprint.h"
+
+real,intent(in) :: prad ! radius of the planet (m)
+real,intent(in) :: pg ! gravitational acceleration (m/s2)
+real,intent(in) :: pr ! ! reduced gas constant R/mu
+real,intent(in) :: pcpp ! specific heat Cp
+real,intent(in) :: punjours ! length (in s) of a standard day
+integer,intent(in) :: nlayer ! number of atmospheric layers
+integer,intent(in) :: ii ! number of atmospheric coulumns along longitudes
+integer,intent(in) :: jj  ! number of atompsheric columns along latitudes
+integer,intent(in) :: nbp ! number of physics columns for this MPI process
+integer,intent(in) :: communicator ! MPI communicator
+real,intent(in) :: rlatudyn(jj+1) ! latitudes of the physics grid
+real,intent(in) :: rlatvdyn(jj) ! latitude boundaries of the physics grid
+real,intent(in) :: rlonvdyn(ii+1) ! longitudes of the physics grid
+real,intent(in) :: rlonudyn(ii+1) ! longitude boundaries of the physics grid
+real,intent(in) :: airedyn(ii+1,jj+1) ! area of the dynamics grid (m2)
+real,intent(in) :: cudyn((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)
+real,intent(in) :: cvdyn((ii+1)*jj) ! cv coeff. (v_covariant = cv * v)
+integer,intent(in) :: pdayref ! reference day of for the simulation
+real,intent(in) :: ptimestep !physics time step (s)
+integer,intent(in) :: iflag_phys ! type of physics to be called
+
+logical :: ok_slab_ocean
+
+  ! the common part for all planetary physics
+  !------------------------------------------
+  ! --> initialize physics distribution, global fields and geometry
+  ! (i.e. things in phy_common or dynphy_lonlat)
+  CALL inigeomphy(ii,jj,nlayer, &
+               nbp, communicator, &
+               rlatudyn,rlatvdyn, &
+               rlonudyn,rlonvdyn, &
+               airedyn,cudyn,cvdyn)
+
+  ! the distinct part for all planetary physics (ie. things in phystd)
+  !------------------------------------------
+
+!$OMP PARALLEL
+
+! copy some fundamental parameters to physics
+! and do some initializations
+
+! Initialize dimphy module
+call init_dimphy(klon_omp,nlayer)
+
+! copy over preff , ap() and bp() 
+call ini_planete_mod(nlayer,preff,ap,bp)
+
+! for slab ocean, copy over some arrays
+ok_slab_ocean=.false. ! default value
+call getin_p("ok_slab_ocean",ok_slab_ocean)
+if (ok_slab_ocean) then
+  call ini_surf_heat_transp(ip1jm,ip1jmp1,unsairez,fext,unsaire, &
+                            cu,cuvsurcv,cv,cvusurcu,aire,apoln,apols, &
+                            aireu,airev)
+endif
+
+call inifis(klon_omp,nlayer,nqtot,pdayref,punjours,nday,ptimestep, &
+            latitude,longitude,cell_area,prad,pg,pr,pcpp)
+
+
+!$OMP END PARALLEL
+
+end subroutine iniphysiq
+
+
+END MODULE iniphysiq_mod
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/interp_vert.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/interp_vert.F	(revision 1644)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/interp_vert.F	(revision 1644)
@@ -0,0 +1,70 @@
+c******************************************************
+      SUBROUTINE   interp_vert(varo,varn,lmo,lmn,apso,bpso,
+     &             aps,bps,ps,Nhoriz)
+c
+c interpolation lineaire pour passer
+c a une nouvelle discretisation verticale pour
+c les variables de GCM
+c Francois Forget (01/1995)
+c Modif pour coordonnees hybrides FF (03/2003)
+c**********************************************************
+
+      IMPLICIT NONE
+
+c   Declarations:
+c ==============
+c
+c  ARGUMENTS
+c  """""""""
+
+       integer lmo ! dimensions ancienne couches (input)
+       integer lmn ! dimensions nouvelle couches (input)
+
+       real apso(lmo),bpso(lmo)! anciennes coord hybrides midlayer (input)
+       real aps(lmn), bps(lmn)! nouvelles coord hybrides (midlayer) (input)
+
+       integer Nhoriz ! nombre de point horizontale (input)
+       real ps(nhoriz) !pression de surface (input)
+
+       real varo(Nhoriz,lmo) ! var dans l''ancienne grille (input)
+       real varn(Nhoriz,lmn) ! var dans la nouvelle grille (output)
+
+c Autres variables
+c """"""""""""""""
+       integer n, ln ,lo 
+       real coef
+       REAL sigmo(lmo) ! niveau sigma des variables dans les anciennes coord
+       REAL sigmn(lmn) ! niveau sigma des variables dans les nouvelles coord
+
+c run
+c ====
+
+      do n=1,Nhoriz
+        do ln=1,lmn
+            sigmn(ln)=aps(ln)/ps(n)+bps(ln)
+        end do
+        do lo=1,lmo
+            sigmo(lo)=apso(lo)/ps(n)+bpso(lo)
+        end do
+
+        do ln=1,lmn
+           if (sigmn(ln).ge.sigmo(1))then
+             varn(n,ln) =  varo(n,1)  
+           else if (sigmn(ln).le.sigmo(lmo)) then
+             varn(n,ln) =  varo(n,lmo)
+           else
+              do lo =1,lmo-1 
+                if ( (sigmn(ln).le.sigmo(lo)).and.
+     &             (sigmn(ln).gt.sigmo(lo+1)) )then
+                  coef = (sigmn(ln)-sigmo(lo))/(sigmo(lo+1)-sigmo(lo))
+                   varn(n,ln)=varo(n,lo) +coef*(varo(n,lo+1)-varo(n,lo))
+                end if
+              end do           
+           end if
+         end do
+
+      end do
+
+
+      return
+      end
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/lect_start_archive.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/lect_start_archive.F	(revision 1644)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/lect_start_archive.F	(revision 1644)
@@ -0,0 +1,1492 @@
+      SUBROUTINE lect_start_archive(ngrid,nlayer,
+     &     date,tsurf,tsoil,emis,q2,
+     &     t,ucov,vcov,ps,h,phisold_newgrid,
+     &     q,qsurf,surfith,nid,
+     &     rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
+
+!      USE surfdat_h
+      USE comsoil_h, ONLY: nsoilmx, layer, mlayer, volcapa, inertiedat
+      USE tracer_h, ONLY: igcm_co2_ice
+      USE infotrac, ONLY: tname, nqtot
+      USE slab_ice_h, ONLY: noceanmx
+!      USE control_mod
+! to use  'getin'
+      USE callkeys_mod, ONLY: ok_slab_ocean
+      USE comvert_mod, ONLY: ap,bp,aps,bps,preff
+      USE comconst_mod, ONLY: kappa,g,pi
+
+c=======================================================================
+c
+c
+c   Auteur:    05/1997 , 12/2003 : coord hybride  FF
+c   ------
+c
+c
+c   Objet:     Lecture des variables d'un fichier "start_archive"
+c              Plus besoin de régler ancienne valeurs grace
+c              a l'allocation dynamique de memoire (Yann Wanherdrick)
+c
+c
+c
+c=======================================================================
+
+      implicit none
+
+#include "dimensions.h"
+!#include "dimphys.h"
+!#include "planete.h"
+#include "paramet.h"
+#include "comgeom2.h"
+!#include "control.h"
+#include "netcdf.inc"
+!#include"advtrac.h"
+c=======================================================================
+c   Declarations
+c=======================================================================
+
+      INTEGER,INTENT(IN) :: ngrid, nlayer
+
+c Old variables dimensions (from file)
+c------------------------------------
+      INTEGER   imold,jmold,lmold,nsoilold,nqold
+
+
+c Variables pour les lectures des fichiers "ini" 
+c--------------------------------------------------
+!      INTEGER sizei,
+      integer timelen,dimid
+!      INTEGER length
+!      parameter (length = 100)
+      INTEGER tab0
+      INTEGER isoil,iq,iqmax
+      CHARACTER*2   str2
+
+!      REAL dimfirst(4) ! tableau contenant les 1ers elements des dimensions
+
+!      REAL dimlast(4) ! tableau contenant les derniers elements des dimensions
+
+!      REAL dimcycl(4) ! tableau contenant les periodes des dimensions
+!      CHARACTER*120 dimsource
+!      CHARACTER*16 dimname
+!      CHARACTER*80 dimtitle
+!      CHARACTER*40 dimunits
+!      INTEGER   dimtype
+
+!      INTEGER dimord(4)  ! tableau contenant l''ordre
+!      data dimord /1,2,3,4/ ! de sortie des dimensions
+
+!      INTEGER vardim(4)
+      REAL date
+      INTEGER   memo
+!      character (len=50) :: tmpname
+
+c Variable histoire 
+c------------------
+      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
+      REAL h(iip1,jjp1,llm),ps(iip1,jjp1)
+      REAL q(iip1,jjp1,llm,nqtot),qtot(iip1,jjp1,llm)
+
+c autre variables dynamique nouvelle grille
+c------------------------------------------
+
+c!-*-
+!      integer klatdat,klongdat
+!      PARAMETER (klatdat=180,klongdat=360)
+
+c Physique sur grille scalaire 
+c----------------------------
+
+c variable physique
+c------------------
+      REAL tsurf(ngrid) ! surface temperature
+      REAL tsoil(ngrid,nsoilmx) ! soil temperature
+      REAL co2ice(ngrid) ! CO2 ice layer
+      REAL emis(ngrid)
+      REAL q2(ngrid,llm+1),qsurf(ngrid,nqtot)
+      REAL tslab(ngrid,noceanmx)
+      REAL rnat(ngrid),pctsrf_sic(ngrid)
+      REAL tsea_ice(ngrid),sea_ice(ngrid)
+c     REAL phisfi(ngrid)
+
+      INTEGER i,j,l
+      INTEGER nid,nvarid
+c     REAL year_day,periheli,aphelie,peri_day
+c     REAL obliquit,z0,emin_turb,lmixmin
+c     REAL emissiv,emisice(2),albedice(2)
+c     REAL iceradius(2) , dtemisice(2)
+
+      integer ierr
+      integer, dimension(4) :: start,count
+
+c Variable nouvelle grille naturelle au point scalaire
+c------------------------------------------------------
+      real us(iip1,jjp1,llm),vs(iip1,jjp1,llm)
+      REAL phisold_newgrid(iip1,jjp1)
+      REAL t(iip1,jjp1,llm)
+      real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx)
+      real inertiedatS(iip1,jjp1,nsoilmx)
+      real co2iceS(iip1,jjp1)
+      real emisS(iip1,jjp1)
+      REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot)
+      real tslabS(iip1,jjp1,noceanmx)
+      real pctsrf_sicS(iip1,jjp1),tsea_iceS(iip1,jjp1)
+      real rnatS(iip1,jjp1), sea_iceS(iip1,jjp1)
+
+      real ptotal, co2icetotal
+
+c Var intermediaires : vent naturel, mais pas coord scalaire
+c-----------------------------------------------------------
+      real vnat(iip1,jjm,llm),unat(iip1,jjp1,llm)
+
+
+c Variable de l'ancienne grille 
+c---------------------------------------------------------
+
+      real, dimension(:), allocatable :: timelist
+      real, dimension(:), allocatable :: rlonuold, rlatvold
+      real, dimension(:), allocatable :: rlonvold, rlatuold
+      real, dimension(:), allocatable :: apsold,bpsold
+      real, dimension(:), allocatable :: mlayerold
+      real, dimension(:,:,:), allocatable :: uold,vold,told,q2old
+      real, dimension(:,:,:), allocatable :: tsoilold,qsurfold
+      real, dimension(:,:,:),allocatable :: tsoiloldnew
+! tsoiloldnew: old soil values, but along new subterranean grid
+      real, dimension(:,:,:), allocatable :: inertiedatold
+! inertiedatoldnew: old inertia values, but along new subterranean grid
+      real, dimension(:,:,:), allocatable :: inertiedatoldnew
+      real, dimension(:,:), allocatable :: psold,phisold
+      real, dimension(:,:), allocatable :: co2iceold
+      real, dimension(:,:), allocatable :: tsurfold
+      real, dimension(:,:), allocatable :: emisold
+      real, dimension(:,:,:,:), allocatable :: qold
+      real, dimension(:,:,:), allocatable :: tslabold
+      real, dimension(:,:), allocatable :: rnatold,pctsrf_sicold
+      real, dimension(:,:), allocatable :: tsea_iceold,sea_iceold
+
+
+      real tab_cntrl(100)
+
+      real ptotalold, co2icetotalold
+
+      logical :: olddepthdef=.false. ! flag
+! olddepthdef=.true. if soil depths are in 'old' (unspecified) format
+      logical :: depthinterpol=.false. ! flag
+! depthinterpol=.true. if interpolation will be requiered
+      logical :: therminertia_3D=.true. ! flag
+! therminertia_3D=.true. if thermal inertia is 3D and read from datafile
+c Variable intermediaires iutilise pour l'extrapolation verticale 
+c----------------------------------------------------------------
+      real, dimension(:,:,:), allocatable :: var,varp1 
+      real, dimension(:), allocatable :: oldgrid, oldval
+      real, dimension(:), allocatable :: newval
+
+      real surfith(iip1,jjp1) ! surface thermal inertia
+      ! surface thermal inertia at old horizontal grid resolution
+      real, dimension(:,:), allocatable :: surfithold 
+
+      character(len=30) :: txt ! to store some text
+
+c=======================================================================
+
+! 0. Preliminary stuff
+
+
+
+!-----------------------------------------------------------------------
+! 1. Read data dimensions (i.e. size and length)
+!-----------------------------------------------------------------------
+
+! 1.2 Read the various dimension lengths of data in file 
+
+      ierr= NF_INQ_DIMID(nid,"Time",dimid)
+      if (ierr.ne.NF_NOERR) then
+         ierr= NF_INQ_DIMID(nid,"temps",dimid)
+      endif
+      ierr= NF_INQ_DIMLEN(nid,dimid,timelen)
+
+      ierr= NF_INQ_DIMID(nid,"latitude",dimid)
+      if (ierr.ne.NF_NOERR) then
+         ierr= NF_INQ_DIMID(nid,"rlatu",dimid)
+      endif
+      ierr= NF_INQ_DIMLEN(nid,dimid,jmold)
+      jmold=jmold-1
+
+      ierr= NF_INQ_DIMID(nid,"longitude",dimid)
+      if (ierr.ne.NF_NOERR) then
+         ierr= NF_INQ_DIMID(nid,"rlonv",dimid)
+      endif
+      ierr= NF_INQ_DIMLEN(nid,dimid,imold)
+      imold=imold-1
+
+      ierr= NF_INQ_DIMID(nid,"altitude",dimid)
+      if (ierr.ne.NF_NOERR) then
+         ierr= NF_INQ_DIMID(nid,"sig_s",dimid)
+      endif
+      ierr= NF_INQ_DIMLEN(nid,dimid,lmold)
+
+      nqold=0
+      do
+         write(str2,'(i2.2)') nqold+1
+         ierr= NF_INQ_VARID(nid,'q'//str2,dimid)
+!        write(*,*) 'q'//str2
+         if (ierr.eq.NF_NOERR) then
+            nqold=nqold+1
+         else
+            exit
+         endif
+      enddo
+
+! 1.2.1 find out the # of subsurface_layers
+      nsoilold=0 !dummy initialisation
+      ierr=NF_INQ_DIMID(nid,"subsurface_layers",dimid)
+      if (ierr.eq.NF_NOERR) then
+        ierr=NF_INQ_DIMLEN(nid,dimid,nsoilold)
+	if (ierr.ne.NF_NOERR) then
+	 write(*,*)'lec_start_archive: ',
+     &              'Failed reading subsurface_layers length'
+	endif
+      else
+        write(*,*)"lec_start_archive: did not find subsurface_layers"
+      endif
+
+      if (nsoilold.eq.0) then ! 'old' archive format;
+      ! must use Tg//str2 fields to compute nsoilold
+      write(*,*)"lec_start_archive: building nsoilold from Tg fields"
+        do
+	 write(str2,'(i2.2)') nsoilold+1
+	 ierr=NF_INQ_VARID(nid,'Tg'//str2,dimid)
+	 if (ierr.eq.NF_NOERR) then
+	  nsoilold=nsoilold+1
+	 else
+	  exit
+	 endif
+	enddo
+      endif
+
+
+      if (nsoilold.ne.nsoilmx) then ! interpolation will be required
+        depthinterpol=.true.
+      endif
+
+! 1.3 Report dimensions
+      
+      write(*,*) "Start_archive dimensions:"
+      write(*,*) "longitude: ",imold
+      write(*,*) "latitude: ",jmold
+      write(*,*) "altitude: ",lmold
+      write(*,*) "tracers: ",nqold
+      write(*,*) "subsurface_layers: ",nsoilold
+      if (depthinterpol) then
+      write(*,*) " => Warning, nsoilmx= ",nsoilmx
+      write(*,*) '    which implies that you want subterranean interpola
+     &tion.'
+      write(*,*) '  Otherwise, set nsoilmx -in dimphys.h- to: ',nsoilold
+      endif
+      write(*,*) "time lenght: ",timelen
+      write(*,*) 
+
+!-----------------------------------------------------------------------
+! 2. Allocate arrays to store datasets
+!-----------------------------------------------------------------------
+
+      allocate(timelist(timelen))
+      allocate(rlonuold(imold+1), rlatvold(jmold))
+      allocate(rlonvold(imold+1), rlatuold(jmold+1))
+      allocate (apsold(lmold),bpsold(lmold))
+      allocate(uold(imold+1,jmold+1,lmold))
+      allocate(vold(imold+1,jmold+1,lmold))
+      allocate(told(imold+1,jmold+1,lmold))
+      allocate(psold(imold+1,jmold+1))
+      allocate(phisold(imold+1,jmold+1))
+      allocate(qold(imold+1,jmold+1,lmold,nqtot))
+      allocate(co2iceold(imold+1,jmold+1))
+      allocate(tsurfold(imold+1,jmold+1))
+      allocate(emisold(imold+1,jmold+1))
+      allocate(q2old(imold+1,jmold+1,lmold+1))
+!      allocate(tsoilold(imold+1,jmold+1,nsoilmx))
+      allocate(tsoilold(imold+1,jmold+1,nsoilold))
+      allocate(tsoiloldnew(imold+1,jmold+1,nsoilmx))
+      allocate(inertiedatold(imold+1,jmold+1,nsoilold)) ! soil thermal inertia
+      allocate(inertiedatoldnew(imold+1,jmold+1,nsoilmx))
+      ! surface thermal inertia at old horizontal grid resolution
+      allocate(surfithold(imold+1,jmold+1))
+      allocate(mlayerold(nsoilold))
+      allocate(qsurfold(imold+1,jmold+1,nqtot))
+      allocate(tslabold(imold+1,jmold+1,noceanmx))
+      allocate(rnatold(imold+1,jmold+1))
+      allocate(pctsrf_sicold(imold+1,jmold+1))
+      allocate(tsea_iceold(imold+1,jmold+1))
+      allocate(sea_iceold(imold+1,jmold+1))
+
+      allocate(var (imold+1,jmold+1,llm))
+      allocate(varp1 (imold+1,jmold+1,llm+1))
+
+      write(*,*) 'q2',ngrid,llm+1
+      write(*,*) 'q2S',iip1,jjp1,llm+1
+      write(*,*) 'q2old',imold+1,jmold+1,lmold+1
+
+!-----------------------------------------------------------------------
+! 3. Read time-independent data
+!-----------------------------------------------------------------------
+
+C-----------------------------------------------------------------------
+c 3.1. Lecture du tableau des parametres du run 
+c     (pour  la lecture ulterieure de "ptotalold" et "co2icetotalold")
+c-----------------------------------------------------------------------
+c
+      ierr = NF_INQ_VARID (nid, "controle", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "Lect_start_archive: champ <controle> est absent"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
+#else
+      ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echoue pour <controle>"
+         CALL abort
+      ENDIF
+c
+      tab0 = 50
+
+c-----------------------------------------------------------------------
+c 3.2 Lecture des longitudes et latitudes
+c-----------------------------------------------------------------------
+c
+      ierr = NF_INQ_VARID (nid, "rlonv", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <rlonv> est absent"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonvold)
+#else
+      ierr = NF_GET_VAR_REAL(nid, nvarid, rlonvold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <rlonv>"
+         CALL abort
+      ENDIF
+c
+      ierr = NF_INQ_VARID (nid, "rlatu", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <rlatu> est absent"
+         CALL abort
+      ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatuold)
+#else
+      ierr = NF_GET_VAR_REAL(nid, nvarid, rlatuold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <rlatu>"
+         CALL abort
+      ENDIF
+c
+      ierr = NF_INQ_VARID (nid, "rlonu", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <rlonu> est absent"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonuold)
+#else
+      ierr = NF_GET_VAR_REAL(nid, nvarid, rlonuold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <rlonu>"
+         CALL abort
+      ENDIF
+c
+      ierr = NF_INQ_VARID (nid, "rlatv", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <rlatv> est absent"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatvold)
+#else
+      ierr = NF_GET_VAR_REAL(nid, nvarid, rlatvold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <rlatv>"
+         CALL abort
+      ENDIF
+c
+
+c-----------------------------------------------------------------------
+c 3.3. Lecture des niveaux verticaux
+c-----------------------------------------------------------------------
+c
+      ierr = NF_INQ_VARID (nid, "aps", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <aps> est absent"
+         apsold=0
+         PRINT*, "<aps> set to 0"
+      ELSE
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, apsold)
+#else
+         ierr = NF_GET_VAR_REAL(nid, nvarid, apsold)
+#endif
+         IF (ierr .NE. NF_NOERR) THEN
+            PRINT*, "lect_start_archive: Lecture echouee pour <aps>"
+         ENDIF
+      ENDIF
+c
+      ierr = NF_INQ_VARID (nid, "bps", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <bps> est absent"
+         PRINT*, "It must be an old start_archive, lets look for sig_s"
+         ierr = NF_INQ_VARID (nid, "sig_s", nvarid)
+         IF (ierr .NE. NF_NOERR) THEN
+            PRINT*, "Nothing to do..."
+            CALL abort
+         ENDIF
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bpsold)
+#else
+      ierr = NF_GET_VAR_REAL(nid, nvarid, bpsold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <bps>"
+         CALL abort
+      END IF
+
+c-----------------------------------------------------------------------
+c 3.4 Read Soil layers depths
+c-----------------------------------------------------------------------
+     
+      ierr=NF_INQ_VARID(nid,"soildepth",nvarid)
+      if (ierr.ne.NF_NOERR) then
+       write(*,*)'lect_start_archive: Could not find <soildepth>'
+       write(*,*)' => Assuming this is an archive in old format'
+       olddepthdef=.true.
+       depthinterpol=.true.
+       ! this is how soil depth was defined in ye old days
+	do isoil=1,nsoilold
+	  mlayerold(isoil)=sqrt(887.75/3.14)*((2.**(isoil-0.5))-1.)
+	enddo
+      else
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayerold)
+#else
+        ierr = NF_GET_VAR_REAL(nid,nvarid,mlayerold)
+#endif
+       if (ierr .NE. NF_NOERR) then
+         PRINT*, "lect_start_archive: Failed reading <soildepth>"
+         CALL abort
+       endif
+
+      endif !of if(ierr.ne.NF_NOERR)
+
+      ! Read (or build) mlayer()
+      if (depthinterpol) then
+       ! Build (default) new soil depths (mlayer(:) is in comsoil.h),
+       ! as in soil_settings.F
+       write(*,*)' => Building default soil depths'
+       do isoil=0,nsoilmx-1
+         mlayer(isoil)=2.e-4*(2.**(isoil-0.5))
+       enddo
+       write(*,*)' => mlayer: ',mlayer
+       ! Also build (default) new soil interlayer depth layer(:)
+       do isoil=1,nsoilmx
+         layer(isoil)=sqrt(mlayer(0)*mlayer(1))*
+     &                      ((mlayer(1)/mlayer(0))**(isoil-1))
+       enddo
+       write(*,*)' =>  layer: ',layer
+      else ! read mlayer() from file
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayer)
+#else
+        ierr = NF_GET_VAR_REAL(nid,nvarid,mlayer)
+#endif
+       if (ierr .NE. NF_NOERR) then
+         PRINT*, "lect_start_archive: Failed reading <soildepth>"
+         CALL abort
+       endif
+      endif ! of if (depthinterpol)
+
+c-----------------------------------------------------------------------
+c 3.5 Read Soil thermal inertia
+c-----------------------------------------------------------------------
+
+      ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
+      if (ierr.ne.NF_NOERR) then
+       write(*,*)'lect_start_archive: Could not find <inertiedat>'
+       write(*,*)' => Assuming this is an archive in old format'
+       therminertia_3D=.false.
+       write(*,*)' => Thermal inertia will be read from reference file'
+       volcapa=1.e6
+       write(*,*)'    and soil volumetric heat capacity is set to ',
+     &           volcapa
+      else
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VAR_DOUBLE(nid,nvarid,inertiedatold)
+#else
+        ierr = NF_GET_VAR_REAL(nid,nvarid,inertiedatold)
+#endif
+       if (ierr .NE. NF_NOERR) then
+         PRINT*, "lect_start_archive: Failed reading <inertiedat>"
+         CALL abort
+       endif
+      endif
+
+c-----------------------------------------------------------------------
+c 3.6 Lecture geopotentiel au sol
+c-----------------------------------------------------------------------
+c
+      ierr = NF_INQ_VARID (nid, "phisinit", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <phisinit> est absent"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisold)
+#else
+      ierr = NF_GET_VAR_REAL(nid, nvarid, phisold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <phisinit>"
+         CALL abort
+      ENDIF
+
+C-----------------------------------------------------------------------
+c   lecture de "ptotalold" et "co2icetotalold"
+c-----------------------------------------------------------------------
+      ptotalold = tab_cntrl(tab0+49)
+      co2icetotalold = tab_cntrl(tab0+50)
+ 
+c-----------------------------------------------------------------------
+c 4. Lecture du temps et choix
+c-----------------------------------------------------------------------
+ 
+c  lecture du temps
+c
+      ierr = NF_INQ_DIMID (nid, "Time", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         ierr = NF_INQ_DIMID (nid, "temps", nvarid)
+         IF (ierr .NE. NF_NOERR) THEN
+            PRINT*, "lect_start_archive: Le champ <Time> est absent"
+            CALL abort
+         endif
+      ENDIF
+
+      ierr = NF_INQ_VARID (nid, "Time", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         ierr = NF_INQ_VARID (nid, "temps", nvarid)
+      endif 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, timelist)
+#else
+      ierr = NF_GET_VAR_REAL(nid, nvarid, timelist)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <Time>"
+         CALL abort
+      ENDIF
+c
+      write(*,*)
+      write(*,*)
+      write(*,*) 'Differentes dates des etats initiaux stockes:'
+      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+      pi=2.*ASIN(1.)
+      do i=1,timelen
+c       call solarlong(timelist(i),sollong(i))
+c       sollong(i) = sollong(i)*180./pi
+        write(*,*) 'etat initial au jour martien' ,int(timelist(i))
+c       write(*,6) nint(timelist(i)),nint(mod(timelist(i),669)),
+c    .    sollong(i)
+      end do
+
+   6  FORMAT(i7,i7,f9.3)
+ 
+      write(*,*)
+      write(*,*) 'Choix de la date'
+ 123  read(*,*,iostat=ierr) date
+      if(ierr.ne.0) goto 123
+      memo = 0
+      do i=1,timelen
+        if (date.eq.int(timelist(i))) then
+            memo = i
+        endif
+      end do
+ 
+      if (memo.eq.0) then
+        write(*,*)
+        write(*,*)
+        write(*,*) 'He alors... Y sait pas lire !?!'
+        write(*,*)
+        write(*,*) 'Differentes dates des etats initiaux stockes:'
+        write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+        do i=1,timelen
+          write(*,*) 'etat initial au jour martien' ,nint(timelist(i))
+c         write(*,6) nint(timelist(i)),nint(mod(timelist(i),669))
+        end do
+        goto 123
+      endif
+
+!-----------------------------------------------------------------------
+! 5. Read (time-dependent) data from datafile
+!-----------------------------------------------------------------------
+
+
+c-----------------------------------------------------------------------
+c 5.1 Lecture des champs 2D (co2ice, emis,ps,tsurf,Tg[10], qsurf)
+c-----------------------------------------------------------------------
+ 
+      start=(/1,1,memo,0/)
+      count=(/imold+1,jmold+1,1,0/)
+       
+! CO2ice is now in qsurf(igcm_co2_ice) ...
+!      ierr = NF_INQ_VARID (nid, "co2ice", nvarid)
+!      IF (ierr .NE. NF_NOERR) THEN
+!         PRINT*, "lect_start_archive: Le champ <co2ice> est absent"
+!         CALL abort
+!      ENDIF
+!#ifdef NC_DOUBLE
+!      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,co2iceold)
+!#else
+!      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,co2iceold)
+!#endif
+!      IF (ierr .NE. NF_NOERR) THEN
+!         PRINT*, "lect_start_archive: Lecture echouee pour <co2ice>"
+!         PRINT*, NF_STRERROR(ierr)
+!         CALL abort
+!      ENDIF
+c
+      ierr = NF_INQ_VARID (nid, "emis", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <emis> est absent"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,emisold)
+#else
+      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,emisold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <emis>"
+         CALL abort
+      ENDIF
+c
+      ierr = NF_INQ_VARID (nid, "ps", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <ps> est absent"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,psold)
+#else
+      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,psold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <ps>"
+         CALL abort
+      ENDIF
+c
+      ierr = NF_INQ_VARID (nid, "tsurf", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <tsurf> est absent"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,tsurfold)
+#else
+      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,tsurfold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <tsurf>"
+         CALL abort
+      ENDIF
+c
+      ierr = NF_INQ_VARID (nid, "q2surf", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <q2surf> est absent"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,q2old)
+#else
+      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,q2old)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <q2surf>"
+         CALL abort
+      ENDIF
+c
+cc Slab ocean
+      if(ok_slab_ocean) then
+      start=(/1,1,1,memo/)
+      count=(/imold+1,jmold+1,noceanmx,1/)
+
+       ierr=NF_INQ_VARID(nid,"tslab",nvarid)
+       IF (ierr.ne.NF_NOERR) then
+          PRINT*,"lect_start_archive: Cannot find <tslab>"
+       ENDIF
+#ifdef NC_DOUBLE
+      ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,tslabold)
+#else
+      ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,tslabold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <tslab>"
+      ENDIF
+
+
+c
+      start=(/1,1,memo,0/)
+      count=(/imold+1,jmold+1,1,0/)
+
+      ierr = NF_INQ_VARID (nid, "rnat", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <rnat> est absent"
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,rnatold)
+#else
+      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,rnatold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <rnat>"
+      ENDIF
+c
+      ierr = NF_INQ_VARID (nid, "pctsrf_sic", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <pctsrf_sic> est absent"
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,pctsrf_sicold)
+#else
+      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,pctsrf_sicold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <pctsrf_sic>"
+      ENDIF
+c
+      ierr = NF_INQ_VARID (nid, "tsea_ice", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <tsea_ice> est absent"
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,tsea_iceold)
+#else
+      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,tsea_iceold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <tsea_ice>"
+      ENDIF
+c
+      ierr = NF_INQ_VARID (nid, "sea_ice", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <sea_ice> est absent"
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,sea_iceold)
+#else
+      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,sea_iceold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <sea_ice>"
+      ENDIF
+  
+      ENDIF! ok_slab_ocean
+c
+      write(*,*)"lect_start_archive: rlonuold:"
+     &           ,rlonuold," rlatvold:",rlatvold
+      write(*,*)
+
+! Surface tracers:      
+      do iq=1,nqtot
+        ! initialize all surface tracers to zero
+        call initial0((jmold+1)*(imold+1), qsurfold(1,1,iq))
+      enddo
+
+
+!      print*,'tname=',tname
+!      print*,'nid',nid
+!      print*,'nvarid',nvarid
+!      stop
+
+      DO iq=1,nqtot
+          txt=trim(tname(iq))//"_surf"
+          if (txt.eq."h2o_vap") then
+            ! There is no surface tracer for h2o_vap;
+            ! "h2o_ice" should be loaded instead
+            txt="h2o_ice_surf"
+            write(*,*) 'lect_start_archive: loading surface tracer',
+     &                     ' h2o_ice instead of h2o_vap'
+          endif
+
+        
+        write(*,*) "lect_start_archive: loading tracer ",trim(txt)
+        ierr = NF_INQ_VARID (nid,txt,nvarid)
+        IF (ierr .NE. NF_NOERR) THEN
+          PRINT*, "lect_start_archive: ",
+     &              " Tracer <",trim(txt),"> not found"
+
+!          print*,'RDW has added hack to let me continue...'
+!          CALL abort
+        ENDIF
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
+     &          qsurfold(1,1,iq))
+#else
+        ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
+     &          qsurfold(1,1,iq))
+#endif
+        IF (ierr .NE. NF_NOERR) THEN
+          PRINT*, "lect_start_archive: ",
+     &             " Failed loading <",trim(txt),">"
+          write (*,*) trim(txt),'    is set to 0'
+!          call initial0((jmold+1)*(imold+1), qsurfold(1,1,iq))
+        ENDIF
+
+      ENDDO ! of DO iq=1,nqtot
+
+
+!-----------------------------------------------------------------------
+! 5.2 Read 3D subterranean fields
+!-----------------------------------------------------------------------
+
+      start=(/1,1,1,memo/)
+      count=(/imold+1,jmold+1,nsoilold,1/)
+!
+! Read soil temperatures
+!
+      if (olddepthdef) then ! tsoil stored using the 'old format'
+         start=(/1,1,memo,0/)
+         count=(/imold+1,jmold+1,1,0/) ! because the "Tg" are 2D datasets
+       do isoil=1,nsoilold
+         write(str2,'(i2.2)') isoil
+c
+         ierr = NF_INQ_VARID (nid, "Tg"//str2, nvarid)
+         IF (ierr .NE. NF_NOERR) THEN
+            PRINT*, "lect_start_archive: ",
+     &              "Field <","Tg"//str2,"> not found"
+            CALL abort
+         ENDIF
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
+     &          tsoilold(1,1,isoil))
+#else
+         ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
+     &          tsoilold(1,1,isoil))
+#endif
+         IF (ierr .NE. NF_NOERR) THEN
+            PRINT*, "lect_start_archive: ",
+     &            "Failed reading <","Tg"//str2,">"
+            CALL abort
+         ENDIF
+c
+       enddo ! of do isoil=1,nsoilold
+      
+      ! reset 'start' and 'count' to "3D" behaviour
+      start=(/1,1,1,memo/)
+      count=(/imold+1,jmold+1,nsoilold,1/)
+      
+      else
+       write(*,*) "lect_start_archive: loading tsoil "
+       ierr=NF_INQ_VARID(nid,"tsoil",nvarid)
+       if (ierr.ne.NF_NOERR) then
+        write(*,*)"lect_start_archive: Cannot find <tsoil>"
+	call abort
+       else
+#ifdef NC_DOUBLE
+      ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,tsoilold)
+#else
+      ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,tsoilold)
+#endif
+       endif ! of if (ierr.ne.NF_NOERR)
+       
+      endif ! of if (olddepthdef)
+
+!
+! Read soil thermal inertias
+!
+!      if (.not.olddepthdef) then ! no thermal inertia data in "old" archives
+!       ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
+!       if (ierr.ne.NF_NOERR) then
+!        write(*,*)"lect_start_archive: Cannot find <inertiedat>"
+!	call abort
+!       else
+!#ifdef NC_DOUBLE
+!      ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,inertiedatold)
+!#else
+!      ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,inertiedatold)
+!#endif
+!       endif ! of if (ierr.ne.NF_NOERR)
+!      endif
+
+c-----------------------------------------------------------------------
+c 5.3	Lecture des champs 3D (t,u,v, q2atm,q)
+c-----------------------------------------------------------------------
+
+      start=(/1,1,1,memo/)
+      count=(/imold+1,jmold+1,lmold,1/)
+
+c
+      ierr = NF_INQ_VARID (nid,"temp", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <temp> est absent"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid, nvarid, start, count, told)
+#else
+      ierr = NF_GET_VARA_REAL(nid, nvarid, start, count, told)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <temp>"
+         CALL abort
+      ENDIF
+c
+      ierr = NF_INQ_VARID (nid,"u", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <u> est absent"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,uold)
+#else
+      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,uold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <u>"
+         CALL abort
+      ENDIF
+c
+      ierr = NF_INQ_VARID (nid,"v", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <v> est absent"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,vold)
+#else
+      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,vold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <v>"
+         CALL abort
+      ENDIF
+c
+      ierr = NF_INQ_VARID (nid,"q2atm", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <q2atm> est absent"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,q2old(1,1,2))
+#else
+      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,q2old(1,1,2))
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <q2atm>"
+         CALL abort
+      ENDIF
+c
+
+! Tracers:      
+      do iq=1,nqtot
+         call initial0((jmold+1)*(imold+1)*lmold,qold(1,1,1,iq) )
+      enddo
+
+      DO iq=1,nqtot
+        txt=tname(iq)
+        write(*,*)"lect_start_archive: loading tracer ",trim(txt)
+        ierr = NF_INQ_VARID (nid,txt,nvarid)
+        IF (ierr .NE. NF_NOERR) THEN
+            PRINT*, "lect_start_archive: ",
+     &              " Tracer <",trim(txt),"> not found"
+!            CALL abort
+        ENDIF
+#ifdef NC_DOUBLE
+        ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start,count,qold(1,1,1,iq))
+#else
+        ierr=NF_GET_VARA_REAL(nid,nvarid,start,count,qold(1,1,1,iq))
+#endif
+        IF (ierr .NE. NF_NOERR) THEN
+          PRINT*, "lect_start_archive: ",
+     &             "  Failed loading <",trim(txt),">"
+          write (*,*) trim(txt),'      set to 1.E-30'
+          do l=1,lmold
+            do j=1,jmold+1
+              do i=1,imold+1
+                 qold(i,j,l,iq)=1.e-30
+              end do
+            end do
+          end do
+        ENDIF
+
+      ENDDO ! of DO iq=1,nqtot
+
+
+!=======================================================================
+! 6. Interpolation from old grid to new grid
+!=======================================================================
+
+c=======================================================================
+c   INTERPOLATION DANS LA NOUVELLE GRILLE et initialisation des variables
+c=======================================================================
+c  Interpolation horizontale puis passage dans la grille physique pour 
+c  les variables physique 
+c  Interpolation verticale puis horizontale pour chaque variable 3D
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c 6.1	Variable 2d :
+c-----------------------------------------------------------------------
+c Relief 
+      call interp_horiz (phisold,phisold_newgrid,imold,jmold,iim,jjm,1,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+
+! CO2 ice is now in qsurf(igcm_co2_ice)
+!      call interp_horiz (co2iceold,co2ices,imold,jmold,iim,jjm,1,
+!     &                   rlonuold,rlatvold,rlonu,rlatv)
+
+c Temperature de surface
+      call interp_horiz (tsurfold,tsurfs,imold,jmold,iim,jjm,1,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tsurfs,tsurf)
+c     write(44,*) 'tsurf', tsurf
+
+c Temperature du sous-sol
+!      call interp_horiz(tsoilold,tsoils,
+!     &                  imold,jmold,iim,jjm,nsoilmx,
+!     &                   rlonuold,rlatvold,rlonu,rlatv)
+!      call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,tsoils,tsoil)
+c     write(45,*) 'tsoil',tsoil
+
+c Emissivite de la surface
+      call interp_horiz (emisold,emiss,imold,jmold,iim,jjm,1,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,emiss,emis)
+c     write(46,*) 'emis',emis
+
+
+
+c-----------------------------------------------------------------------
+c 6.1.2	Traitement special de la pression au sol :
+c-----------------------------------------------------------------------
+
+c  Extrapolation la pression dans la nouvelle grille
+      call interp_horiz(psold,ps,imold,jmold,iim,jjm,1,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+
+c-----------------------------------------------------------------------
+c	On assure la conservation de la masse de l'atmosphere + calottes
+c-----------------------------------------------------------------------
+
+      ptotal =  0.
+      DO j=1,jjp1
+         DO i=1,iim
+            ptotal=ptotal+ps(i,j)*aire(i,j)/g
+         END DO
+      END DO
+      co2icetotal = 0.
+      if (igcm_co2_ice.ne.0) then
+        ! recast surface CO2 ice on new grid
+        call interp_horiz(qsurfold(1,1,igcm_co2_ice),
+     &                  qsurfs(1,1,igcm_co2_ice),
+     &                  imold,jmold,iim,jjm,1,
+     &                  rlonuold,rlatvold,rlonu,rlatv)
+       DO j=1,jjp1
+         DO i=1,iim
+           !co2icetotal = co2icetotal + co2iceS(i,j)*aire(i,j)
+           co2icetotal=co2icetotal+qsurfS(i,j,igcm_co2_ice)*aire(i,j)
+         END DO
+       END DO
+      else
+        write(*,*) "Warning: No co2_ice surface tracer"
+      endif
+
+      write(*,*)
+      write(*,*)'Ancienne grille: masse de l atm :',ptotalold
+      write(*,*)'Nouvelle grille: masse de l atm :',ptotal
+      write (*,*) 'Ratio new atm./ old atm =', ptotal/ptotalold 
+      write(*,*)
+      write(*,*)'Ancienne grille: masse de la glace CO2:',co2icetotalold
+      write(*,*)'Nouvelle grille: masse de la glace CO2:',co2icetotal
+      if (co2icetotalold.ne.0.) then
+      write(*,*)'Ratio new ice./old ice =',co2icetotal/co2icetotalold
+      endif
+      write(*,*)
+
+
+      DO j=1,jjp1
+         DO i=1,iip1
+            ps(i,j)=ps(i,j) * ptotalold/ptotal
+         END DO
+      END DO
+
+      if ( co2icetotalold.gt.0.) then 
+!         DO j=1,jjp1
+!            DO i=1,iip1
+!               co2iceS(i,j)=co2iceS(i,j) * co2icetotalold/co2icetotal
+!            END DO
+!         END DO
+        write(*,*) "Not enforcing conservation of surface CO2 ice"
+        write(*,*) " (should be OK when CO2 ice is a tracer)"
+      end if
+
+c-----------------------------------------------------------------------
+c 6.2 Subterranean 3d variables:
+c-----------------------------------------------------------------------
+
+c-----------------------------------------------------------------------
+c 6.2.1 Thermal Inertia
+c       Note: recall that inertiedat is a common in "comsoil.h"
+c-----------------------------------------------------------------------
+
+      ! depth-wise interpolation, if required
+      if (depthinterpol.and.(.not.olddepthdef)) then
+        allocate(oldval(nsoilold))
+	allocate(newval(nsoilmx))
+        write(*,*)'lect_start_archive: WARNING: vertical interpolation o
+     &f soil thermal inertia; might be wiser to reset it.'
+        write(*,*)
+       
+        do i=1,imold+1
+         do j=1,jmold+1
+	   !copy old values
+	   oldval(1:nsoilold)=inertiedatold(i,j,1:nsoilold)
+	   !interpolate
+	   call interp_line(mlayerold,oldval,nsoilold,
+     &                     mlayer,newval,nsoilmx)
+           !copy interpolated values
+           inertiedatoldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
+	 enddo
+        enddo
+        ! cleanup
+	deallocate(oldval)
+	deallocate(newval)
+      endif !of if (depthinterpol)
+
+      if (therminertia_3D) then
+        ! We have inertiedatold
+       if((imold.ne.iim).or.(jmold.ne.jjm)) then
+       write(*,*)'lect_start_archive: WARNING: horizontal interpolation 
+     &of thermal inertia; might be better to reset it.'
+       write(*,*)
+       endif
+       
+        ! Do horizontal interpolation
+	if (depthinterpol) then
+	  call interp_horiz(inertiedatoldnew,inertiedatS,
+     &                  imold,jmold,iim,jjm,nsoilmx,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+	else
+          call interp_horiz(inertiedatold,inertiedatS,
+     &                  imold,jmold,iim,jjm,nsoilold,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+        endif ! of if (depthinterpol)
+
+      else ! no 3D thermal inertia data
+       write(*,*)'lect_start_archive: using reference surface inertia'
+        ! Use surface inertia (and extend it to all depths)
+        do i=1,nsoilmx
+         inertiedatS(1:iip1,1:jjp1,i)=surfith(1:iip1,1:jjp1)
+        enddo
+	! Build an old resolution surface thermal inertia
+	! (will be needed for tsoil interpolation)
+	call interp_horiz(surfith,surfithold,
+     &                    iim,jjm,imold,jmold,1,
+     &                    rlonu,rlatv,rlonuold,rlatvold)
+      endif
+
+
+      ! Reshape inertiedatS to scalar grid as inertiedat
+      call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,
+     &                  inertiedatS,inertiedat)
+      
+c-----------------------------------------------------------------------
+c 6.2.2 Soil temperature
+c-----------------------------------------------------------------------
+!      write(*,*) 'Soil'
+
+      !print*,'Problem in lect_start_archive interpolating'
+      !print*,'to new resolution!!'
+
+      ! Recast temperatures along soil depth, if necessary
+      if (olddepthdef) then
+        allocate(oldgrid(nsoilold+1))
+        allocate(oldval(nsoilold+1))
+	allocate(newval(nsoilmx))
+        do i=1,imold+1
+	 do j=1,jmold+1
+
+            !if(i.gt.iip1 .or. j.gt.jjp1)then
+               !print*,'Problem in lect_start_archive interpolating'
+               !print*,'to new resolution!!'
+               !call abort
+            !endif
+
+	   ! copy values
+	   oldval(1)=tsurfold(i,j)
+!	   oldval(1)=tsurfS(i,j)
+	   oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
+	   ! build vertical coordinate
+	   oldgrid(1)=0. ! ground
+	   oldgrid(2:nsoilold+1)=mlayerold(1:nsoilold)*
+     &                (surfithold(i,j)/1.e6)
+          ! Note; at this stage, we impose volcapa=1.e6 above
+	  ! since volcapa isn't set in old soil definitions
+
+	  ! interpolate
+	  call interp_line(oldgrid,oldval,nsoilold+1,
+     &                     mlayer,newval,nsoilmx)
+	 ! copy result in tsoilold
+	 tsoiloldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
+	 enddo
+	enddo
+        ! cleanup
+	deallocate(oldgrid)
+	deallocate(oldval)
+	deallocate(newval)
+
+      else
+       if (depthinterpol) then ! if vertical interpolation is required
+        allocate(oldgrid(nsoilold+1))
+        allocate(oldval(nsoilold+1))
+	allocate(newval(nsoilmx))
+        ! build vertical coordinate
+	oldgrid(1)=0. ! ground
+	oldgrid(2:nsoilold+1)=mlayerold(1:nsoilold)
+        do i=1,imold+1
+	 do j=1,jmold+1
+	   ! copy values
+	   oldval(1)=tsurfold(i,j)
+!	   oldval(1)=tsurfS(i,j)
+	   oldval(2:nsoilold+1)=tsoilold(i,j,1:nsoilold)
+	  ! interpolate
+	  call interp_line(oldgrid,oldval,nsoilold+1,
+     &                     mlayer,newval,nsoilmx)
+	 ! copy result in tsoilold
+	 tsoiloldnew(i,j,1:nsoilmx)=newval(1:nsoilmx)
+	 enddo
+	enddo
+!	write(*,*)'tsoiloldnew(1,1,1):',tsoiloldnew(1,1,1)
+        ! cleanup
+	deallocate(oldgrid)
+	deallocate(oldval)
+	deallocate(newval)
+       
+       else
+        tsoiloldnew(:,:,:)=tsoilold(:,:,:)
+       endif ! of if (depthinterpol)
+      endif ! of if (olddepthdef)
+
+      ! Do the horizontal interpolation
+       call interp_horiz(tsoiloldnew,tsoilS,
+     &                  imold,jmold,iim,jjm,nsoilmx,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+
+      ! Reshape tsoilS to scalar grid as tsoil
+       call gr_dyn_fi (nsoilmx,iim+1,jjm+1,ngrid,tsoilS,tsoil)
+
+c-----------------------------------------------------------------------
+c 6.3   Slab Ocean :
+c-----------------------------------------------------------------------
+      call interp_horiz (tslabold,tslabs,imold,jmold,iim,jjm,noceanmx,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (noceanmx,iim+1,jjm+1,ngrid,tslabs,tslab)
+
+      call interp_horiz (rnatold,rnats,imold,jmold,iim,jjm,1,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,rnats,rnat)
+
+      call interp_horiz (pctsrf_sicold,pctsrf_sics,imold,jmold,iim,
+     &                   jjm,1,rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,pctsrf_sics,pctsrf_sic)
+
+      call interp_horiz (tsea_iceold,tsea_ices,imold,jmold,iim,jjm,1,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tsea_ices,tsea_ice)
+
+      call interp_horiz (sea_iceold,sea_ices,imold,jmold,iim,jjm,1,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,sea_ices,sea_ice)
+
+c-----------------------------------------------------------------------
+c 6.4 Variable 3d :
+c-----------------------------------------------------------------------
+      
+c temperatures atmospheriques
+      write (*,*) 'lect_start_archive: told ', told (1,jmold+1,1)  ! INFO
+      call interp_vert
+     &    (told,var,lmold,llm,apsold,bpsold,aps,bps,
+     &     psold,(imold+1)*(jmold+1))
+      write (*,*) 'lect_start_archive: var ', var (1,jmold+1,1)  ! INFO
+      call interp_horiz(var,t,imold,jmold,iim,jjm,llm,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+      write (*,*) 'lect_start_archive: t ', t(1,jjp1,1)  ! INFO
+
+c q2 : pbl wind variance
+      write (*,*) 'lect_start_archive: q2old ', q2old (1,2,1)  ! INFO
+      call interp_vert (q2old,varp1,lmold+1,llm+1,
+     &     apsold,bpsold,ap,bp,psold,(imold+1)*(jmold+1))
+      write (*,*) 'lect_start_archive: varp1 ', varp1 (1,2,1)  ! INFO
+      call interp_horiz(varp1,q2s,imold,jmold,iim,jjm,llm+1,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+      write (*,*) 'lect_start_archive: q2s ', q2s (1,2,1)  ! INFO
+      call gr_dyn_fi (llm+1,iim+1,jjm+1,ngrid,q2s,q2)
+      write (*,*) 'lect_start_archive: q2 ', q2 (1,2)  ! INFO
+c     write(47,*) 'q2',q2
+
+c calcul des champ de vent; passage en vent covariant
+      write (*,*) 'lect_start_archive: uold ', uold (1,2,1)  ! INFO
+      call interp_vert
+     & (uold,var,lmold,llm,apsold,bpsold,aps,bps,
+     &  psold,(imold+1)*(jmold+1))
+      write (*,*) 'lect_start_archive: var ', var (1,2,1)  ! INFO
+      call interp_horiz(var,us,imold,jmold,iim,jjm,llm,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+      write (*,*) 'lect_start_archive: us ', us (1,2,1)   ! INFO
+
+      call interp_vert
+     & (vold,var,lmold,llm,
+     &  apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
+      call interp_horiz(var,vs,imold,jmold,iim,jjm,llm,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+      call scal_wind(us,vs,unat,vnat)
+      write (*,*) 'lect_start_archive: unat ', unat (1,2,1)    ! INFO
+      do l=1,llm
+        do j = 1, jjp1
+          do i=1,iip1
+            ucov( i,j,l ) = unat( i,j,l ) * cu(i,j)
+c           ucov( i,j,l ) = 0
+          end do
+        end do
+      end do 
+      write (*,*) 'lect_start_archive: ucov ', ucov (1,2,1)  ! INFO
+c     write(48,*) 'ucov',ucov
+      do l=1,llm
+        do j = 1, jjm
+          do i=1,iim
+            vcov( i,j,l ) = vnat( i,j,l ) * cv(i,j)
+c           vcov( i,j,l ) = 0
+          end do
+          vcov( iip1,j,l ) = vcov( 1,j,l )
+        end do
+      end do
+c     write(49,*) 'ucov',vcov
+
+      if (nqtot .gt. 0) then
+c traceurs surface
+      do iq = 1, nqtot
+            call interp_horiz(qsurfold(1,1,iq) ,qsurfs(1,1,iq),
+     &                  imold,jmold,iim,jjm,1,
+     &                  rlonuold,rlatvold,rlonu,rlatv)
+      enddo
+
+      call gr_dyn_fi (nqtot,iim+1,jjm+1,ngrid,qsurfs,qsurf)
+
+c traceurs 3D
+      do  iq = 1, nqtot
+            call interp_vert(qold(1,1,1,iq),var,lmold,llm,
+     &        apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
+            call interp_horiz(var,q(1,1,1,iq),imold,jmold,iim,jjm,llm,
+     &                  rlonuold,rlatvold,rlonu,rlatv)
+      enddo
+cccccccccccccccccccccccccccccc      
+c  make sure that sum of q = 1      
+c dominent species is = 1 - sum(all other species)      
+cccccccccccccccccccccccccccccc      
+c      iqmax=1
+c      
+c      if (nqold.gt.10) then
+c       do l=1,llm
+c        do j=1,jjp1
+c          do i=1,iip1
+c           do iq=1,nqold
+c            if (q(i,j,l,iq).gt.q(i,j,l,iqmax)) then
+c              iqmax=iq
+c            endif
+c           enddo
+c           q(i,j,l,iqmax)=1.
+c           qtot(i,j,l)=0
+c           do iq=1,nqold
+c            if (iq.ne.iqmax) then        
+c              q(i,j,l,iqmax)=q(i,j,l,iqmax)-q(i,j,l,iq)        
+c            endif
+c           enddo !iq
+c           do iq=1,nqold
+c            qtot(i,j,l)=qtot(i,j,l)+q(i,j,l,iq)
+c            if (i.eq.1.and.j.eq.1.and.l.Eq.1) write(*,*)' qtot(i,j,l)',
+c     $    qtot(i,j,l)
+c           enddo !iq
+c          enddo !i   
+c         enddo !j   
+c       enddo !l  
+c      endif
+ccccccccccccccccccccccccccccccc
+
+c     Periodicite :
+      do  iq = 1, nqtot
+         do l=1, llm
+            do j = 1, jjp1
+               q(iip1,j,l,iq) = q(1,j,l,iq)
+            end do
+         end do
+      enddo
+      
+!      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,co2ices,co2ice)
+! no need to transfer "co2ice" any more; it is in qsurf(igcm_co2_ice)
+
+      endif !! if nqtot .ne. 0
+
+c-----------------------------------------------------------------------
+c   Initialisation  h:	(passage de t -> h)
+c-----------------------------------------------------------------------
+
+      DO l=1,llm
+         DO j=1,jjp1
+            DO i=1,iim
+               h(i,j,l) = t(i,j,l)*((ps(i,j)/preff)**kappa)
+            END DO
+            h(iip1,j,l) =  h(1,j,l)
+         END DO
+      END DO
+
+
+c***********************************************************************
+c***********************************************************************
+c     Fin subroutine lecture ini
+c***********************************************************************
+c***********************************************************************
+
+      deallocate(timelist)
+      deallocate(rlonuold, rlatvold)
+      deallocate(rlonvold, rlatuold)
+      deallocate(apsold,bpsold)
+      deallocate(uold)
+      deallocate(vold)
+      deallocate(told)
+      deallocate(psold)
+      deallocate(phisold)
+      deallocate(qold)
+      deallocate(co2iceold)
+      deallocate(tsurfold)
+      deallocate(emisold)
+      deallocate(q2old)
+      deallocate(tsoilold)
+      deallocate(tsoiloldnew)
+      deallocate(inertiedatold)
+      deallocate(inertiedatoldnew)
+      deallocate(surfithold)
+      deallocate(mlayerold)
+      deallocate(qsurfold)
+      deallocate(var,varp1)
+      deallocate(tslabold)
+      deallocate(rnatold)
+      deallocate(pctsrf_sicold)
+      deallocate(tsea_iceold)
+      deallocate(sea_iceold)
+
+!      write(*,*)'lect_start_archive: END'
+      return
+      end
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/newstart.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/newstart.F	(revision 1644)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/newstart.F	(revision 1644)
@@ -0,0 +1,1667 @@
+C======================================================================
+      PROGRAM newstart
+c=======================================================================
+c
+c
+c   Auteur:   Christophe Hourdin/Francois Forget/Yann Wanherdrick
+c   ------
+c             Derniere modif : 12/03
+c
+c
+c   Objet:  Create or modify the initial state for the LMD Mars GCM
+c   -----           (fichiers NetCDF start et startfi)
+c
+c
+c=======================================================================
+
+      use mod_phys_lmdz_para, only: is_parallel, is_sequential,
+     &                              is_mpi_root, is_omp_root,
+     &                              is_master
+      use infotrac, only: infotrac_init, nqtot, tname
+      USE tracer_h, ONLY: igcm_co2_ice, igcm_h2o_vap, igcm_h2o_ice
+      USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat
+      USE surfdat_h, ONLY: phisfi, albedodat,
+     &                     zmea, zstd, zsig, zgam, zthe
+      use datafile_mod, only: datadir, surfdir
+      use ioipsl_getin_p_mod, only: getin_p
+      use control_mod, only: day_step, iphysiq, anneeref, planet_type
+      use phyredem, only: physdem0, physdem1
+      use iostart, only: open_startphy
+      use slab_ice_h, only:noceanmx
+      use filtreg_mod, only: inifilr
+      USE mod_const_mpi, ONLY: COMM_LMDZ
+      USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff
+      USE comconst_mod, ONLY: lllm,daysec,dtvr,dtphys,cpp,kappa,
+     .                        rad,omeg,g,r,pi
+      USE serre_mod, ONLY: alphax
+      USE temps_mod, ONLY: day_ini
+      USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+      use iniphysiq_mod, only: iniphysiq
+      implicit none
+
+      include "dimensions.h"
+      integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm) 
+      include "paramet.h"
+      include "comgeom2.h"
+      include "comdissnew.h"
+      include "netcdf.inc"
+
+c=======================================================================
+c   Declarations
+c=======================================================================
+
+c Variables dimension du fichier "start_archive"
+c------------------------------------
+      CHARACTER        relief*3
+
+
+c Variables pour les lectures NetCDF des fichiers "start_archive" 
+c--------------------------------------------------
+      INTEGER nid_dyn, nid_fi,nid,nvarid
+      INTEGER length
+      parameter (length = 100)
+      INTEGER tab0
+      INTEGER   NB_ETATMAX
+      parameter (NB_ETATMAX = 100)
+
+      REAL  date
+      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
+
+c Variable histoire 
+c------------------
+      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
+      REAL phis(iip1,jjp1)
+      REAL,ALLOCATABLE :: q(:,:,:,:)               ! champs advectes
+
+c autre variables dynamique nouvelle grille
+c------------------------------------------
+      REAL pks(iip1,jjp1)
+      REAL w(iip1,jjp1,llm+1)
+      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
+!      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
+!      REAL dh(ip1jmp1,llm),dp(ip1jmp1)
+      REAL phi(iip1,jjp1,llm)
+
+      integer klatdat,klongdat
+      PARAMETER (klatdat=180,klongdat=360)
+
+c Physique sur grille scalaire 
+c----------------------------
+      real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
+      real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
+
+c variable physique
+c------------------
+      REAL tsurf(ngridmx)        ! surface temperature
+      REAL,ALLOCATABLE :: tsoil(:,:) ! soil temperature
+!      REAL co2ice(ngridmx)        ! CO2 ice layer
+      REAL emis(ngridmx)        ! surface emissivity
+      real emisread             ! added by RW
+      REAL,ALLOCATABLE :: qsurf(:,:)
+      REAL q2(ngridmx,llm+1)
+!      REAL rnaturfi(ngridmx)
+      real alb(iip1,jjp1),albfi(ngridmx) ! albedos
+      real,ALLOCATABLE :: ith(:,:,:),ithfi(:,:) ! thermal inertia (3D)
+      real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D)
+      REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
+
+      INTEGER i,j,l,isoil,ig,idum
+      real mugaz ! molar mass of the atmosphere
+
+      integer ierr 
+
+      REAL rnat(ngridmx)
+      REAL,ALLOCATABLE :: tslab(:,:) ! slab ocean temperature
+      REAL pctsrf_sic(ngridmx) ! sea ice cover
+      REAL tsea_ice(ngridmx) ! temperature sea_ice
+      REAL sea_ice(ngridmx) ! mass of sea ice (kg/m2)
+
+c Variables on the new grid along scalar points 
+c------------------------------------------------------
+!      REAL p(iip1,jjp1)
+      REAL t(iip1,jjp1,llm)
+      REAL tset(iip1,jjp1,llm)
+      real phisold_newgrid(iip1,jjp1)
+      REAL :: teta(iip1, jjp1, llm)
+      REAL :: pk(iip1,jjp1,llm)
+      REAL :: pkf(iip1,jjp1,llm)
+      REAL :: ps(iip1, jjp1)
+      REAL :: masse(iip1,jjp1,llm)
+      REAL :: xpn,xps,xppn(iim),xpps(iim)
+      REAL :: p3d(iip1, jjp1, llm+1)
+      REAL :: beta(iip1,jjp1,llm)
+!      REAL dteta(ip1jmp1,llm)
+
+c Variable de l'ancienne grille 
+c------------------------------
+      real time
+      real tab_cntrl(100)
+      real tab_cntrl_bis(100)
+
+c variables diverses
+c-------------------
+      real choix_1,pp
+      character*80      fichnom
+      character*250  filestring
+      integer Lmodif,iq
+      character modif*20
+      real z_reel(iip1,jjp1)
+      real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove
+      real ptoto,pcap,patm,airetot,ptotn,patmn,psea
+!      real ssum
+      character*1 yes
+      logical :: flagtset=.false. ,  flagps0=.false.
+      real val, val2, val3, val4 ! to store temporary variables
+      real :: iceith=2000 ! thermal inertia of subterranean ice
+
+      INTEGER :: itau
+      
+      character(len=20) :: txt ! to store some text
+      character(len=50) :: surfacefile ! "surface.nc" (or equivalent file)
+      character(len=150) :: longline
+      integer :: count
+      real :: profile(llm+1) ! to store an atmospheric profile + surface value
+
+!     added by BC for equilibrium temperature startup
+      real teque
+
+!     added by BC for cloud fraction setup
+      REAL hice(ngridmx),cloudfrac(ngridmx,llm)
+      REAL totalfrac(ngridmx)
+
+
+!     added by RW for nuketharsis
+      real fact1
+      real fact2
+
+
+c sortie visu pour les champs dynamiques
+c---------------------------------------
+!      INTEGER :: visuid
+!      real :: time_step,t_ops,t_wrt
+!      CHARACTER*80 :: visu_file
+
+      cpp    = 0.
+      preff  = 0.
+      pa     = 0. ! to ensure disaster rather than minor error if we don`t
+                  ! make deliberate choice of these values elsewhere.
+
+      planet_type="titan"
+
+! initialize "serial/parallel" related stuff
+! (required because we call tabfi() below, before calling iniphysiq)
+      is_sequential=.true.
+      is_parallel=.false.
+      is_mpi_root=.true.
+      is_omp_root=.true.
+      is_master=.true.
+      
+! Load tracer number and names:
+      call infotrac_init
+! allocate arrays
+      allocate(q(iip1,jjp1,llm,nqtot))
+      allocate(qsurf(ngridmx,nqtot))
+      
+! get value of nsoilmx and allocate corresponding arrays
+      ! a default value of nsoilmx is set in comsoil_h
+      call getin_p("nsoilmx",nsoilmx)
+      
+      allocate(tsoil(ngridmx,nsoilmx))
+      allocate(ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx))
+      allocate(tslab(ngridmx,nsoilmx))
+      
+c=======================================================================
+c   Choice of the start file(s) to use
+c=======================================================================
+      write(*,*) 'From which kind of files do you want to create new',
+     .  'start and startfi files'
+      write(*,*) '    0 - from a file start_archive'
+      write(*,*) '    1 - from files start and startfi'
+ 
+c-----------------------------------------------------------------------
+c   Open file(s) to modify (start or start_archive)
+c-----------------------------------------------------------------------
+
+      DO
+         read(*,*,iostat=ierr) choix_1
+         if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT
+      ENDDO
+
+c     Open start_archive
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~
+      if (choix_1.eq.0) then
+
+        write(*,*) 'Creating start files from:'
+        write(*,*) './start_archive.nc'
+        write(*,*)
+        fichnom = 'start_archive.nc'
+        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
+        IF (ierr.NE.NF_NOERR) THEN
+          write(6,*)' Problem opening file:',fichnom
+          write(6,*)' ierr = ', ierr
+          CALL ABORT
+        ENDIF
+        tab0 = 50 
+        Lmodif = 1
+
+c     OR open start and startfi files
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      else
+        write(*,*) 'Creating start files from:'
+        write(*,*) './start.nc and ./startfi.nc'
+        write(*,*)
+        fichnom = 'start.nc'
+        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn)
+        IF (ierr.NE.NF_NOERR) THEN
+          write(6,*)' Problem opening file:',fichnom
+          write(6,*)' ierr = ', ierr
+          CALL ABORT
+        ENDIF
+ 
+        fichnom = 'startfi.nc'
+        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi)
+        IF (ierr.NE.NF_NOERR) THEN
+          write(6,*)' Problem opening file:',fichnom
+          write(6,*)' ierr = ', ierr
+          CALL ABORT
+        ENDIF
+
+        tab0 = 0 
+        Lmodif = 0
+
+      endif
+
+
+c=======================================================================
+c  INITIALISATIONS DIVERSES
+c=======================================================================
+
+! Initialize global tracer indexes (stored in tracer.h)
+! ... this has to be done before phyetat0
+      call initracer(ngridmx,nqtot,tname)
+
+! Take care of arrays in common modules
+      ! ALLOCATE ARRAYS in surfdat_h (if not already done, e.g. when using start_archive)
+      IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngridmx))
+      IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngridmx))
+      IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngridmx))
+      IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngridmx))
+      IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngridmx))
+      IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngridmx))
+      IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngridmx))
+
+c-----------------------------------------------------------------------
+c Lecture du tableau des parametres du run (pour la dynamique)
+c-----------------------------------------------------------------------
+      if (choix_1.eq.0) then
+
+        write(*,*) 'reading tab_cntrl START_ARCHIVE'
+c
+        ierr = NF_INQ_VARID (nid, "controle", nvarid)
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
+#else
+        ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
+#endif
+c
+      else if (choix_1.eq.1) then
+
+        write(*,*) 'reading tab_cntrl START'
+c
+        ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid)
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl)
+#else
+        ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl)
+#endif
+c
+        write(*,*) 'reading tab_cntrl STARTFI'
+c
+        ierr = NF_INQ_VARID (nid_fi, "controle", nvarid)
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis)
+#else
+        ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis)
+#endif
+c
+        do i=1,50
+          tab_cntrl(i+50)=tab_cntrl_bis(i)
+        enddo
+        write(*,*) 'printing tab_cntrl', tab_cntrl
+        do i=1,100
+          write(*,*) i,tab_cntrl(i)
+        enddo
+        
+        write(*,*) 'Reading file START'
+        fichnom = 'start.nc'
+        CALL dynetat0(fichnom,vcov,ucov,teta,q,masse,
+     .       ps,phis,time)
+
+        CALL iniconst
+        CALL inigeom
+
+! Initialize the physics
+         CALL iniphysiq(iim,jjm,llm,
+     &                  (jjm-1)*iim+2,comm_lmdz,
+     &                  daysec,day_ini,dtphys,
+     &                  rlatu,rlatv,rlonu,rlonv,
+     &                  aire,cu,cv,rad,g,r,cpp,
+     &                  1)
+
+        ! Lmodif set to 0 to disable modifications possibility in phyeta0                           
+        Lmodif=0
+        write(*,*) 'Reading file STARTFI'
+        fichnom = 'startfi.nc'
+        CALL phyetat0 (ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx,
+     .        nqtot,day_ini,time,
+     .        tsurf,tsoil,emis,q2,qsurf,   !) ! temporary modif by RDW
+     .        cloudfrac,totalfrac,hice,rnat,pctsrf_sic,tslab,tsea_ice,
+     .        sea_ice)
+
+        ! copy albedo and soil thermal inertia on (local) physics grid
+        do i=1,ngridmx
+          albfi(i) = albedodat(i)
+          do j=1,nsoilmx
+           ithfi(i,j) = inertiedat(i,j)
+          enddo
+        ! build a surfithfi(:) using 1st layer of ithfi(:), which might
+        ! be needed later on if reinitializing soil thermal inertia
+          surfithfi(i)=ithfi(i,1)
+        enddo
+        ! also copy albedo and soil thermal inertia on (local) dynamics grid
+        ! so that options below can manipulate either (but must then ensure
+        ! to correctly recast things on physics grid)
+        call gr_fi_dyn(1,ngridmx,iip1,jjp1,albfi,alb)
+        call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
+        call gr_fi_dyn(1,ngridmx,iip1,jjp1,surfithfi,surfith)
+      
+      endif
+c-----------------------------------------------------------------------
+c                Initialisation des constantes dynamique
+c-----------------------------------------------------------------------
+
+      kappa = tab_cntrl(9) 
+      etot0 = tab_cntrl(12)
+      ptot0 = tab_cntrl(13)
+      ztot0 = tab_cntrl(14)
+      stot0 = tab_cntrl(15)
+      ang0 = tab_cntrl(16)
+      write(*,*) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0"
+      write(*,*) kappa,etot0,ptot0,ztot0,stot0,ang0
+
+      ! for vertical coordinate
+      preff=tab_cntrl(18) ! reference surface pressure
+      pa=tab_cntrl(17)  ! reference pressure at which coord is purely pressure
+      !NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0
+      write(*,*) "Newstart: preff=",preff," pa=",pa
+      yes=' '
+      do while ((yes.ne.'y').and.(yes.ne.'n'))
+        write(*,*) "Change the values of preff and pa? (y/n)"
+        read(*,fmt='(a)') yes
+      end do
+      if (yes.eq.'y') then
+        write(*,*)"New value of reference surface pressure preff?"
+        write(*,*)"   (for Mars, typically preff=610)"
+        read(*,*) preff
+        write(*,*)"New value of reference pressure pa for purely"
+        write(*,*)"pressure levels (for hybrid coordinates)?"
+        write(*,*)"   (for Mars, typically pa=20)"
+        read(*,*) pa
+      endif
+c-----------------------------------------------------------------------
+c   Lecture du tab_cntrl et initialisation des constantes physiques
+c  - pour start:  Lmodif = 0 => pas de modifications possibles
+c                  (modif dans le tabfi de readfi + loin)
+c  - pour start_archive:  Lmodif = 1 => modifications possibles
+c-----------------------------------------------------------------------
+      if (choix_1.eq.0) then
+         ! tabfi requires that input file be first opened by open_startphy(fichnom)
+         fichnom = 'start_archive.nc'
+         call open_startphy(fichnom)
+         call tabfi (ngridmx,nid,Lmodif,tab0,day_ini,lllm,p_rad,
+     .            p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
+      else if (choix_1.eq.1) then
+         fichnom = 'startfi.nc'
+         call open_startphy(fichnom)
+         Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0                           
+         call tabfi (ngridmx,nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
+     .            p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
+      endif
+
+      if (p_omeg .eq. -9999.) then
+        p_omeg = 8.*atan(1.)/p_daysec
+        print*,"new value of omega",p_omeg
+      endif
+
+      rad = p_rad
+      omeg = p_omeg
+      g = p_g
+      cpp = p_cpp
+      mugaz = p_mugaz
+      daysec = p_daysec
+
+      if (p_omeg .eq. -9999.) then
+        p_omeg = 8.*atan(1.)/p_daysec
+        print*,"new value of omega",p_omeg
+      endif
+
+      kappa = 8.314*1000./(p_mugaz*p_cpp) ! added by RDW
+
+c=======================================================================
+c  INITIALISATIONS DIVERSES
+c=======================================================================
+! Load parameters from run.def file
+      CALL defrun_new( 99, .TRUE. )
+! Initialize dynamics geometry and co. (which, when using start.nc may 
+! have changed because of modifications (tabi, preff, pa) above)
+      CALL iniconst 
+      CALL inigeom
+      idum=-1
+      idum=0
+
+! Initialize the physics for start_archive only
+      IF (choix_1.eq.0) THEN
+         CALL iniphysiq(iim,jjm,llm,
+     &                  (jjm-1)*iim+2,comm_lmdz,
+     &                  daysec,day_ini,dtphys,
+     &                  rlatu,rlatv,rlonu,rlonv,
+     &                  aire,cu,cv,rad,g,r,cpp,
+     &                  1)
+      ENDIF
+
+c=======================================================================
+c   lecture topographie, albedo, inertie thermique, relief sous-maille
+c=======================================================================
+
+      if (choix_1.eq.0) then  ! for start_archive files, 
+                              ! where an external "surface.nc" file is needed
+
+c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla'))
+c       write(*,*)
+c       write(*,*) 'choix du relief (mola,pla)'
+c       write(*,*) '(Topographie MGS MOLA, plat)'
+c       read(*,fmt='(a3)') relief
+        relief="mola"
+c     enddo
+
+!    First get the correct value of datadir, if not already done:
+        ! default 'datadir' is set in "datafile_mod"
+        call getin_p("datadir",datadir)
+        write(*,*) 'Available surface data files are:'
+        filestring='ls '//trim(datadir)//'/'//
+     &                    trim(surfdir)//' | grep .nc'
+        call system(filestring)
+        ! but in ye old days these files were in datadir, so scan it as well
+        ! for the sake of retro-compatibility
+        filestring='ls '//trim(datadir)//'/'//' | grep .nc'
+        call system(filestring)
+
+        write(*,*) ''
+        write(*,*) 'Please choose the relevant file',
+     &  ' (e.g. "surface_mars.nc")'
+        write(*,*) ' or "none" to not use any (and set planetary'
+        write(*,*) '  albedo and surface thermal inertia)'
+        read(*,fmt='(a50)') surfacefile
+
+        if (surfacefile.ne."none") then
+
+          CALL datareadnc(relief,surfacefile,phis,alb,surfith,
+     &          zmeaS,zstdS,zsigS,zgamS,ztheS)
+        else
+        ! specific case when not using a "surface.nc" file
+          phis(:,:)=0
+          zmeaS(:,:)=0
+          zstdS(:,:)=0
+          zsigS(:,:)=0
+          zgamS(:,:)=0
+          ztheS(:,:)=0
+          
+          write(*,*) "Enter value of albedo of the bare ground:"
+          read(*,*) alb(1,1)
+          alb(:,:)=alb(1,1)
+          
+          write(*,*) "Enter value of thermal inertia of soil:"
+          read(*,*) surfith(1,1)
+          surfith(:,:)=surfith(1,1)
+        
+        endif ! of if (surfacefile.ne."none")
+
+        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
+        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
+        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
+        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
+        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig)
+        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
+        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
+
+      endif ! of if (choix_1.eq.0)
+
+
+c=======================================================================
+c  Lecture des fichiers (start ou start_archive)
+c=======================================================================
+
+      if (choix_1.eq.0) then
+
+        write(*,*) 'Reading file START_ARCHIVE'
+        CALL lect_start_archive(ngridmx,llm,
+     &   date,tsurf,tsoil,emis,q2,
+     &   t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
+     &   surfith,nid,
+     &   rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
+        write(*,*) "OK, read start_archive file"
+        ! copy soil thermal inertia
+        ithfi(:,:)=inertiedat(:,:)
+        
+        ierr= NF_CLOSE(nid)
+
+      else if (choix_1.eq.1) then 
+         !do nothing, start and startfi have already been read
+      else 
+        CALL exit(1)
+      endif
+
+      dtvr   = daysec/FLOAT(day_step)
+      dtphys   = dtvr * FLOAT(iphysiq)
+
+c=======================================================================
+c 
+c=======================================================================
+
+      do ! infinite loop on list of changes
+
+      write(*,*)
+      write(*,*)
+      write(*,*) 'List of possible changes :'
+      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+      write(*,*)
+      write(*,*) 'flat : no topography ("aquaplanet")'
+      write(*,*) 'set_ps_to_preff : used if changing preff with topo'
+      write(*,*) 'nuketharsis : no Tharsis bulge'
+      write(*,*) 'bilball : uniform albedo and thermal inertia'
+      write(*,*) 'coldspole : cold subsurface and high albedo at S.pole'
+      write(*,*) 'qname : change tracer name'
+      write(*,*) 't=profile  : read temperature profile in profile.in'
+      write(*,*) 'q=0 : ALL tracer =zero'
+      write(*,*) 'q=x : give a specific uniform value to one tracer'
+      write(*,*) 'q=profile    : specify a profile for a tracer'
+!      write(*,*) 'ini_q : tracers initialisation for chemistry, water an
+!     $d ice   '
+!      write(*,*) 'ini_q-H2O : tracers initialisation for chemistry and 
+!     $ice '
+!      write(*,*) 'ini_q-iceH2O : tracers initialisation for chemistry on
+!     $ly '
+      write(*,*) 'noglacier : Remove tropical H2O ice if |lat|<45'
+      write(*,*) 'watercapn : H20 ice on permanent N polar cap '
+      write(*,*) 'watercaps : H20 ice on permanent S polar cap '
+      write(*,*) 'noacglac  : H2O ice across Noachis Terra'
+      write(*,*) 'oborealis : H2O ice across Vastitas Borealis'
+      write(*,*) 'iceball   : Thick ice layer all over surface'
+      write(*,*) 'supercontinent: Create a continent of given Ab and TI'
+      write(*,*) 'wetstart  : start with a wet atmosphere'
+      write(*,*) 'isotherm  : Isothermal Temperatures, wind set to zero'
+      write(*,*) 'radequi   : Earth-like radiative equilibrium temperature
+     $ profile (lat-alt) and winds set to zero'
+      write(*,*) 'coldstart : Start X K above the CO2 frost point and 
+     $set wind to zero (assumes 100% CO2)'
+      write(*,*) 'co2ice=0 : remove CO2 polar cap'
+      write(*,*) 'ptot : change total pressure'
+      write(*,*) 'emis : change surface emissivity'
+      write(*,*) 'therm_ini_s : Set soil thermal inertia to reference sur
+     &face values'
+      write(*,*) 'slab_ocean_0 : initialisation of slab 
+     $ocean variables'
+
+        write(*,*)
+        write(*,*) 'Change to perform ?'
+        write(*,*) '   (enter keyword or return to end)'
+        write(*,*)
+
+        read(*,fmt='(a20)') modif
+        if (modif(1:1) .eq. ' ')then
+         exit ! exit loop on changes
+        endif
+
+        write(*,*)
+        write(*,*) trim(modif) , ' : '
+
+c       'flat : no topography ("aquaplanet")'
+c       -------------------------------------
+        if (trim(modif) .eq. 'flat') then
+c         set topo to zero 
+          z_reel(1:iip1,1:jjp1)=0
+          phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1)
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+          write(*,*) 'topography set to zero.'
+          write(*,*) 'WARNING : the subgrid topography parameters',
+     &    ' were not set to zero ! => set calllott to F'                    
+
+c        Choice of surface pressure
+         yes=' '
+         do while ((yes.ne.'y').and.(yes.ne.'n'))
+            write(*,*) 'Do you wish to choose homogeneous surface',
+     &                 'pressure (y) or let newstart interpolate ',
+     &                 ' the previous field  (n)?'
+             read(*,fmt='(a)') yes
+         end do
+         if (yes.eq.'y') then
+           flagps0=.true.
+           write(*,*) 'New value for ps (Pa) ?'
+ 201       read(*,*,iostat=ierr) patm
+            if(ierr.ne.0) goto 201
+            write(*,*) patm
+            if (patm.eq.-9999.) then
+              patm = preff
+              write(*,*) "we set patm = preff", preff
+            endif
+             write(*,*)
+             write(*,*) ' new ps everywhere (Pa) = ', patm
+             write(*,*)
+             do j=1,jjp1
+               do i=1,iip1
+                 ps(i,j)=patm
+               enddo
+             enddo
+         end if
+
+c       'set_ps_to_preff' : used if changing preff with topo  
+c       ----------------------------------------------------
+        else if (trim(modif) .eq. 'set_ps_to_preff') then
+          do j=1,jjp1
+           do i=1,iip1
+             ps(i,j)=preff
+           enddo
+          enddo
+
+c       'nuketharsis : no tharsis bulge for Early Mars'
+c       ---------------------------------------------
+        else if (trim(modif) .eq. 'nuketharsis') then
+
+           DO j=1,jjp1        
+              DO i=1,iim
+                 ig=1+(j-2)*iim +i
+                 if(j.eq.1) ig=1
+                 if(j.eq.jjp1) ig=ngridmx
+
+                 fact1=(((rlonv(i)*180./pi)+100)**2 + 
+     &                (rlatu(j)*180./pi)**2)/65**2 
+                 fact2=exp( -fact1**2.5 )
+
+                 phis(i,j) = phis(i,j) - (phis(i,j)+4000.*g)*fact2
+
+!                 if(phis(i,j).gt.2500.*g)then
+!                    if(rlatu(j)*180./pi.gt.-80.)then ! avoid chopping south polar cap
+!                       phis(i,j)=2500.*g
+!                    endif
+!                 endif
+
+              ENDDO
+           ENDDO
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+
+
+c       bilball : uniform albedo, thermal inertia
+c       -----------------------------------------
+        else if (trim(modif) .eq. 'bilball') then
+          write(*,*) 'constante albedo and iner.therm:'
+          write(*,*) 'New value for albedo (ex: 0.25) ?'
+ 101      read(*,*,iostat=ierr) alb_bb
+          if(ierr.ne.0) goto 101
+          write(*,*)
+          write(*,*) ' uniform albedo (new value):',alb_bb
+          write(*,*)
+
+          write(*,*) 'New value for thermal inertia (eg: 247) ?'
+ 102      read(*,*,iostat=ierr) ith_bb
+          if(ierr.ne.0) goto 102
+          write(*,*) 'uniform thermal inertia (new value):',ith_bb
+          DO j=1,jjp1
+             DO i=1,iip1
+                alb(i,j) = alb_bb        ! albedo
+                do isoil=1,nsoilmx
+                  ith(i,j,isoil) = ith_bb        ! thermal inertia
+                enddo
+             END DO
+          END DO
+!          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi)
+          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
+
+c       coldspole : sous-sol de la calotte sud toujours froid
+c       -----------------------------------------------------
+        else if (trim(modif) .eq. 'coldspole') then
+          write(*,*)'new value for the subsurface temperature',
+     &   ' beneath the permanent southern polar cap ? (eg: 141 K)'
+ 103      read(*,*,iostat=ierr) tsud
+          if(ierr.ne.0) goto 103
+          write(*,*)
+          write(*,*) ' new value of the subsurface temperature:',tsud
+c         nouvelle temperature sous la calotte permanente
+          do l=2,nsoilmx
+               tsoil(ngridmx,l) =  tsud
+          end do
+
+
+          write(*,*)'new value for the albedo',
+     &   'of the permanent southern polar cap ? (eg: 0.75)'
+ 104      read(*,*,iostat=ierr) albsud
+          if(ierr.ne.0) goto 104
+          write(*,*)
+
+c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c         Option 1:  only the albedo of the pole is modified :    
+          albfi(ngridmx)=albsud
+          write(*,*) 'ig=',ngridmx,'   albedo perennial cap ',
+     &    albfi(ngridmx)
+
+c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
+c          Option 2 A haute resolution : coordonnee de la vrai calotte ~    
+c           DO j=1,jjp1
+c             DO i=1,iip1
+c                ig=1+(j-2)*iim +i
+c                if(j.eq.1) ig=1
+c                if(j.eq.jjp1) ig=ngridmx
+c                if ((rlatu(j)*180./pi.lt.-84.).and.
+c     &            (rlatu(j)*180./pi.gt.-91.).and.
+c     &            (rlonv(i)*180./pi.gt.-91.).and.
+c     &            (rlonv(i)*180./pi.lt.0.))         then
+cc    albedo de la calotte permanente fixe a albsud
+c                   alb(i,j)=albsud
+c                   write(*,*) 'lat=',rlatu(j)*180./pi,
+c     &                      ' lon=',rlonv(i)*180./pi
+cc     fin de la condition sur les limites de la calotte permanente
+c                end if
+c             ENDDO
+c          ENDDO
+c      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+c         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
+
+
+c       ptot : Modification of the total pressure: ice + current atmosphere 
+c       -------------------------------------------------------------------
+        else if (trim(modif).eq.'ptot') then
+
+          ! check if we have a co2_ice surface tracer:
+          if (igcm_co2_ice.eq.0) then
+            write(*,*) " No surface CO2 ice !"
+            write(*,*) " only atmospheric pressure will be considered!"
+          endif
+c         calcul de la pression totale glace + atm actuelle
+          patm=0.
+          airetot=0.
+          pcap=0.
+          DO j=1,jjp1
+             DO i=1,iim
+                ig=1+(j-2)*iim +i
+                if(j.eq.1) ig=1
+                if(j.eq.jjp1) ig=ngridmx
+                patm = patm + ps(i,j)*aire(i,j)
+                airetot= airetot + aire(i,j)
+                if (igcm_co2_ice.ne.0) then
+                  !pcap = pcap + aire(i,j)*co2ice(ig)*g
+                  pcap = pcap + aire(i,j)*qsurf(ig,igcm_co2_ice)*g
+                endif
+             ENDDO
+          ENDDO
+          ptoto = pcap + patm
+
+          print*,'Current total pressure at surface (co2 ice + atm) ',
+     &       ptoto/airetot
+
+          print*,'new value?'
+          read(*,*) ptotn
+          ptotn=ptotn*airetot
+          patmn=ptotn-pcap
+          print*,'ptoto,patm,ptotn,patmn'
+          print*,ptoto,patm,ptotn,patmn
+          print*,'Mult. factor for pressure (atm only)', patmn/patm
+          do j=1,jjp1
+             do i=1,iip1
+                ps(i,j)=ps(i,j)*patmn/patm
+             enddo
+          enddo
+
+
+
+c        Correction pour la conservation des traceurs
+         yes=' '
+         do while ((yes.ne.'y').and.(yes.ne.'n'))
+            write(*,*) 'Do you wish to conserve tracer total mass (y)',
+     &              ' or tracer mixing ratio (n) ?'
+             read(*,fmt='(a)') yes
+         end do
+
+         if (yes.eq.'y') then
+           write(*,*) 'OK : conservation of tracer total mass'
+           DO iq =1, nqtot
+             DO l=1,llm
+               DO j=1,jjp1
+                  DO i=1,iip1
+                    q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn
+                  ENDDO
+               ENDDO
+             ENDDO
+           ENDDO
+          else
+            write(*,*) 'OK : conservation of tracer mixing ratio'
+          end if
+
+c        Correction pour la pression au niveau de la mer
+         yes=' '
+         do while ((yes.ne.'y').and.(yes.ne.'n'))
+            write(*,*) 'Do you wish fix pressure at sea level (y)',
+     &              ' or not (n) ?'
+             read(*,fmt='(a)') yes
+         end do
+
+         if (yes.eq.'y') then
+           write(*,*) 'Value?'
+                read(*,*,iostat=ierr) psea
+             DO i=1,iip1
+               DO j=1,jjp1
+                    ps(i,j)=psea
+
+                  ENDDO
+               ENDDO
+                write(*,*) 'psea=',psea
+          else
+            write(*,*) 'no'
+          end if
+
+
+c       emis : change surface emissivity (added by RW)
+c       ----------------------------------------------
+        else if (trim(modif).eq.'emis') then
+
+           print*,'new value?'
+           read(*,*) emisread
+
+           do i=1,ngridmx
+              emis(i)=emisread
+           enddo
+
+c       qname : change tracer name
+c       --------------------------
+        else if (trim(modif).eq.'qname') then
+         yes='y'
+         do while (yes.eq.'y')
+          write(*,*) 'Which tracer name do you want to change ?'
+          do iq=1,nqtot
+            write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
+          enddo
+          write(*,'(a35,i3)')
+     &            '(enter tracer number; between 1 and ',nqtot
+          write(*,*)' or any other value to quit this option)'
+          read(*,*) iq
+          if ((iq.ge.1).and.(iq.le.nqtot)) then
+            write(*,*)'Change tracer name ',trim(tname(iq)),' to ?'
+            read(*,*) txt
+            tname(iq)=txt
+            write(*,*)'Do you want to change another tracer name (y/n)?'
+            read(*,'(a)') yes 
+          else
+! inapropiate value of iq; quit this option
+            yes='n'
+          endif ! of if ((iq.ge.1).and.(iq.le.nqtot))
+         enddo ! of do while (yes.ne.'y')
+
+c       q=0 : set tracers to zero
+c       -------------------------
+        else if (trim(modif).eq.'q=0') then
+c          mise a 0 des q (traceurs)
+          write(*,*) 'Tracers set to 0 (1.E-30 in fact)'
+           DO iq =1, nqtot
+             DO l=1,llm
+               DO j=1,jjp1
+                  DO i=1,iip1
+                    q(i,j,l,iq)=1.e-30
+                  ENDDO
+               ENDDO
+             ENDDO
+           ENDDO
+
+c          set surface tracers to zero
+           DO iq =1, nqtot
+             DO ig=1,ngridmx
+                 qsurf(ig,iq)=0.
+             ENDDO
+           ENDDO
+
+c       q=x : initialise tracer manually 
+c       --------------------------------
+        else if (trim(modif).eq.'q=x') then
+             write(*,*) 'Which tracer do you want to modify ?'
+             do iq=1,nqtot
+               write(*,*)iq,' : ',trim(tname(iq))
+             enddo
+             write(*,*) '(choose between 1 and ',nqtot,')'
+             read(*,*) iq 
+             write(*,*)'mixing ratio of tracer ',trim(tname(iq)),
+     &                 ' ? (kg/kg)'
+             read(*,*) val
+             DO l=1,llm
+               DO j=1,jjp1
+                  DO i=1,iip1
+                    q(i,j,l,iq)=val
+                  ENDDO
+               ENDDO
+             ENDDO
+             write(*,*) 'SURFACE value of tracer ',trim(tname(iq)),
+     &                   ' ? (kg/m2)'
+             read(*,*) val
+             DO ig=1,ngridmx
+                 qsurf(ig,iq)=val
+             ENDDO
+
+c       t=profile : initialize temperature with a given profile
+c       -------------------------------------------------------
+        else if (trim(modif) .eq. 't=profile') then
+             write(*,*) 'Temperature profile from ASCII file'
+             write(*,*) "'profile.in' e.g. 1D output"
+             write(*,*) "(one value per line in file; starting with"
+             write(*,*) "surface value, the 1st atmospheric layer"
+             write(*,*) "followed by 2nd, etc. up to top of atmosphere)"
+             txt="profile.in"
+             open(33,file=trim(txt),status='old',form='formatted',
+     &            iostat=ierr)
+             if (ierr.eq.0) then
+               ! OK, found file 'profile_...', load the profile
+               do l=1,llm+1
+                 read(33,*,iostat=ierr) profile(l)
+                 write(*,*) profile(l)
+                 if (ierr.ne.0) then ! something went wrong
+                   exit ! quit loop
+                 endif
+               enddo
+               if (ierr.eq.0) then
+                 tsurf(1:ngridmx)=profile(1)
+                 tsoil(1:ngridmx,1:nsoilmx)=profile(1)
+                 do l=1,llm
+                   Tset(1:iip1,1:jjp1,l)=profile(l+1)
+                   flagtset=.true.
+                 enddo
+                 ucov(1:iip1,1:jjp1,1:llm)=0.
+                 vcov(1:iip1,1:jjm,1:llm)=0.
+                 q2(1:ngridmx,1:llm+1)=0.
+               else
+                 write(*,*)'problem reading file ',trim(txt),' !'
+                 write(*,*)'No modifications to temperature'
+               endif
+             else
+               write(*,*)'Could not find file ',trim(txt),' !'
+               write(*,*)'No modifications to temperature'
+             endif
+
+c       q=profile : initialize tracer with a given profile
+c       --------------------------------------------------
+        else if (trim(modif) .eq. 'q=profile') then
+             write(*,*) 'Tracer profile will be sought in ASCII file'
+             write(*,*) "'profile_tracer' where 'tracer' is tracer name"
+             write(*,*) "(one value per line in file; starting with"
+             write(*,*) "surface value, the 1st atmospheric layer"
+             write(*,*) "followed by 2nd, etc. up to top of atmosphere)"
+             write(*,*) 'Which tracer do you want to set?'
+             do iq=1,nqtot
+               write(*,*)iq,' : ',trim(tname(iq))
+             enddo
+             write(*,*) '(choose between 1 and ',nqtot,')'
+             read(*,*) iq 
+             if ((iq.lt.1).or.(iq.gt.nqtot)) then
+               ! wrong value for iq, go back to menu
+               write(*,*) "wrong input value:",iq
+               cycle
+             endif
+             ! look for input file 'profile_tracer'
+             txt="profile_"//trim(tname(iq))
+             open(41,file=trim(txt),status='old',form='formatted',
+     &            iostat=ierr)
+             if (ierr.eq.0) then
+               ! OK, found file 'profile_...', load the profile
+               do l=1,llm+1
+                 read(41,*,iostat=ierr) profile(l)
+                 if (ierr.ne.0) then ! something went wrong
+                   exit ! quit loop
+                 endif
+               enddo
+               if (ierr.eq.0) then
+                 ! initialize tracer values
+                 qsurf(:,iq)=profile(1)
+                 do l=1,llm
+                   q(:,:,l,iq)=profile(l+1)
+                 enddo
+                 write(*,*)'OK, tracer ',trim(tname(iq)),' initialized '
+     &                    ,'using values from file ',trim(txt)
+               else
+                 write(*,*)'problem reading file ',trim(txt),' !'
+                 write(*,*)'No modifications to tracer ',trim(tname(iq))
+               endif
+             else
+               write(*,*)'Could not find file ',trim(txt),' !'
+               write(*,*)'No modifications to tracer ',trim(tname(iq))
+             endif
+             
+
+c      wetstart : wet atmosphere with a north to south gradient
+c      --------------------------------------------------------
+       else if (trim(modif) .eq. 'wetstart') then
+        ! check that there is indeed a water vapor tracer
+        if (igcm_h2o_vap.eq.0) then
+          write(*,*) "No water vapour tracer! Can't use this option"
+          stop
+        endif
+          DO l=1,llm
+            DO j=1,jjp1
+              DO i=1,iip1
+                q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi
+              ENDDO
+            ENDDO
+          ENDDO
+
+         write(*,*) 'Water mass mixing ratio at north pole='
+     *               ,q(1,1,1,igcm_h2o_vap)
+         write(*,*) '---------------------------south pole='
+     *               ,q(1,jjp1,1,igcm_h2o_vap)
+
+c      noglacier : remove tropical water ice (to initialize high res sim)
+c      --------------------------------------------------
+        else if (trim(modif) .eq. 'noglacier') then
+           if (igcm_h2o_ice.eq.0) then
+             write(*,*) "No water ice tracer! Can't use this option"
+             stop
+           endif
+           do ig=1,ngridmx
+             j=(ig-2)/iim +2
+              if(ig.eq.1) j=1
+              write(*,*) 'OK: remove surface ice for |lat|<45'
+              if (abs(rlatu(j)*180./pi).lt.45.) then
+                   qsurf(ig,igcm_h2o_ice)=0.
+              end if
+           end do
+
+
+c      watercapn : H20 ice on permanent northern cap
+c      --------------------------------------------------
+        else if (trim(modif) .eq. 'watercapn') then
+           if (igcm_h2o_ice.eq.0) then
+             write(*,*) "No water ice tracer! Can't use this option"
+             stop
+           endif
+
+           DO j=1,jjp1        
+              DO i=1,iim
+                 ig=1+(j-2)*iim +i
+                 if(j.eq.1) ig=1
+                 if(j.eq.jjp1) ig=ngridmx
+
+                 if (rlatu(j)*180./pi.gt.80.) then
+                    qsurf(ig,igcm_h2o_ice)=3.4e3
+                    !do isoil=1,nsoilmx
+                    !   ith(i,j,isoil) = 18000. ! thermal inertia
+                    !enddo
+                   write(*,*)'     ==> Ice mesh North boundary (deg)= ',
+     &                   rlatv(j-1)*180./pi
+                 end if
+              ENDDO
+           ENDDO
+           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
+
+c$$$           do ig=1,ngridmx
+c$$$             j=(ig-2)/iim +2
+c$$$              if(ig.eq.1) j=1
+c$$$              if (rlatu(j)*180./pi.gt.80.) then
+c$$$
+c$$$                   qsurf(ig,igcm_h2o_ice)=1.e5
+c$$$                   qsurf(ig,igcm_h2o_vap)=0.0!1.e5
+c$$$
+c$$$                   write(*,*) 'ig=',ig,'    H2O ice mass (kg/m2)= ',
+c$$$     &              qsurf(ig,igcm_h2o_ice)
+c$$$
+c$$$                   write(*,*)'     ==> Ice mesh South boundary (deg)= ',
+c$$$     &              rlatv(j)*180./pi
+c$$$                end if
+c$$$           enddo
+
+c      watercaps : H20 ice on permanent southern cap
+c      -------------------------------------------------
+        else if (trim(modif) .eq. 'watercaps') then
+           if (igcm_h2o_ice.eq.0) then
+              write(*,*) "No water ice tracer! Can't use this option"
+              stop
+           endif
+
+           DO j=1,jjp1        
+              DO i=1,iim
+                 ig=1+(j-2)*iim +i
+                 if(j.eq.1) ig=1
+                 if(j.eq.jjp1) ig=ngridmx
+
+                 if (rlatu(j)*180./pi.lt.-80.) then
+                    qsurf(ig,igcm_h2o_ice)=3.4e3
+                    !do isoil=1,nsoilmx
+                    !   ith(i,j,isoil) = 18000. ! thermal inertia
+                    !enddo
+                   write(*,*)'     ==> Ice mesh South boundary (deg)= ',
+     &                   rlatv(j-1)*180./pi
+                 end if
+              ENDDO
+           ENDDO
+           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
+
+c$$$           do ig=1,ngridmx
+c$$$               j=(ig-2)/iim +2
+c$$$               if(ig.eq.1) j=1
+c$$$               if (rlatu(j)*180./pi.lt.-80.) then
+c$$$                  qsurf(ig,igcm_h2o_ice)=1.e5
+c$$$                  qsurf(ig,igcm_h2o_vap)=0.0 !1.e5
+c$$$
+c$$$                  write(*,*) 'ig=',ig,'   H2O ice mass (kg/m2)= ',
+c$$$     &                 qsurf(ig,igcm_h2o_ice)
+c$$$                  write(*,*)'     ==> Ice mesh North boundary (deg)= ',
+c$$$     &                 rlatv(j-1)*180./pi
+c$$$               end if
+c$$$           enddo
+
+
+c       noacglac : H2O ice across highest terrain
+c       --------------------------------------------
+        else if (trim(modif) .eq. 'noacglac') then
+           if (igcm_h2o_ice.eq.0) then
+             write(*,*) "No water ice tracer! Can't use this option"
+             stop
+           endif
+          DO j=1,jjp1        
+             DO i=1,iim
+                ig=1+(j-2)*iim +i
+                if(j.eq.1) ig=1
+                if(j.eq.jjp1) ig=ngridmx
+
+                if(phis(i,j).gt.1000.*g)then
+                    alb(i,j) = 0.5 ! snow value
+                    do isoil=1,nsoilmx
+                       ith(i,j,isoil) = 18000. ! thermal inertia
+                       ! this leads to rnat set to 'ocean' in physiq.F90
+                       ! actually no, because it is soil not surface
+                    enddo
+                endif
+             ENDDO
+          ENDDO
+          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+
+
+
+c       oborealis : H2O oceans across Vastitas Borealis
+c       -----------------------------------------------
+        else if (trim(modif) .eq. 'oborealis') then
+           if (igcm_h2o_ice.eq.0) then
+             write(*,*) "No water ice tracer! Can't use this option"
+             stop
+           endif
+          DO j=1,jjp1        
+             DO i=1,iim
+                ig=1+(j-2)*iim +i
+                if(j.eq.1) ig=1
+                if(j.eq.jjp1) ig=ngridmx
+
+                if(phis(i,j).lt.-4000.*g)then
+!                if( (phis(i,j).lt.-4000.*g)
+!     &               .and. (rlatu(j)*180./pi.lt.0.) )then ! south hemisphere only
+!                    phis(i,j)=-8000.0*g ! proper ocean
+                    phis(i,j)=-4000.0*g ! scanty ocean
+
+                    alb(i,j) = 0.07 ! oceanic value
+                    do isoil=1,nsoilmx
+                       ith(i,j,isoil) = 18000. ! thermal inertia
+                       ! this leads to rnat set to 'ocean' in physiq.F90
+                       ! actually no, because it is soil not surface
+                    enddo
+                endif
+             ENDDO
+          ENDDO
+          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+
+c       iborealis : H2O ice in Northern plains
+c       --------------------------------------
+        else if (trim(modif) .eq. 'iborealis') then
+           if (igcm_h2o_ice.eq.0) then
+             write(*,*) "No water ice tracer! Can't use this option"
+             stop
+           endif
+          DO j=1,jjp1        
+             DO i=1,iim
+                ig=1+(j-2)*iim +i
+                if(j.eq.1) ig=1
+                if(j.eq.jjp1) ig=ngridmx
+
+                if(phis(i,j).lt.-4000.*g)then
+                   !qsurf(ig,igcm_h2o_ice)=1.e3
+                   qsurf(ig,igcm_h2o_ice)=241.4 ! to make total 33 kg m^-2
+                endif
+             ENDDO
+          ENDDO
+          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+
+
+c       oceanball : H2O liquid everywhere
+c       ----------------------------
+        else if (trim(modif) .eq. 'oceanball') then
+           if (igcm_h2o_ice.eq.0) then
+             write(*,*) "No water ice tracer! Can't use this option"
+             stop
+           endif
+          DO j=1,jjp1        
+             DO i=1,iim
+                ig=1+(j-2)*iim +i
+                if(j.eq.1) ig=1
+                if(j.eq.jjp1) ig=ngridmx
+
+                qsurf(ig,igcm_h2o_vap)=0.0    ! for ocean, this is infinite source
+                qsurf(ig,igcm_h2o_ice)=0.0
+                alb(i,j) = 0.07 ! ocean value
+
+                do isoil=1,nsoilmx
+                   ith(i,j,isoil) = 18000. ! thermal inertia
+                   !ith(i,j,isoil) = 50. ! extremely low for test
+                   ! this leads to rnat set to 'ocean' in physiq.F90
+                enddo
+
+             ENDDO
+          ENDDO
+          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
+
+c       iceball : H2O ice everywhere
+c       ----------------------------
+        else if (trim(modif) .eq. 'iceball') then
+           if (igcm_h2o_ice.eq.0) then
+             write(*,*) "No water ice tracer! Can't use this option"
+             stop
+           endif
+          DO j=1,jjp1        
+             DO i=1,iim
+                ig=1+(j-2)*iim +i
+                if(j.eq.1) ig=1
+                if(j.eq.jjp1) ig=ngridmx
+
+                qsurf(ig,igcm_h2o_vap)=-50.    ! for ocean, this is infinite source
+                qsurf(ig,igcm_h2o_ice)=50.     ! == to 0.05 m of oceanic ice
+                alb(i,j) = 0.6 ! ice albedo value
+
+                do isoil=1,nsoilmx
+                   !ith(i,j,isoil) = 18000. ! thermal inertia
+                   ! this leads to rnat set to 'ocean' in physiq.F90
+                enddo
+
+             ENDDO
+          ENDDO
+          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
+
+c       supercontinent : H2O ice everywhere
+c       ----------------------------
+        else if (trim(modif) .eq. 'supercontinent') then
+          write(*,*) 'Minimum longitude (-180,180)'
+          read(*,*) val
+          write(*,*) 'Maximum longitude (-180,180)'
+          read(*,*) val2
+          write(*,*) 'Minimum latitude (-90,90)'
+          read(*,*) val3
+          write(*,*) 'Maximum latitude (-90,90)'
+          read(*,*) val4
+
+          do j=1,jjp1
+            do i=1,iip1
+              ig=1+(j-2)*iim +i
+              if(j.eq.1) ig=1
+              if(j.eq.jjp1) ig=ngridmx
+
+c             Supercontinent:
+              if (((rlatu(j)*180./pi.le.val4).and.
+     &            (rlatu(j)*180./pi.ge.val3).and.
+     &            (rlonv(i)*180./pi.le.val2).and.
+     &            (rlonv(i)*180./pi.ge.val))) then
+
+                rnat(ig)=1.
+                alb(i,j) = 0.3
+                do isoil=1,nsoilmx
+                  ith(i,j,isoil) = 2000.
+                enddo
+c             Ocean:
+              else
+                rnat(ig)=0.
+                alb(i,j) = 0.07
+                do isoil=1,nsoilmx
+                  ith(i,j,isoil) = 18000.
+                enddo
+              end if
+
+            enddo
+          enddo
+          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
+
+c       isotherm : Isothermal temperatures and no winds
+c       -----------------------------------------------
+        else if (trim(modif) .eq. 'isotherm') then
+
+          write(*,*)'Isothermal temperature of the atmosphere, 
+     &           surface and subsurface'
+          write(*,*) 'Value of this temperature ? :'
+ 203      read(*,*,iostat=ierr) Tiso
+          if(ierr.ne.0) goto 203
+
+          tsurf(1:ngridmx)=Tiso
+          
+          tsoil(1:ngridmx,1:nsoilmx)=Tiso
+          
+          Tset(1:iip1,1:jjp1,1:llm)=Tiso
+          flagtset=.true.
+
+          t(1:iip1,1:jjp1,1:llm)=Tiso
+          !! otherwise hydrost. integrations below
+          !! use the wrong temperature
+          !! -- NB: Tset might be useless
+        
+          ucov(1:iip1,1:jjp1,1:llm)=0
+          vcov(1:iip1,1:jjm,1:llm)=0
+          q2(1:ngridmx,1:llm+1)=0
+
+c       radequi : Radiative equilibrium profile of temperatures and no winds
+c       --------------------------------------------------------------------
+        else if (trim(modif) .eq. 'radequi') then
+
+          write(*,*)'radiative equilibrium temperature profile'       
+
+          do ig=1, ngridmx
+             teque= 335.0-60.0*sin(latfi(ig))*sin(latfi(ig))-
+     &            10.0*cos(latfi(ig))*cos(latfi(ig))
+             tsurf(ig) = MAX(220.0,teque)
+          end do
+          do l=2,nsoilmx
+             do ig=1, ngridmx
+                tsoil(ig,l) = tsurf(ig)
+             end do
+          end do
+          DO j=1,jjp1
+             DO i=1,iim
+                Do l=1,llm
+                   teque=335.-60.0*sin(rlatu(j))*sin(rlatu(j))-
+     &                  10.0*cos(rlatu(j))*cos(rlatu(j))
+                   Tset(i,j,l)=MAX(220.0,teque)
+                end do
+             end do
+          end do
+          flagtset=.true.
+          ucov(1:iip1,1:jjp1,1:llm)=0
+          vcov(1:iip1,1:jjm,1:llm)=0
+          q2(1:ngridmx,1:llm+1)=0
+
+c       coldstart : T set 1K above CO2 frost point and no winds
+c       ------------------------------------------------
+        else if (trim(modif) .eq. 'coldstart') then
+
+          write(*,*)'set temperature of the atmosphere,' 
+     &,'surface and subsurface how many degrees above CO2 frost point?'
+ 204      read(*,*,iostat=ierr) Tabove
+          if(ierr.ne.0) goto 204
+
+            DO j=1,jjp1
+             DO i=1,iim
+                ig=1+(j-2)*iim +i
+                if(j.eq.1) ig=1
+                if(j.eq.jjp1) ig=ngridmx
+            tsurf(ig) = (-3167.8)/(log(.01*ps(i,j))-23.23)+Tabove
+             END DO
+            END DO
+          do l=1,nsoilmx
+            do ig=1, ngridmx
+              tsoil(ig,l) = tsurf(ig)
+            end do
+          end do
+          DO j=1,jjp1
+           DO i=1,iim
+            Do l=1,llm
+               pp = aps(l) +bps(l)*ps(i,j) 
+               Tset(i,j,l)=(-3167.8)/(log(.01*pp)-23.23)+Tabove
+            end do
+           end do
+          end do
+
+          flagtset=.true.
+          ucov(1:iip1,1:jjp1,1:llm)=0
+          vcov(1:iip1,1:jjm,1:llm)=0
+          q2(1:ngridmx,1:llm+1)=0
+
+
+c       co2ice=0 : remove CO2 polar ice caps'
+c       ------------------------------------------------
+        else if (trim(modif) .eq. 'co2ice=0') then
+         ! check that there is indeed a co2_ice tracer ...
+          if (igcm_co2_ice.ne.0) then
+           do ig=1,ngridmx
+              !co2ice(ig)=0
+              qsurf(ig,igcm_co2_ice)=0
+              emis(ig)=emis(ngridmx/2)
+           end do
+          else
+            write(*,*) "Can't remove CO2 ice!! (no co2_ice tracer)"
+          endif
+        
+!       therm_ini_s: (re)-set soil thermal inertia to reference surface values
+!       ----------------------------------------------------------------------
+
+        else if (trim(modif) .eq. 'therm_ini_s') then
+!          write(*,*)"surfithfi(1):",surfithfi(1)
+          do isoil=1,nsoilmx
+            inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
+          enddo
+          write(*,*)'OK: Soil thermal inertia has been reset to referenc
+     &e surface values'
+!          write(*,*)"inertiedat(1,1):",inertiedat(1,1)
+          ithfi(:,:)=inertiedat(:,:)
+         ! recast ithfi() onto ith()
+         call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
+! Check:
+!         do i=1,iip1
+!           do j=1,jjp1
+!             do isoil=1,nsoilmx
+!               write(77,*) i,j,isoil,"  ",ith(i,j,isoil)
+!             enddo
+!           enddo
+!         enddo
+
+
+
+c       slab_ocean_initialisation
+c       ------------------------------------------------
+        else if (trim(modif) .eq. 'slab_ocean_0') then
+        write(*,*)'OK: initialisation of slab ocean' 
+
+      DO ig=1, ngridmx
+         rnat(ig)=1.
+         tslab(ig,1)=0.
+         tslab(ig,2)=0.
+         tsea_ice(ig)=0.
+         sea_ice(ig)=0.
+         pctsrf_sic(ig)=0.
+         
+         if(ithfi(ig,1).GT.10000.)then
+           rnat(ig)=0.
+           phisfi(ig)=0.
+           tsea_ice(ig)=273.15-1.8
+           tslab(ig,1)=tsurf(ig)
+           tslab(ig,2)=tsurf(ig)!*2./3.+(273.15-1.8)/3.
+          endif
+
+      ENDDO
+          CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis)
+
+
+
+        else
+          write(*,*) '       Unknown (misspelled?) option!!!'
+        end if ! of if (trim(modif) .eq. '...') elseif ...
+
+
+
+       enddo ! of do ! infinite loop on liste of changes
+
+ 999  continue
+
+ 
+c=======================================================================
+c   Initialisation for cloud fraction and oceanic ice (added by BC 2010)
+c=======================================================================
+      DO ig=1, ngridmx
+         totalfrac(ig)=0.5
+         DO l=1,llm
+            cloudfrac(ig,l)=0.5
+         ENDDO
+! Ehouarn, march 2012: also add some initialisation for hice
+         hice(ig)=0.0
+      ENDDO
+      
+c========================================================
+
+!      DO ig=1,ngridmx
+!         IF(tsurf(ig) .LT. 273.16-1.8) THEN
+!            hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1.
+!            hice(ig)=min(hice(ig),1.0)
+!         ENDIF
+!      ENDDO
+
+
+
+
+c=======================================================================
+c   Correct pressure on the new grid (menu 0)
+c=======================================================================
+
+
+      if ((choix_1.eq.0).and.(.not.flagps0)) then
+        r = 1000.*8.31/mugaz
+
+        do j=1,jjp1
+          do i=1,iip1
+             ps(i,j) = ps(i,j) *
+     .            exp((phisold_newgrid(i,j)-phis(i,j)) /
+     .                                  (t(i,j,1) * r))
+          end do
+        end do
+
+c periodicite de ps en longitude
+        do j=1,jjp1
+          ps(1,j) = ps(iip1,j)
+        end do
+      end if
+
+         
+c=======================================================================
+c=======================================================================
+
+c=======================================================================
+c    Initialisation de la physique / ecriture de newstartfi :
+c=======================================================================
+
+
+      CALL inifilr 
+      CALL pression(ip1jmp1, ap, bp, ps, p3d)
+
+c-----------------------------------------------------------------------
+c   Initialisation  pks:
+c-----------------------------------------------------------------------
+
+      CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
+! Calcul de la temperature potentielle teta
+
+      if (flagtset) then
+          DO l=1,llm
+             DO j=1,jjp1
+                DO i=1,iim
+                   teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l)
+                ENDDO
+                teta (iip1,j,l)= teta (1,j,l)
+             ENDDO
+          ENDDO
+      else if (choix_1.eq.0) then
+         DO l=1,llm
+            DO j=1,jjp1
+               DO i=1,iim
+                  teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l)
+               ENDDO
+               teta (iip1,j,l)= teta (1,j,l)
+            ENDDO
+         ENDDO
+      endif
+
+C Calcul intermediaire
+c
+      if (choix_1.eq.0) then
+         CALL massdair( p3d, masse  )
+c
+         print *,' ALPHAX ',alphax
+
+         DO  l = 1, llm
+           DO  i    = 1, iim
+             xppn(i) = aire( i, 1   ) * masse(  i     ,  1   , l )
+             xpps(i) = aire( i,jjp1 ) * masse(  i     , jjp1 , l )
+           ENDDO
+             xpn      = SUM(xppn)/apoln
+             xps      = SUM(xpps)/apols
+           DO i   = 1, iip1
+             masse(   i   ,   1     ,  l )   = xpn
+             masse(   i   ,   jjp1  ,  l )   = xps
+           ENDDO
+         ENDDO
+      endif
+      phis(iip1,:) = phis(1,:)
+
+      itau=0
+      if (choix_1.eq.0) then
+         day_ini=int(date)
+      endif
+c
+      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
+
+      CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis ,
+     *                phi,w, pbaru,pbarv,day_ini+time )
+
+          
+      CALL dynredem0("restart.nc",day_ini,phis)
+      CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,masse,ps) 
+C
+C Ecriture etat initial physique
+C
+
+      call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
+     &              nqtot,dtphys,real(day_ini),0.0,
+     &              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
+      call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
+     &                dtphys,real(day_ini),
+     &                tsurf,tsoil,emis,q2,qsurf,
+     &                cloudfrac,totalfrac,hice,
+     &                rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
+
+c=======================================================================
+c        Formats 
+c=======================================================================
+
+   1  FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema
+     *rrage est differente de la valeur parametree iim =',i4//)
+   2  FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema
+     *rrage est differente de la valeur parametree jjm =',i4//)
+   3  FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar
+     *rage est differente de la valeur parametree llm =',i4//)
+
+      write(*,*) "newstart: All is well that ends well."
+
+      end
+
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/scal_wind.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/scal_wind.F	(revision 1644)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/scal_wind.F	(revision 1644)
@@ -0,0 +1,55 @@
+      SUBROUTINE scal_wind(xus,xvs,xu,xv)
+c=======================================================================
+c
+c
+c   Subject:
+c   ------
+c On passe  les variable xus, xvs  aux points de vent u et v (xu et xv)
+c
+c=======================================================================
+      IMPLICIT NONE
+c-----------------------------------------------------------------------
+c   Declararations:
+c   ---------------
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comgeom.h"
+
+c   Arguments:
+c   ----------
+
+      REAL xu(iip1,jjp1,llm),xv(iip1,jjm,llm)
+      REAL xus(iip1,jjp1,llm), xvs (iip1,jjp1,llm)
+
+c   Local:
+c   ------
+
+      INTEGER i,j,l
+
+c-----------------------------------------------------------------------
+
+c   transport zonal:
+c   ----------------
+      DO l=1,llm
+        Do j=1,jjp1
+	      DO i=1,iim
+            xu(i,j,l)=0.5*(xus(i,j,l)+xus(i+1,j,l))
+	      ENDDO
+          xu(iip1,j,l)=xu(1,j,l)
+	    ENDDO
+      ENDDO
+
+
+c   Transport meridien:
+c   -------------------
+      DO l=1,llm
+         DO j=1,jjm
+           do i=1 ,iip1
+	         xv(i,j,l)=.5*(xvs(i,j,l)+xvs(i,j+1,l))
+           end do
+	     ENDDO
+	  ENDDO
+
+      RETURN
+      END
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/start2archive.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/start2archive.F	(revision 1644)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/start2archive.F	(revision 1644)
@@ -0,0 +1,544 @@
+c=======================================================================
+      PROGRAM start2archive
+c=======================================================================
+c
+c
+c   Date:    01/1997
+c   ----
+c
+c
+c   Objet:   Passage des  fichiers netcdf d'etat initial "start" et
+c   -----    "startfi" a un fichier netcdf unique "start_archive" 
+c
+c  "start_archive" est une banque d'etats initiaux:
+c  On peut stocker plusieurs etats initiaux dans un meme fichier "start_archive"
+c    (Veiller dans ce cas avoir un day_ini different pour chacun des start)
+c 
+c
+c
+c=======================================================================
+
+      use infotrac, only: infotrac_init, nqtot, tname
+      USE comsoil_h
+      
+!      USE comgeomfi_h, ONLY: lati, long, area
+!      use control_mod
+!      use comgeomphy, only: initcomgeomphy
+      use slab_ice_h, only: noceanmx
+! to use  'getin'
+      USE ioipsl_getincom
+      USE planete_mod, only: year_day
+      USE mod_const_mpi, ONLY: COMM_LMDZ
+      USE control_mod, only: planet_type
+      USE callkeys_mod, ONLY: ok_slab_ocean
+      use filtreg_mod, only: inifilr
+      USE comvert_mod, ONLY: ap,bp
+      USE comconst_mod, ONLY: daysec,dtphys,rad,g,r,cpp
+      USE temps_mod, ONLY: day_ini
+      USE iniphysiq_mod, ONLY: iniphysiq
+      implicit none
+
+      include "dimensions.h"
+      integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm) 
+      include "paramet.h"
+      include "comdissip.h"
+      include "comgeom.h"
+!#include "control.h"
+
+!#include "dimphys.h"
+!#include "planete.h"
+!#include"advtrac.h"
+      include "netcdf.inc"
+c-----------------------------------------------------------------------
+c   Declarations
+c-----------------------------------------------------------------------
+
+c variables dynamiques du GCM
+c -----------------------------
+      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
+      REAL teta(ip1jmp1,llm)                    ! temperature potentielle 
+      REAL,ALLOCATABLE :: q(:,:,:)   ! champs advectes
+      REAL pks(ip1jmp1)                      ! exner (f pour filtre)
+      REAL pk(ip1jmp1,llm)
+      REAL pkf(ip1jmp1,llm)
+      REAL beta(iip1,jjp1,llm)
+      REAL phis(ip1jmp1)                     ! geopotentiel au sol
+      REAL masse(ip1jmp1,llm)                ! masse de l'atmosphere
+      REAL ps(ip1jmp1)                       ! pression au sol
+      REAL p3d(iip1, jjp1, llm+1)            ! pression aux interfaces
+      
+c Variable Physiques (grille physique)
+c ------------------------------------
+      REAL tsurf(ngridmx)        ! Surface temperature
+      REAL,ALLOCATABLE :: tsoil(:,:) ! Soil temperature
+      REAL co2ice(ngridmx)        ! CO2 ice layer
+      REAL q2(ngridmx,llm+1)
+      REAL,ALLOCATABLE :: qsurf(:,:)
+      REAL emis(ngridmx)
+      INTEGER start,length
+      PARAMETER (length = 100)
+      REAL tab_cntrl_fi(length) ! tableau des parametres de startfi
+      REAL tab_cntrl_dyn(length) ! tableau des parametres de start
+      INTEGER*4 day_ini_fi
+
+!     added by FF for cloud fraction setup
+      REAL hice(ngridmx)
+      REAL cloudfrac(ngridmx,llm),totalcloudfrac(ngridmx)
+
+!     added by BC for slab ocean
+      REAL rnat(ngridmx),pctsrf_sic(ngridmx),sea_ice(ngridmx)
+      REAL tslab(ngridmx,noceanmx),tsea_ice(ngridmx)
+
+
+c Variable naturelle / grille scalaire
+c ------------------------------------
+      REAL T(ip1jmp1,llm),us(ip1jmp1,llm),vs(ip1jmp1,llm)
+      REAL tsurfS(ip1jmp1)
+      REAL,ALLOCATABLE :: tsoilS(:,:)
+      REAL,ALLOCATABLE :: ithS(:,:) ! Soil Thermal Inertia
+      REAL co2iceS(ip1jmp1)
+      REAL q2S(ip1jmp1,llm+1)
+      REAL,ALLOCATABLE :: qsurfS(:,:)
+      REAL emisS(ip1jmp1)
+
+!     added by FF for cloud fraction setup
+      REAL hiceS(ip1jmp1)
+      REAL cloudfracS(ip1jmp1,llm),totalcloudfracS(ip1jmp1)
+
+!     added by BC for slab ocean
+      REAL rnatS(ip1jmp1),pctsrf_sicS(ip1jmp1),sea_iceS(ip1jmp1)
+      REAL tslabS(ip1jmp1,noceanmx),tsea_iceS(ip1jmp1)
+
+
+c Variables intermediaires : vent naturel, mais pas coord scalaire
+c----------------------------------------------------------------
+      REAL vn(ip1jm,llm),un(ip1jmp1,llm)
+
+c Autres  variables
+c -----------------
+      LOGICAL startdrs
+      INTEGER Lmodif
+
+      REAL ptotal, co2icetotal
+      REAL timedyn,timefi !fraction du jour dans start, startfi
+      REAL date
+
+      CHARACTER*2 str2
+      CHARACTER*80 fichier 
+      data  fichier /'startfi'/
+
+      INTEGER ij, l,i,j,isoil,iq
+      character*80      fichnom
+      integer :: ierr,ntime
+      integer :: nq,numvanle
+      character(len=30) :: txt ! to store some text
+
+c Netcdf
+c-------
+      integer varid,dimid,timelen 
+      INTEGER nid,nid1
+
+c-----------------------------------------------------------------------
+c   Initialisations 
+c-----------------------------------------------------------------------
+
+      CALL defrun_new(99, .TRUE. )
+
+      planet_type="titan"
+
+c=======================================================================
+c Lecture des donnees
+c=======================================================================
+! Load tracer number and names:
+      call infotrac_init
+
+! allocate arrays:
+      allocate(q(ip1jmp1,llm,nqtot))
+      allocate(qsurf(ngridmx,nqtot))
+      allocate(qsurfS(ip1jmp1,nqtot))
+! other array allocations:
+!      call ini_comsoil_h(ngridmx) ! done via iniphysiq
+
+      fichnom = 'start.nc'
+      CALL dynetat0(fichnom,vcov,ucov,teta,q,masse,
+     .       ps,phis,timedyn)
+
+! load 'controle' array from dynamics start file
+
+       ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
+       IF (ierr.NE.NF_NOERR) THEN
+         write(6,*)' Pb d''ouverture du fichier'//trim(fichnom)
+        CALL ABORT
+       ENDIF
+                                                
+      ierr = NF_INQ_VARID (nid1, "controle", varid)
+      IF (ierr .NE. NF_NOERR) THEN
+       PRINT*, "start2archive: Le champ <controle> est absent"
+       CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+       ierr = NF_GET_VAR_DOUBLE(nid1, varid, tab_cntrl_dyn)
+#else
+      ierr = NF_GET_VAR_REAL(nid1, varid, tab_cntrl_dyn)
+#endif
+       IF (ierr .NE. NF_NOERR) THEN
+          PRINT*, "start2archive: Lecture echoue pour <controle>"
+          CALL abort
+       ENDIF
+
+      ierr = NF_CLOSE(nid1)
+
+! Get value of the "subsurface_layers" dimension from physics start file
+      fichnom = 'startfi.nc'
+      ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
+       IF (ierr.NE.NF_NOERR) THEN
+         write(6,*)' Pb d''ouverture du fichier'//trim(fichnom)
+        CALL ABORT
+       ENDIF
+      ierr = NF_INQ_DIMID(nid1,"subsurface_layers",varid)
+      IF (ierr .NE. NF_NOERR) THEN
+       PRINT*, "start2archive: No subsurface_layers dimension!!"
+       CALL abort
+      ENDIF
+      ierr = NF_INQ_DIMLEN(nid1,varid,nsoilmx)
+      IF (ierr .NE. NF_NOERR) THEN
+       PRINT*, "start2archive: Failed reading subsurface_layers value!!"
+       CALL abort
+      ENDIF
+      ierr = NF_CLOSE(nid1)
+      
+      ! allocate arrays of nsoilmx size
+      allocate(tsoil(ngridmx,nsoilmx))
+      allocate(tsoilS(ip1jmp1,nsoilmx))
+      allocate(ithS(ip1jmp1,nsoilmx))
+
+c-----------------------------------------------------------------------
+c   Initialisations 
+c-----------------------------------------------------------------------
+
+      CALL defrun_new(99, .FALSE. )
+      call iniconst
+      call inigeom
+      call inifilr
+
+! Initialize the physics
+         CALL iniphysiq(iim,jjm,llm,
+     &                  (jjm-1)*iim+2,comm_lmdz,
+     &                  daysec,day_ini,dtphys,
+     &                  rlatu,rlatv,rlonu,rlonv,
+     &                  aire,cu,cv,rad,g,r,cpp,
+     &                  1)
+
+      fichnom = 'startfi.nc'
+      Lmodif=0
+
+
+      CALL phyetat0 (ngridmx,llm,fichnom,0,Lmodif,nsoilmx,nqtot,
+     .      day_ini_fi,timefi,
+     .      tsurf,tsoil,emis,q2,qsurf,
+!       change FF 05/2011
+     .       cloudfrac,totalcloudfrac,hice,
+!       change BC 05/2014
+     .       rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
+
+
+
+
+! load 'controle' array from physics start file
+
+       ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
+       IF (ierr.NE.NF_NOERR) THEN
+         write(6,*)' Pb d''ouverture du fichier'//trim(fichnom)
+        CALL ABORT
+       ENDIF
+                                                
+      ierr = NF_INQ_VARID (nid1, "controle", varid)
+      IF (ierr .NE. NF_NOERR) THEN
+       PRINT*, "start2archive: Le champ <controle> est absent"
+       CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+       ierr = NF_GET_VAR_DOUBLE(nid1, varid, tab_cntrl_fi)
+#else
+      ierr = NF_GET_VAR_REAL(nid1, varid, tab_cntrl_fi)
+#endif
+       IF (ierr .NE. NF_NOERR) THEN
+          PRINT*, "start2archive: Lecture echoue pour <controle>"
+          CALL abort
+       ENDIF
+
+      ierr = NF_CLOSE(nid1)
+
+
+c-----------------------------------------------------------------------
+c Controle de la synchro
+c-----------------------------------------------------------------------
+!mars a voir      if ((day_ini_fi.ne.day_ini).or.(abs(timefi-timedyn).gt.1.e-10)) 
+      if ((day_ini_fi.ne.day_ini)) 
+     &  stop ' Probleme de Synchro entre start et startfi !!!'
+
+
+c *****************************************************************
+c    Option : Reinitialisation des dates dans la premieres annees :
+       do while (day_ini.ge.year_day)
+          day_ini=day_ini-year_day
+       enddo
+c *****************************************************************
+
+      CALL pression(ip1jmp1, ap, bp, ps, p3d)
+      call exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
+
+c=======================================================================
+c Transformation EN VARIABLE NATURELLE / GRILLE SCALAIRE si necessaire
+c=======================================================================
+c  Les variables modeles dependent de la resolution. Il faut donc
+c  eliminer les facteurs responsables de cette dependance
+c  (pour utiliser newstart)
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c Vent   (depend de la resolution horizontale) 
+c-----------------------------------------------------------------------
+c
+c ucov --> un  et  vcov --> vn
+c un --> us  et   vn --> vs
+c
+c-----------------------------------------------------------------------
+
+      call covnat(llm,ucov, vcov, un, vn) 
+      call wind_scal(un,vn,us,vs) 
+
+c-----------------------------------------------------------------------
+c Temperature  (depend de la resolution verticale => de "sigma.def")
+c-----------------------------------------------------------------------
+c
+c h --> T
+c
+c-----------------------------------------------------------------------
+
+      DO l=1,llm
+         DO ij=1,ip1jmp1
+            T(ij,l)=teta(ij,l)*pk(ij,l)/cpp !mars deduit de l'equation dans newstart
+         ENDDO
+      ENDDO
+
+c-----------------------------------------------------------------------
+c Variable physique 
+c-----------------------------------------------------------------------
+c
+c tsurf --> tsurfS
+c co2ice --> co2iceS
+c tsoil --> tsoilS
+c emis --> emisS
+c q2 --> q2S
+c qsurf --> qsurfS
+c
+c-----------------------------------------------------------------------
+
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,tsurf,tsurfS)
+!      call gr_fi_dyn(1,ngridmx,iip1,jjp1,co2ice,co2iceS)
+      call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,tsoil,tsoilS)
+      ! Note: thermal inertia "inertiedat" is in comsoil.h
+      call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,inertiedat,ithS)
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,emis,emisS)
+      call gr_fi_dyn(llm+1,ngridmx,iip1,jjp1,q2,q2S)
+      call gr_fi_dyn(nqtot,ngridmx,iip1,jjp1,qsurf,qsurfS)
+      call gr_fi_dyn(llm,ngridmx,iip1,jjp1,cloudfrac,cloudfracS)
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,hice,hiceS)
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,totalcloudfrac,totalcloudfracS)
+
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,rnat,rnatS)
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,pctsrf_sic,pctsrf_sicS)
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,tsea_ice,tsea_iceS)
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,sea_ice,sea_iceS)
+      call gr_fi_dyn(noceanmx,ngridmx,iip1,jjp1,tslab,tslabS)
+
+c=======================================================================
+c Info pour controler
+c=======================================================================
+
+      ptotal =  0.
+      co2icetotal = 0.
+      DO j=1,jjp1
+         DO i=1,iim
+           ptotal=ptotal+aire(i+(iim+1)*(j-1))*ps(i+(iim+1)*(j-1))/g
+!           co2icetotal = co2icetotal + 
+!     &            co2iceS(i+(iim+1)*(j-1))*aire(i+(iim+1)*(j-1))
+         ENDDO
+      ENDDO
+      write(*,*)'Ancienne grille : masse de l''atm :',ptotal
+!      write(*,*)'Ancienne grille : masse de la glace CO2 :',co2icetotal
+
+c-----------------------------------------------------------------------
+c Passage de "ptotal" et "co2icetotal" par tab_cntrl_fi
+c-----------------------------------------------------------------------
+
+      tab_cntrl_fi(49) = ptotal
+      tab_cntrl_fi(50) = co2icetotal
+
+c=======================================================================
+c Ecriture dans le fichier  "start_archive"
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c Ouverture de "start_archive" 
+c-----------------------------------------------------------------------
+
+      ierr = NF_OPEN ('start_archive.nc', NF_WRITE,nid)
+ 
+c-----------------------------------------------------------------------
+c  si "start_archive" n'existe pas:
+c    1_ ouverture
+c    2_ creation de l'entete dynamique ("ini_archive")
+c-----------------------------------------------------------------------
+c ini_archive:
+c On met dans l'entete le tab_cntrl dynamique (1 a 16) 
+c  On y ajoute les valeurs du tab_cntrl_fi (a partir de 51)
+c  En plus les deux valeurs ptotal et co2icetotal (99 et 100)
+c-----------------------------------------------------------------------
+
+      if (ierr.ne.NF_NOERR) then
+         write(*,*)'OK, Could not open file "start_archive.nc"'
+         write(*,*)'So let s create a new "start_archive"'
+         ierr = NF_CREATE('start_archive.nc', NF_CLOBBER, nid)
+         call ini_archive(nid,day_ini,phis,ithS,tab_cntrl_fi,
+     &                                          tab_cntrl_dyn)
+      endif
+
+c-----------------------------------------------------------------------
+c Ecriture de la coordonnee temps (date en jours)
+c-----------------------------------------------------------------------
+
+      date = day_ini
+      ierr= NF_INQ_VARID(nid,"Time",varid)
+      ierr= NF_INQ_DIMID(nid,"Time",dimid)
+      ierr= NF_INQ_DIMLEN(nid,dimid,timelen)
+      ntime=timelen+1
+
+      write(*,*) "******************"
+      write(*,*) "ntime",ntime
+      write(*,*) "******************"
+#ifdef NC_DOUBLE
+      ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
+#else
+      ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
+#endif
+      if (ierr.ne.NF_NOERR) then
+         write(*,*) "time matter ",NF_STRERROR(ierr)
+         stop
+      endif
+
+c-----------------------------------------------------------------------
+c Ecriture des champs  (co2ice,emis,ps,Tsurf,T,u,v,q2,q,qsurf)
+c-----------------------------------------------------------------------
+c ATTENTION: q2 a une couche de plus!!!!
+c    Pour creer un fichier netcdf lisible par grads,
+c    On passe donc une des couches de q2 a part
+c    comme une variable 2D (la couche au sol: "q2surf")
+c    Les lmm autres couches sont nommees "q2atm" (3D) 
+c-----------------------------------------------------------------------
+
+!      call write_archive(nid,ntime,'co2ice','couche de glace co2',
+!     &  'kg/m2',2,co2iceS)
+      call write_archive(nid,ntime,'emis','grd emis',' ',2,emisS)
+      call write_archive(nid,ntime,'ps','Psurf','Pa',2,ps)
+      call write_archive(nid,ntime,'tsurf','surf T','K',2,tsurfS)
+      call write_archive(nid,ntime,'temp','temperature','K',3,t)
+      call write_archive(nid,ntime,'u','Vent zonal','m.s-1',3,us)
+      call write_archive(nid,ntime,'v','Vent merid','m.s-1',3,vs)
+      call write_archive(nid,ntime,'q2surf','wind variance','m2.s-2',2,
+     .              q2S)
+      call write_archive(nid,ntime,'q2atm','wind variance','m2.s-2',3,
+     .              q2S(1,2))
+
+c-----------------------------------------------------------------------
+c Ecriture du champs  q  ( q[1,nqtot] )
+c-----------------------------------------------------------------------
+      do iq=1,nqtot
+        call write_archive(nid,ntime,tname(iq),'tracer','kg/kg',
+     &         3,q(1,1,iq))
+      end do
+c-----------------------------------------------------------------------
+c Ecriture du champs  qsurf  ( qsurf[1,nqtot] )
+c-----------------------------------------------------------------------
+      do iq=1,nqtot
+        txt=trim(tname(iq))//"_surf"
+        call write_archive(nid,ntime,txt,'Tracer on surface',
+     &  'kg.m-2',2,qsurfS(1,iq))
+      enddo
+
+
+c-----------------------------------------------------------------------
+c Ecriture du champs  tsoil  ( Tg[1,10] )
+c-----------------------------------------------------------------------
+c "tsoil" Temperature au sol definie dans 10 couches dans le sol
+c   Les 10 couches sont lues comme 10 champs 
+c  nommees Tg[1,10]
+
+c      do isoil=1,nsoilmx
+c       write(str2,'(i2.2)') isoil
+c       call write_archive(nid,ntime,'Tg'//str2,'Ground Temperature ',
+c     .   'K',2,tsoilS(1,isoil))
+c      enddo
+
+! Write soil temperatures tsoil
+      call write_archive(nid,ntime,'tsoil','Soil temperature',
+     &     'K',-3,tsoilS)
+
+! Write soil thermal inertia
+      call write_archive(nid,ntime,'inertiedat',
+     &     'Soil thermal inertia',
+     &     'J.s-1/2.m-2.K-1',-3,ithS)
+
+! Write (0D) volumetric heat capacity (stored in comsoil.h)
+!      call write_archive(nid,ntime,'volcapa',
+!     &     'Soil volumetric heat capacity',
+!     &     'J.m-3.K-1',0,volcapa)
+! Note: no need to write volcapa, it is stored in "controle" table
+
+c-----------------------------------------------------------------------
+c Ecriture du champs  cloudfrac,hice,totalcloudfrac
+c-----------------------------------------------------------------------
+      call write_archive(nid,ntime,'hice',
+     &         'Height of oceanic ice','m',2,hiceS)
+      call write_archive(nid,ntime,'totalcloudfrac',
+     &        'Total cloud Fraction','',2,totalcloudfracS)
+      call write_archive(nid,ntime,'cloudfrac'
+     &        ,'Cloud fraction','',3,cloudfracS)
+
+c-----------------------------------------------------------------------
+c Slab ocean
+c-----------------------------------------------------------------------
+      OPEN(99,file='callphys.def',status='old',form='formatted'
+     &     ,iostat=ierr)
+      CLOSE(99)
+
+      IF(ierr.EQ.0) THEN
+
+
+         write(*,*) "Use slab-ocean ?"
+         ok_slab_ocean=.false.         ! default value
+         call getin("ok_slab_ocean",ok_slab_ocean)
+         write(*,*) "ok_slab_ocean = ",ok_slab_ocean
+
+      if(ok_slab_ocean) then
+      call write_archive(nid,ntime,'rnat'
+     &        ,'rnat','',2,rnatS)
+      call write_archive(nid,ntime,'pctsrf_sic'
+     &        ,'pctsrf_sic','',2,pctsrf_sicS)
+      call write_archive(nid,ntime,'sea_ice'
+     &        ,'sea_ice','',2,sea_iceS)
+      call write_archive(nid,ntime,'tslab'
+     &        ,'tslab','',-2,tslabS)
+      call write_archive(nid,ntime,'tsea_ice'
+     &        ,'tsea_ice','',2,tsea_iceS)
+      endif !ok_slab_ocean
+      ENDIF
+c-----------------------------------------------------------------------
+c Fin 
+c-----------------------------------------------------------------------
+      ierr=NF_CLOSE(nid)
+
+      write(*,*) "start2archive: All is well that ends well."
+
+      end 
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/wind_scal.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/wind_scal.F	(revision 1644)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/wind_scal.F	(revision 1644)
@@ -0,0 +1,55 @@
+      SUBROUTINE wind_scal(pbaru,pbarv,us,vs)
+c=======================================================================
+c
+c
+c   Subject:
+c   ------
+c   On ramene les flux de masse /vents  aux points scalaires.
+c
+c=======================================================================
+      IMPLICIT NONE
+c-----------------------------------------------------------------------
+c   Declararations:
+c   ---------------
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comgeom.h"
+
+c   Arguments:
+c   ----------
+
+      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
+      REAL us(ip1jmp1,llm), vs (ip1jmp1,llm)
+
+c   Local:
+c   ------
+
+      INTEGER ij,l
+
+c-----------------------------------------------------------------------
+
+c   transport zonal:
+c   ----------------
+      DO l=1,llm
+	 DO ij=2,ip1jmp1
+            us(ij,l)=.5*(pbaru(ij,l)+pbaru(ij-1,l))
+	 ENDDO
+      ENDDO
+      CALL SCOPY(jjp1*llm,us(iip1,1),iip1,us(1,1),iip1)
+
+
+c   Transport meridien:
+c   -------------------
+      DO l=1,llm
+         DO ij=iip2,ip1jm
+	    vs(ij,l)=.5*(pbarv(ij,l)+pbarv(ij-iip1,l))
+	 ENDDO
+	 DO ij=1,iip1
+	    vs(ij,l)=0.
+	    vs(ip1jm+ij,l)=0.
+	 ENDDO
+      ENDDO
+
+      RETURN
+      END
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/write_archive.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/write_archive.F	(revision 1644)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/write_archive.F	(revision 1644)
@@ -0,0 +1,293 @@
+c=======================================================================
+      subroutine write_archive(nid,ntime,nom,titre,unite,dim,px)
+c=======================================================================
+c
+c
+c   Date:    01/1997
+c   ----
+c
+c   Objet:   Ecriture de champs sur grille scalaire (iip1*jjp1)
+c   -----    dans un fichier DRS nomme "start_archive"
+c
+c    Il faut au prealable avoir cree un entete avec un "call ini_archive".
+c    Ces variables peuvent etre 3d (ex: temperature), 2d (ex: temperature
+c    de surface), ou 0d (pour un scalaire qui ne depend que du temps)
+c    (ex: la longitude solaire)
+c
+c
+c   Arguments: 
+c   ----------
+c
+c     Inputs:
+c     ------
+c
+c		  nid      Unite logique du fichier "start_archive"
+c         nom      nom du champ a ecrire dans le fichier "start_archive"
+c         titre    titre de la variable dans le fichier DRS "start_archive"
+c         unite    unite de la variable ....
+c         dim      dimension de la variable a ecrire
+c         px       tableau contenant la variable a ecrire
+c
+c
+c=======================================================================
+
+      use comsoil_h, only: nsoilmx
+      use slab_ice_h, only: noceanmx
+
+      implicit none
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comgeom.h"
+#include "netcdf.inc"
+
+c-----------------------------------------------------------------------
+c	Declarations   
+c-----------------------------------------------------------------------
+
+c Arguments:
+
+      INTEGER nid
+      integer ntime ! time index
+      integer dim 
+      REAL px(iip1,jjp1,llm) 
+
+      CHARACTER*(*) nom, titre, unite
+
+      integer ierr
+
+
+c local
+      integer, dimension(4) :: edges,corner,id
+      integer :: varid,i,j,l
+c-----------------------------------------------------------------------
+c      Ecriture du champs dans le fichier            (3 cas)      
+c-----------------------------------------------------------------------
+
+! For an atmospheric 3D Variable
+!--------------------------------
+        if (dim.eq.3) then
+
+!         Ecriture du champs
+
+! nom de la variable
+           ierr= NF_INQ_VARID(nid,nom,varid)
+           if (ierr /= NF_NOERR) then
+! choix du nom des coordonnees
+              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
+              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
+              ierr= NF_INQ_DIMID(nid,"altitude",id(3))
+              ierr= NF_INQ_DIMID(nid,"Time",id(4))
+
+! Creation de la variable si elle n'existait pas
+
+              write (*,*) "====================="
+              write (*,*) "creation de ",nom
+              call def_var(nid,nom,titre,unite,4,id,varid,ierr)
+
+           endif
+
+! mars s'arranger pour qu'il n'y ai plus besoin de ca
+
+c          do l=1,llm
+c             do j=1,jjp1
+c                do i=1,iip1
+c                   pxbis(i,j,l)=px(i,j,llm-l+1)
+c                enddo
+c             enddo
+c          enddo
+           corner(1)=1
+           corner(2)=1
+           corner(3)=1
+           corner(4)=ntime
+
+           edges(1)=iip1
+           edges(2)=jjp1
+           edges(3)=llm
+           edges(4)=1
+#ifdef NC_DOUBLE
+           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,px)
+#else
+           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,px)
+#endif
+
+           if (ierr.ne.NF_NOERR) then
+              write(*,*) "***** PUT_VAR matter in write_archive"
+              write(*,*) "***** with ",nom," ",nf_STRERROR(ierr)
+              call abort
+           endif
+
+
+! For a subterranean 3D Variable
+!-------------------------------
+
+        else if (dim.eq.-3) then
+	! get variables' ID, if it exists
+	ierr=NF_INQ_VARID(nid,nom,varid)
+	
+	 if (ierr.ne.NF_NOERR) then ! variable not defined yet
+	  ! build related coordinates
+	  ierr=NF_INQ_DIMID(nid,"longitude",id(1))
+	  ierr=NF_INQ_DIMID(nid,"latitude",id(2))
+	  ierr=NF_INQ_DIMID(nid,"subsurface_layers",id(3))
+	  if (ierr.ne.NF_NOERR) then
+	   write(*,*)"write_archive: dimension <subsurface_layers>",
+     &               " is missing !!!"
+	   call abort
+	  endif
+          ierr=NF_INQ_DIMID(nid,"Time",id(4))
+	  
+	  ! define the variable
+	  write(*,*)"====================="
+	  write(*,*)"defining ",nom
+	  call def_var(nid,nom,titre,unite,4,id,varid,ierr)
+	  
+	 endif
+
+        ! build cedges and corners
+        corner(1)=1
+        corner(2)=1
+        corner(3)=1
+        corner(4)=ntime
+
+        edges(1)=iip1
+        edges(2)=jjp1
+        edges(3)=nsoilmx
+        edges(4)=1
+#ifdef NC_DOUBLE
+           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,px)
+#else
+           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,px)
+#endif
+
+! For a 3D ocean temperature Variable
+!-------------------------------
+
+        else if (dim.eq.-2) then
+        ! get variables' ID, if it exists
+        ierr=NF_INQ_VARID(nid,nom,varid)
+
+         if (ierr.ne.NF_NOERR) then ! variable not defined yet
+          ! build related coordinates
+          ierr=NF_INQ_DIMID(nid,"longitude",id(1))
+          ierr=NF_INQ_DIMID(nid,"latitude",id(2))
+          ierr=NF_INQ_DIMID(nid,"ocean_layers",id(3))
+          if (ierr.ne.NF_NOERR) then
+           write(*,*)"write_archive: dimension <ocean_layers>",
+     &               " is missing !!!"
+           call abort
+          endif
+          ierr=NF_INQ_DIMID(nid,"Time",id(4))
+
+          ! define the variable
+          write(*,*)"====================="
+          write(*,*)"defining ",nom
+          call def_var(nid,nom,titre,unite,4,id,varid,ierr)
+
+         endif
+
+        ! build cedges and corners
+        corner(1)=1
+        corner(2)=1
+        corner(3)=1
+        corner(4)=ntime
+
+        edges(1)=iip1
+        edges(2)=jjp1
+        edges(3)=noceanmx
+        edges(4)=1
+#ifdef NC_DOUBLE
+           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,px)
+#else
+           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,px)
+#endif
+
+
+
+! For a surface 2D Variable
+!--------------------------
+
+        else if (dim.eq.2) then
+
+!         Ecriture du champs
+
+           ierr= NF_INQ_VARID(nid,nom,varid)
+           if (ierr /= NF_NOERR) then
+!  choix du nom des coordonnees
+              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
+              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
+              ierr= NF_INQ_DIMID(nid,"Time",id(3))
+
+! Creation de la variable si elle n'existait pas
+
+              write (*,*) "====================="
+              write (*,*) "creation de ",nom
+
+              call def_var(nid,nom,titre,unite,3,id,varid,ierr)
+
+           endif
+
+           corner(1)=1
+           corner(2)=1
+           corner(3)=ntime
+           edges(1)=iip1
+           edges(2)=jjp1
+           edges(3)=1
+
+
+#ifdef NC_DOUBLE
+           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,px)
+#else         
+           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,px)
+#endif     
+
+           if (ierr.ne.NF_NOERR) then
+              write(*,*) "***** PUT_VAR matter in write_archive"
+              write(*,*) "***** with ",nom,nf_STRERROR(ierr)
+              call abort
+           endif
+
+
+!Cas Variable 0D (scalaire dependant du temps)
+!---------------------------------------------
+
+        else if (dim.eq.0) then
+
+!         Ecriture du champs
+
+           ierr= NF_INQ_VARID(nid,nom,varid)
+           if (ierr /= NF_NOERR) then
+!  choix du nom des coordonnees
+              ierr= NF_INQ_DIMID(nid,"Time",id(1))
+
+! Creation de la variable si elle n'existait pas
+
+              write (*,*) "====================="
+              write (*,*) "creation de ",nom
+
+              call def_var(nid,nom,titre,unite,1,id,varid,ierr)
+
+           endif
+
+           corner(1)=ntime
+           edges(1)=1
+
+#ifdef NC_DOUBLE
+           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,px)
+#else
+           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,px)
+#endif
+           if (ierr.ne.NF_NOERR) then
+              write(*,*) "***** PUT_VAR matter in write_archive"
+              write(*,*) "***** with ",nom,nf_STRERROR(ierr)
+              call abort
+           endif
+
+        else
+	  write(*,*) "write_archive: dim=",dim," ?!?"
+	  call abort
+        endif ! of if (dim.eq.3) else if (dim.eq.-3) ....
+
+      return
+      end
+
