Index: trunk/LMDZ.TITAN/libf/phytitan/profile.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/profile.F	(revision 893)
+++ trunk/LMDZ.TITAN/libf/phytitan/profile.F	(revision 894)
@@ -1,3 +1,3 @@
-      SUBROUTINE profile(unit,nlev,dzst,temp)
+      SUBROUTINE profile(unit,nlev,dzst,pres,temp)
       IMPLICIT NONE
 c=======================================================================
@@ -10,7 +10,8 @@
 c   differents profils d'atmospheres. T=f(z)
 c   entree:
-c     unit    unite de lecture de "testphys1d.def"
+c     unit    unite de lecture de "rcm1d.def"
 c     nlev    nombre de niveaux (nlev=llm+1, surf + couches 1 a llm)
 c     dzst    dz/T (avec dz = epaisseur de la couche en m)
+c     pres    pressure profile
 c     ichoice choix de l'atmosphere:
 c             1 Temperature constante
@@ -40,5 +41,5 @@
 
        INTEGER nlev, unit
-       REAL dzst(nlev),temp(nlev)
+       REAL dzst(nlev),pres(nlev),temp(nlev)
 
 c   local:
@@ -51,4 +52,5 @@
       REAL hauteur,tmp
       REAL zkm(nlev)    ! altitude en km
+      real a1,b1,c1,a2,b2,c2
 
       isin = 0
@@ -86,17 +88,18 @@
 
       ELSE IF(ichoice.EQ.2) THEN
-       print*,"Profil T a faire..."
-       stop
-c        temp(1) = 735.
-c        zkm(1)  = 0.0
-c        DO il=2,nlev
-c           zkm(il) = zkm(il-1)+temp(il-1)*dzst(il)/1000. ! approx avec T(l-1)
-c           if(zkm(il).lt.60.) then
-c             temp(il)=735.-7.95*zkm(il)
-c           else
-c             temp(il)=AMAX1(258.-3.*(zkm(il)-60.),168.)
-c           endif
-c           zkm(il) = zkm(il-1)+(temp(il-1)+temp(il))/2.*dzst(il)/1000.
-c        ENDDO
+       a1 =       142.1 
+       b1 =      -21.45 
+       c1 =       40.11 
+       a2 =       106.3 
+       b2 =       3183. 
+       c2 =       4737. 
+       DO il=1,nlev
+         temp(il)=a1*exp(-((pres(il)-b1)/c1)^2) 
+     .          + a2*exp(-((pres(il)-b2)/c2)^2)
+       ENDDO
+       zkm(1)  = 0.0
+       DO il=2,nlev
+          zkm(il) = zkm(il-1)+(temp(il-1)+temp(il))/2.*dzst(il)/1000.
+       ENDDO
 
 c-----------------------------------------------------------------------
Index: trunk/LMDZ.TITAN/libf/phytitan/rcm1d.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/rcm1d.F	(revision 893)
+++ trunk/LMDZ.TITAN/libf/phytitan/rcm1d.F	(revision 894)
@@ -78,5 +78,5 @@
       REAL pk(llm),pks, w(llm)
       INTEGER l, ierr, aslun
-      REAL tmp1(0:llm),tmp2(0:llm)                        
+      REAL tmp1(0:llm),tmp2(0:llm),tmp3(0:llm)                        
 
       character*2 str2
@@ -85,4 +85,5 @@
       COMMON/cpdetvenus/cppdyn,nu_venus,t0_venus
       REAL cppdyn,nu_venus,t0_venus
+      real pi
 
 c=======================================================================
@@ -279,17 +280,20 @@
         tmp1(ilayer)=-log(play(ilayer)/play(ilayer-1))*r/g
       ENDDO
-      call profile(unit,nlayer+1,tmp1,tmp2)
+      DO ilayer=0,nlayer
+        tmp2(ilayer)=plev(ilayer+1)
+      ENDDO
+      call profile(unit,nlayer+1,tmp1,tmp2,tmp3)
       CLOSE(unit)
 
       print*,"               Pression        Altitude     Temperature"
       ilayer=1
-      tsurf(1)=tmp2(0)
-       temp(1)=tmp2(1)
-       zlay(1)=tmp2(1)*tmp1(1)
+      tsurf(1)=tmp3(0)
+       temp(1)=tmp3(1)
+       zlay(1)=tmp3(1)*tmp1(1)
       print*,"           0",tsurf(1)
       print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer)
       DO ilayer=2,nlayer
-        temp(ilayer)=tmp2(ilayer)
-        zlay(ilayer)=zlay(ilayer-1)+tmp2(ilayer)*tmp1(ilayer)
+        temp(ilayer)=tmp3(ilayer)
+        zlay(ilayer)=zlay(ilayer-1)+tmp3(ilayer)*tmp1(ilayer)
         print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer)
       ENDDO
Index: trunk/LMDZ.VENUS/libf/phyvenus/profile.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/profile.F	(revision 893)
+++ trunk/LMDZ.VENUS/libf/phyvenus/profile.F	(revision 894)
@@ -10,5 +10,5 @@
 c   differents profils d'atmospheres. T=f(z)
 c   entree:
-c     unit    unite de lecture de "testphys1d.def"
+c     unit    unite de lecture de "rcm1d.def"
 c     nlev    nombre de niveaux (nlev=llm+1, surf + couches 1 a llm)
 c     dzst    dz/T (avec dz = epaisseur de la couche en m)
Index: trunk/LMDZ.VENUS/libf/phyvenus/rcm1d.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/rcm1d.F	(revision 893)
+++ trunk/LMDZ.VENUS/libf/phyvenus/rcm1d.F	(revision 894)
@@ -52,5 +52,5 @@
       REAL gru,grv   ! prescribed "geostrophic" background wind
       REAL temp(llm)   ! temperature at the middle of the layers
-      REAL q(llm,nqtot) ! tracer mixing ratio (e.g. kg/kg)
+      REAL,allocatable :: q(:,:) ! tracer mixing ratio (e.g. kg/kg)
       REAL tsoil(nsoilmx)   ! subsurface soik temperature (K)
       REAL zlay(llm)   ! altitude estimee dans les couches (km)
@@ -67,6 +67,5 @@
       REAL dudyn(llm),dvdyn(llm),dtempdyn(llm)
       REAL dpsurf    
-      REAL dq(llm,nqtot)
-      REAL dqdyn(llm,nqtot)
+      REAL,allocatable :: dq(:,:)
 
 c   Various intermediate variables
@@ -82,4 +81,5 @@
       COMMON/cpdetvenus/cppdyn,nu_venus,t0_venus
       REAL cppdyn,nu_venus,t0_venus
+      real pi
 
 c=======================================================================
@@ -113,4 +113,15 @@
       r= 8.314511E+0 *1000.E+0/mugaz
       rcp= r/cpp
+
+c-----------------------------------------------------------------------
+c   Initialisation des traceurs
+c   ---------------------------
+c  Choix du nombre de traceurs et du schema pour l'advection
+c  dans fichier traceur.def
+      call infotrac_init
+
+c Allocation de la tableau q : champs advectes   
+      allocate(q(llm,nqtot))
+      allocate(dq(llm,nqtot))
 
 c ------------------------------------------------------
Index: trunk/LMDZ.VENUS/libf/phyvenus/testphys1d.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/testphys1d.F	(revision 893)
+++ 	(revision )
@@ -1,481 +1,0 @@
-
-      PROGRAM testphys1d
-      
-      use control_mod
-      IMPLICIT NONE
-
-c=======================================================================
-c   subject:
-c   --------
-c   PROGRAM useful to run physical part of the venusian GCM in a 1D column
-c       
-c Can be compiled with a command like (e.g. for 50 layers)
-c  "makegcm -p venus -d 50 testphys1d"
-
-c A REVOIR POUR VENUS....
-c It requires the files "testphys1d.def" "callphys.def"
-c      and a file describing the sigma layers (e.g. "z2sig.def")
-c
-c   author: Frederic Hourdin, R.Fournier,F.Forget (original Mars version)
-c   ------- Sebastien Lebonnois (Venus version)
-c   
-c=======================================================================
-
-c VENUS: revoir tous les include...... par rapport a calfis.F
-c        voir aussi les arguments de calfis: clesphy0
-c        voir aussi dans les include de physiq.F les variables qui 
-c             doivent etre initialisees...
-#include "dimensions.h"
-#include "dimsoil.h"
-#include "comcstfi.h"
-!#include "control.h"
-#include "comvert.h"
-#include "netcdf.inc"
-#include "comg1d.h"
-#include "logic.h"
-#include "clesphys.h"
-#include "iniprint.h"
-
-c --------------------------------------------------------------
-c  Declarations
-c --------------------------------------------------------------
-c
-      INTEGER unit           ! unite de lecture de "testphys1d.def"
-      INTEGER unitstart      ! unite d'ecriture de "startphy.nc"
-      INTEGER nlayer,nlevel,nsoil,ndt
-      INTEGER ilayer,ilevel,isoil,idt,iq
-      LOGICAl firstcall,lastcall
-c
-      INTEGER day0          ! date initial (sol ; =0 a Ls=0)
-      REAL day              ! date durant le run
-      REAL time             ! time (0<time<1 ; time=0.5 a midi)
-      REAL play(llm)   ! Pressure at the middle of the layers (Pa)
-      REAL plev(llm+1) ! intermediate pressure levels (pa)
-      REAL psurf,tsurf(1)      
-      REAL u(llm),v(llm)  ! zonal, meridional wind
-      REAL gru,grv   ! prescribed "geostrophic" background wind
-      REAL temp(llm)   ! temperature at the middle of the layers
-      REAL q(llm,nqtot) ! tracer mixing ratio (e.g. kg/kg)
-      REAL tsoil(nsoilmx)   ! subsurface soik temperature (K)
-      REAL zlay(llm)   ! altitude estimee dans les couches (km)
-      REAL long(1),lati(1),area(1)
-      REAL cufi(1),cvfi(1)
-      REAL phisfi(1),albedo(1)
-      REAL solsw(1),sollwdown(1),dlw(1),radsol(1)
-
-c    Physical and dynamical tandencies (e.g.  m.s-2, K/s, Pa/s)
-      REAL du(llm),dv(llm),dtemp(llm)
-      REAL dudyn(llm),dvdyn(llm),dtempdyn(llm)
-      REAL dpsurf    
-      REAL dq(llm,nqtot)
-      REAL dqdyn(llm,nqtot)
-
-c   Various intermediate variables
-      REAL zls
-      REAL phi(llm),s(llm),aps(llm),bps(llm)
-      REAL pk(llm),pks, w(llm)
-      INTEGER l, ierr, aslun,radpas
-      REAL tmp1(0:llm),tmp2(0:llm)
-      INTEGER        longcles
-      PARAMETER    ( longcles = 20 )
-      REAL clesphy0( longcles      )
-                        
-
-      character*2 str2
-
-c=======================================================================
-
-c=======================================================================
-c INITIALISATION
-c=======================================================================
-
-      lunout = 6
-
-c ------------------------------------------------------
-c  Constantes prescrites ICI
-c ------------------------------------------------------
-
-      pi=2.E+0*asin(1.E+0)
-
-c     Constante de la Planete Venus
-c     -----------------------------
-      rad=6051300.               ! rayon de Venus (m)  ~6051300 m
-      daysec= 1.0087e7           ! duree du sol (s)  ~1.e7 s
-      omeg=4.*asin(1.)/19.4141e6 ! vitesse de rotation (rad.s-1)
-      g= 8.87                    ! gravite (m.s-2) ~8.87
-      mugaz=43.44                ! Masse molaire de l'atm (g.mol-1) ~43.44
-! ADAPTATION GCM POUR CP(T)
-! VENUS: Cp(T) = cpp*(T/T0)^nu 
-! avec T0=460. et nu=0.35
-!     cpp=1.0e3
-      cpp=9.0e2      ! version constante
-      r= 8.314511E+0 *1000.E+0/mugaz
-      rcp= r/cpp
-
-c ------------------------------------------------------
-c  Lecture des parametres dans "testphys1d.def" 
-c ------------------------------------------------------
-
-c   Opening parameters file "testphys1d.def"
-c   ---------------------------------------
-      unit =97
-      OPEN(unit,file='testphys1d.def',status='old',form='formatted'
-     .     ,iostat=ierr)
-
-      IF(ierr.ne.0) THEN
-        write(*,*) 'Problem to open "testphys1d.def'
-        write(*,*) 'Is it there ?'
-        stop
-      END IF
-
-c  Date et heure locale du debut du run
-c  ------------------------------------
-c    Date (en sols depuis le solstice de printemps) du debut du run
-      day0 = 0
-      PRINT *,'date de depart ?'
-      READ(unit,*) day0
-      day=float(day0)
-      PRINT *,day0
-c  Heure de demarrage
-      PRINT *,'heure de debut de simulation (entre 0 et 24) ?'
-      READ(unit,*) time
-      time=time/24.E+0
-
-c  Discretisation (Definition de la grille et des pas de temps)
-c  --------------
-c
-      nlayer=llm
-      nlevel=nlayer+1
-      nsoil=nsoilmx
-      PRINT *,'nombre de pas de temps par jour ?'
-      READ(unit,*) day_step
-      print*,day_step
-
-      PRINT *,'nombre d appel au rayonnement par jour ?'
-      READ(unit,*) nbapp_rad
-      print*,nbapp_rad
-
-      PRINT *,'nombre de pas physique entre chaque sortie ?'
-      READ(unit,*) nbphypers
-      print*,nbphypers
-
-      PRINT *,'nombre de jours simules ?'
-      READ(unit,*) ndt
-      print*,ndt
-
-      ndt=ndt*day_step     
-      dtphys=daysec/day_step  
-
-c Pression de surface sur la planete
-c ------------------------------------
-c
-      PRINT *,'pression au sol'
-      READ(unit,*) psurf
-      PRINT *,psurf
-c Pression de reference  ! voir dyn3d/etat0_venus
-c     pa     =  5.e4 
-      pa     =  1.e6 
-      preff  = 9.2e6 ! 92 bars
-c     preff  = psurf
- 
-c  latitude/longitude 
-c  -------------------
-      PRINT *,'latitude en degres ?'
-      READ(unit,*) lati(1)
-      PRINT *,lati(1)
-      long(1)=0.E+0
-
-c  Initialisation albedo 
-c  ----------------------
-c
-      PRINT *,'Albedo du sol nu ?'
-      READ(unit,*) albedo(1)
-      PRINT *,albedo(1)
-
-c   Initialisation speciales "physiq"
-c   ---------------------------------
-c   la surface de chaque maille est inutile en 1D --->
-      area(1)=1.E+0
-c de meme ?
-      cufi(1)=1.E+0
-      cvfi(1)=1.E+0
-
-c   le geopotentiel au sol est inutile en 1D car tout est controle
-c   par la pression de surface --->
-      phisfi(1)=0.E+0
-
-      CALL iniphysiq(1,llm,daysec,day0,dtphys,
-     .            lati,long,area,cufi,cvfi,rad,g,r,cpp)
-
-c   Initialisation pour prendre en compte les vents en 1-D
-c   ------------------------------------------------------
- 
-c    vent geostrophique
-      PRINT *,'composante vers l est du vent geostrophique (U) ?'
-      READ(unit,*) gru
-      PRINT *,'composante vers le nord du vent geostrophique (V) ?'
-      READ(unit,*) grv
-
-c     Initialisation des vents  au premier pas de temps
-      DO ilayer=1,nlayer
-         u(ilayer)=gru
-         v(ilayer)=grv
-      ENDDO
-
-c  calcul des pressions et altitudes en utilisant les niveaux sigma
-c  ----------------------------------------------------------------
-
-c    Vertical Coordinates  (hybrids)
-c    """"""""""""""""""""
-      CALL  disvert_noterre
-      
-c     Calcul au milieu des couches : Vient de la version Mars
-c     WARNING : le choix de placer le milieu des couches au niveau de
-c     pression intermédiaire est arbitraire et pourrait etre modifié.
-c     C'est fait de la meme facon dans disvert
-
-      DO l = 1, llm
-       aps(l) =  0.5 *( ap(l) +ap(l+1))
-       bps(l) =  0.5 *( bp(l) +bp(l+1))
-      ENDDO
-
-      DO ilevel=1,nlevel
-        plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
-      ENDDO
-
-      DO ilayer=1,nlayer
-        play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
-        pk(ilayer)  =cpp*(play(ilayer)/preff)**rcp
-c       write(120,*) ilayer,plev(ilayer),play(ilayer)
-      ENDDO
-c     write(120,*) nlevel,plev(nlevel)
-c     stop
-      
-      pks=cpp*(psurf/preff)**rcp
-
-c  profil de temperature et altitude au premier appel
-c  --------------------------------------------------
-
-c modif par rapport a Mars: 
-c   on envoie dz/T=-log(play/psurf)*r/g dans profile
-      tmp1(0)=0.0
-      tmp1(1)= -log(play(1)/psurf)*r/g
-      DO ilayer=2,nlayer
-        tmp1(ilayer)=-log(play(ilayer)/play(ilayer-1))*r/g
-      ENDDO
-      call profile(unit,nlayer+1,tmp1,tmp2)
-c file testphys1d.def closed in profile!
-
-      print*,"               Pression        Altitude     Temperature"
-      ilayer=1
-      tsurf(1)=tmp2(0)
-       temp(1)=tmp2(1)
-       zlay(1)=tmp2(1)*tmp1(1)
-      print*,"           0",tsurf(1)
-      print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer)
-      DO ilayer=2,nlayer
-        temp(ilayer)=tmp2(ilayer)
-        zlay(ilayer)=zlay(ilayer-1)+tmp2(ilayer)*tmp1(ilayer)
-        print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer)
-      ENDDO
-      
-c     temperature du sous-sol
-c     ~~~~~~~~~~~~~~~~~~~~~~~
-      DO isoil=1,nsoil
-         tsoil(isoil)=tsurf(1)
-      ENDDO
-
-c    Initialisation des traceurs
-c    ---------------------------
-
-      DO iq=1,nqtot
-        DO ilayer=1,nlayer
-           q(ilayer,iq) = 0.
-        ENDDO
-      ENDDO
-
-c    Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl"
-c    ----------------------------------------------------------------
-c    (effectuee avec "writeg1d" a partir de "physiq.F"
-c    ou tout autre subroutine physique, y compris celle ci).
-
-        g1d_nlayer=nlayer
-        g1d_nomfich='g1d.dat'
-        g1d_unitfich=40
-        g1d_nomctl='g1d.ctl'
-        g1d_unitctl=41
-        g1d_premier=.true.
-        g2d_premier=.true.
-
-c  Ecriture de "startphy.nc"
-c  -------------------------
-c  (Ce fichier sera aussitot relu au premier
-c   appel de "physiq", mais il est necessaire pour passer
-c   les variables purement physiques a "physiq"...
-
-      clesphy0(1) = 0.
-      clesphy0(2) = nbapp_rad
-      clesphy0(3) = 0. ! cycle diurne
-      clesphy0(4) = 1. ! soil_model
-      clesphy0(5) = 1. ! new_oliq
-      clesphy0(6) = 0.
-      clesphy0(7) = 0.
-      clesphy0(8) = 0.
-
-      solsw(1)    = 0.
-      sollwdown(1)= 0.
-      dlw(1)      = 0.
-      radsol(1)   = 0.
-     
-      radpas      = NINT(day_step/clesphy0(2))
-      soil_model  = .true.
-      new_oliq    = .true.
-
-      call phyredem("startphy.nc  ",dtphys,radpas,
-     .              lati,long,
-     .              tsurf,tsoil,albedo,
-     .              solsw,sollwdown,dlw,radsol,
-     .              temp,q)
-
-c=======================================================================
-c  BOUCLE TEMPORELLE DU MODELE 1D 
-c=======================================================================
-c
-      firstcall=.true.
-      lastcall=.false.
-
-      DO idt=1,ndt
-c        IF (idt.eq.ndt) lastcall=.true.
-        IF (idt.eq.ndt-day_step-1) then       !test
-         lastcall=.true.
-c toujours nulle dans le cas de Venus, pour l'instant...
-         zls = 0.0
-c        write(103,*) 'Ls=',zls*180./pi
-c        write(103,*) 'Lat=', lati(1)
-c        write(103,*) 'RunEnd - Atmos. Temp. File'
-c        write(103,*) 'RunEnd - Atmos. Temp. File'
-c        write(104,*) 'Ls=',zls*180./pi
-c        write(104,*) 'Lat=', lati(1)
-c        write(104,*) 'RunEnd - Atmos. Temp. File'
-        ENDIF
-
-c    calcul du geopotentiel 
-c     ~~~~~~~~~~~~~~~~~~~~~
-! ADAPTATION GCM POUR CP(T)
-      DO ilayer=1,nlayer
-        s(ilayer)=(play(ilayer)/psurf)**rcp
-      ENDDO
-      phi(1)=cpp*temp(1)*(1.E+0-s(1))
-      DO ilayer=2,nlayer
-         phi(ilayer)=phi(ilayer-1)+
-     &     cpp*(temp(ilayer-1)/s(ilayer-1)+temp(ilayer)/s(ilayer))*0.5
-     &        *(s(ilayer-1)-s(ilayer))
-
-      ENDDO
-
-c       appel de la physique
-c       --------------------
-
-      CALL physiq (1,llm,nqtot,
-     ,     firstcall,lastcall,
-     ,     day,time,dtphys,
-     ,     plev,play,pk,phi,phisfi,
-     ,     presnivs,clesphy0,
-     ,     u,v,temp,q,  
-     ,     w,
-C - sorties
-     s     du,dv,dtemp,dq,dpsurf)
-
-c     print*,"DT APRES PHYSIQ=",day,time
-c     print*,dtemp
-c     print*,temp
-c     print*," "
-c     stop
-
-c       evolution du vent : modele 1D
-c       -----------------------------
- 
-c       la physique calcule les derivees temporelles de u et v.
-c       Pas de coriolis 
-          DO ilayer=1,nlayer
-             du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
-             dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
-          ENDDO
-c      
-c         call writeg1d(1,nlayer,du,"du","du")
-c         call writeg1d(1,nlayer,dv,"dv","dv")
-c
-c       Calcul du temps au pas de temps suivant
-c       ---------------------------------------
-        firstcall=.false.
-        time=time+dtphys/daysec
-        IF (time.gt.1.E+0) then
-            time=time-1.E+0
-            day=day+1
-        ENDIF
-
-c       calcul des vitesses et temperature au pas de temps suivant
-c       ----------------------------------------------------------
-
-        DO ilayer=1,nlayer
-           u(ilayer)=u(ilayer)+dtphys*du(ilayer)
-           v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
-           temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
-        ENDDO
-c         call writeg1d(1,nlayer,u,"u","Vent zonal")
-c         call writeg1d(1,nlayer,v,"v","Vent meridien")
-
-c       calcul des pressions au pas de temps suivant
-c       ----------------------------------------------------------
-
-           psurf=psurf+dtphys*dpsurf   ! evolution de la pression de surface
-           DO ilevel=1,nlevel
-             plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
-           ENDDO
-           DO ilayer=1,nlayer
-             play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
-           ENDDO
-
-           call writeg1d(1,nlayer,play,"press","Pressure")
-           call writeg1d(1,nlayer,phi,"phi","Geopot")
-
-      ENDDO   ! fin de la boucle temporelle
-
-c    ========================================================
-c    GESTION DES SORTIE
-c    ========================================================
-
-c    fermeture pour conclure les sorties format grads dans "g1d.dat"
-c    et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F"
-c    ou tout autre subroutine physique, y compris celle ci).
-
-        CALL endg1d(1,nlayer,zlay/1000.,ndt)
-        print*,"Temperature finale:"
-        print*,temp
-        
-c stabilite
-      DO ilayer=1,nlayer
-        zlay(ilayer) = phi(ilayer)/8870.  !en km
-      ENDDO
-      DO ilayer=2,nlayer
-        tmp1(ilayer) =
-     .    (temp(ilayer)-temp(ilayer-1))/(zlay(ilayer)-zlay(ilayer-1)) 
-     .   + 8870./cpp
-      ENDDO
-
-      OPEN(11,file='profile.new')
-      write (11,*) tsurf
-      DO ilayer=1,nlayer
-        write (11,*) zlay(ilayer),temp(ilayer),tmp1(ilayer)
-      ENDDO
-
-c    ========================================================
-      END
- 
-c***********************************************************************
-c***********************************************************************
-
-#include "../dyn3d/disvert_noterre.F"
-#include "../dyn3d/abort_gcm.F"
-#include "../dyn3d/dump2d.F"
-
