Index: trunk/LMDZ.TITAN/libf/phytitan/effg.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/effg.F	(revision 855)
+++ trunk/LMDZ.TITAN/libf/phytitan/effg.F	(revision 887)
@@ -6,8 +6,8 @@
 !    (et dans physique, attention a clmain et autres...)
 
-!      EFFG = RG * (RA/(RA + Z ) )**2
+      EFFG = RG * (RA/(RA + Z ) )**2
 
 ! Pour l'instant:
-      EFFG = RG 
+!     EFFG = RG 
       RETURN
       END
Index: trunk/LMDZ.TITAN/libf/phytitan/physiq.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/physiq.F	(revision 855)
+++ trunk/LMDZ.TITAN/libf/phytitan/physiq.F	(revision 887)
@@ -780,7 +780,7 @@
       DO l=1,klev
          DO i=1,klon
-            zzlay(i,l)=zphi(i,l)/RG
+c           zzlay(i,l)=zphi(i,l)/RG
 c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
-c           zzlay(i,l)=RG*RA*RA/(RG*RA-zphi(i,l))-RA
+            zzlay(i,l)=RG*RA*RA/(RG*RA-zphi(i,l))-RA
          ENDDO
       ENDDO
@@ -1548,38 +1548,4 @@
 #endif
 
-c ---------- Sorties testphys1d -------------
-      
-      if (klon.eq.1) then
-        call writeg1d(klon,klev,t_seri,"Temp","Temperature")
-        call writeg1d(klon,1,ftsol,"tsurf","Surface Temp")
-      DO k = 1, klev
-      DO i = 1, klon
-c        tmpout(i,k) = real(heat(i,k))
-         tmpout(i,k) = heat(i,k)
-      ENDDO
-      ENDDO
-        call writeg1d(klon,klev,tmpout,
-     .                 "heat","Solar heating")
-      DO k = 1, klev
-      DO i = 1, klon
-c        tmpout(i,k) = real(dtrad(i,k))
-         tmpout(i,k) = dtrad(i,k)
-      ENDDO
-      ENDDO
-        call writeg1d(klon,klev,tmpout,
-     .                 "dtrad","IR cooling")
-        call writeg1d(klon,klev,lwnet,"lwnet","Net LW flux")
-        call writeg1d(klon,klev,swnet,"swnet","Net SW flux")
-        call writeg1d(klon,klev,fluxt,"flux_vdf","Turbulent flux")
-        call writeg1d(klon,klev,flux_ajs,"flux_ajs","Dry adjust. flux")
-        call writeg1d(klon,klev,flux_ec,"flux_ec","Ec flux")
-c       call writeg1d(klon,1,solsw,"surfsw","Net SW flux at surface")
-c       call writeg1d(klon,1,sollw,"surflw","Net LW flux at surface")
-        call writeg1d(klon,1,radsol,"surfnet","Net flux at surface")
-        call writeg1d(klon,klev,d_t_vdf,"dt_vdf","DT from clmain")
-        call writeg1d(klon,klev,d_t_ajs,"dt_ajs","DT from ajsec")
-        call writeg1d(klon,klev,d_t_ec,"dt_ec","DT from Ec")
-      endif
-      
 c====================================================================
 c Si c'est la fin, il faut conserver l'etat de redemarrage
Index: trunk/LMDZ.TITAN/libf/phytitan/profile.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/profile.F	(revision 887)
+++ trunk/LMDZ.TITAN/libf/phytitan/profile.F	(revision 887)
@@ -0,0 +1,205 @@
+      SUBROUTINE profile(unit,nlev,dzst,temp)
+      IMPLICIT NONE
+c=======================================================================
+c     Subroutine utilisee dans le modele 1-D  "rcm1d"
+c     pour l'initialisation du profil atmospherique
+c=======================================================================
+c
+c   VERSION VENUS
+c
+c   differents profils d'atmospheres. T=f(z)
+c   entree:
+c     unit    unite de lecture de "testphys1d.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     ichoice choix de l'atmosphere:
+c             1 Temperature constante
+c             2 profil Huygens lisse
+c             3 
+c             4 
+c             5 
+c             6 T constante + perturbation gauss (level) (christophe 10/98)
+c             7 T constante + perturbation gauss   (km)  (christophe 10/98)
+c             8 Lecture du profile dans un fichier ASCII (profile)
+c     tref    temperature de reference
+c     isin    ajout d'une perturbation (isin=1)
+c     pic     pic perturbation gauss pour ichoice = 6 ou 7
+c     largeur largeur de la perturbation gauss pour ichoice = 6 ou 7
+c     hauteur hauteur de la perturbation gauss pour ichoice = 6 ou 7
+c
+c   sortie:
+c     temp    temperatures en K
+c     
+c=======================================================================
+c-----------------------------------------------------------------------
+c   declarations:
+c   -------------
+
+c   arguments:
+c   ----------
+
+       INTEGER nlev, unit
+       REAL dzst(nlev),temp(nlev)
+
+c   local:
+c   ------
+
+      INTEGER il,ichoice,isin,iter
+      REAL pi
+      REAL tref,t1,t2,t3,ww
+      REAL pic,largeur
+      REAL hauteur,tmp
+      REAL zkm(nlev)    ! altitude en km
+
+      isin = 0
+
+c-----------------------------------------------------------------------
+c   choix du profil:
+c   ----------------
+
+c la lecture se fait dans le rcm1d.def, ouvert par rcm1d.F
+      READ(unit,*) 
+      READ(unit,*) 
+      READ(unit,*) 
+      READ(unit,*) ichoice
+      READ(unit,*) tref
+      READ(unit,*) isin
+      READ(unit,*) pic
+      READ(unit,*) largeur
+      READ(unit,*) hauteur
+
+c-----------------------------------------------------------------------
+c   ichoice=1 temperature constante:
+c   --------------------------------
+
+      IF(ichoice.EQ.1) THEN
+         temp(1) = tref
+         zkm(1)  = 0.0
+         DO il=2,nlev
+            temp(il)= tref
+            zkm(il) = zkm(il-1)+temp(il)*dzst(il)/1000.
+         ENDDO
+
+c-----------------------------------------------------------------------
+c   ichoice=2 Huygens lisse:
+c   ------------------------
+
+      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
+
+c-----------------------------------------------------------------------
+c   ichoice=3
+c   ----------------------------
+
+      ELSE IF(ichoice.EQ.3) THEN
+       print*,"Profil T a faire..."
+       stop
+
+c-----------------------------------------------------------------------
+c   ichoice=4 :
+c   ------------------
+
+      ELSE IF(ichoice.EQ.4) THEN
+         print*,"Cas non defini..."
+         print*,"Stop dans profile.F"
+         STOP
+
+c-----------------------------------------------------------------------
+c   ichoice=5 :
+c   ----------------
+
+      ELSE IF(ichoice.EQ.5) THEN
+         print*,"Cas non defini..."
+         print*,"Stop dans profile.F"
+         STOP
+
+c-----------------------------------------------------------------------
+c   ichoice=6 
+c   ---------
+
+      ELSE IF(ichoice.EQ.6) THEN
+      temp(1) = tref
+      zkm(1)  = 0.0
+      DO il=2,nlev
+        tmp=il-pic
+        temp(il)= tref + hauteur*exp(-tmp*tmp/largeur/largeur)
+        zkm(il) = zkm(il-1)+temp(il)*dzst(il)/1000.
+      ENDDO
+
+
+c-----------------------------------------------------------------------
+c   ichoice=7
+c   ---------
+
+      ELSE IF(ichoice.EQ.7) THEN
+      temp(1) = tref
+      zkm(1)  = 0.0
+      DO il=2,nlev
+        zkm(il) = zkm(il-1)+tref*dzst(il)/1000. ! approx
+        tmp=zkm(il)-pic
+        temp(il)= tref + hauteur*exp(-tmp*tmp*4/largeur/largeur)
+        zkm(il) = zkm(il-1)+(temp(il-1)+temp(il))/2.*dzst(il)/1000.
+      ENDDO
+
+c-----------------------------------------------------------------------
+c   ichoice=8
+c   ---------
+
+      ELSE IF(ichoice.GE.8) THEN
+      OPEN(11,file='profile',status='old',form='formatted',err=101)
+      DO il=1,nlev
+        READ (11,*) temp(il)
+      ENDDO
+      zkm(1) = 0.0
+      DO il=2,nlev
+        zkm(il) = zkm(il-1)+(temp(il-1)+temp(il))/2.*dzst(il)/1000.
+      ENDDO
+
+      GOTO 201
+101   STOP'fichier profile inexistant'
+201   CONTINUE
+      CLOSE(10)
+
+c-----------------------------------------------------------------------
+
+      ENDIF
+
+c-----------------------------------------------------------------------
+c   rajout eventuel d'une perturbation:
+c   -----------------------------------
+
+      IF(isin.EQ.1) THEN
+	 pi=2.*ASIN(1.)
+	 DO il=1,nlev
+        temp(il)=temp(il)+(1.-1000./(1000+zkm(il)*zkm(il)))*(
+     s      6.*SIN(zkm(il)*pi/6.)+9.*SIN(zkm(il)*pi/10.3) )
+	 ENDDO
+      ENDIF
+
+
+c-----------------------------------------------------------------------
+c   Ecriture du profil de temperature dans un fichier profile.out
+c   -------------------------------------------------------------
+
+
+c     OPEN(12,file='profile.out')
+c         DO il=1,nlev
+c            write(12,*) temp(il)
+c           write(12,*) temp(il),zkm(il)
+c         ENDDO
+c     CLOSE(12)
+
+      RETURN
+      END
Index: trunk/LMDZ.TITAN/libf/phytitan/rcm1d.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/rcm1d.F	(revision 887)
+++ trunk/LMDZ.TITAN/libf/phytitan/rcm1d.F	(revision 887)
@@ -0,0 +1,492 @@
+
+      PROGRAM rcm1d
+      
+      USE infotrac
+      use control_mod
+      use comgeomphy
+      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 55 layers)
+c  "makelmdz -p titan -d 55 rcm1d"
+
+c It requires the files "rcm1d.def" "physiq.def" "traceur.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 Version TITAN a tester et verifier  
+c  - verifier pour Ls... 
+c  - faire un profile.F ...
+
+#include "dimensions.h"
+#include "dimsoil.h"
+#include "comcstfi.h"
+#include "comvert.h"
+#include "netcdf.inc"
+#include "logic.h"
+#include "clesphys.h"
+#include "iniprint.h"
+#include "tabcontrol.h"
+
+c --------------------------------------------------------------
+c  Declarations
+c --------------------------------------------------------------
+c
+      INTEGER unit           ! unite de lecture de "rcm1d.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,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)
+      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)
+      REAL zmea(1), zstd(1)
+      REAL zsig(1), zgam(1), zthe(1)
+      REAL zpic(1), zval(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,allocatable :: dq(:,:)
+
+c   Various intermediate variables
+      REAL zls
+      REAL phi(llm),s(llm)
+      REAL pk(llm),pks, w(llm)
+      INTEGER l, ierr, aslun
+      REAL tmp1(0:llm),tmp2(0:llm)                        
+
+      character*2 str2
+
+c normalement dans dyn3d/comconst.h
+      COMMON/cpdetvenus/cppdyn,nu_venus,t0_venus
+      REAL cppdyn,nu_venus,t0_venus
+
+c=======================================================================
+
+c=======================================================================
+c INITIALISATION
+c=======================================================================
+
+      lunout = 6
+
+c ------------------------------------------------------
+c  Constantes prescrites ICI
+c ------------------------------------------------------
+
+      pi=2.E+0*asin(1.E+0)
+
+c     Constante de Titan
+c     -------------------
+      planet_type = "titan"
+      rad=2575000.               ! rayon de Venus (m) 
+      daysec=1.37889e6           ! duree du sol (s)  
+      omeg=2*pi/daysec           ! vitesse de rotation (rad.s-1)
+      g= 1.35                    ! gravite (m.s-2)
+      mugaz=28.                  ! Masse molaire de l'atm (g.mol-1) 
+      cpp=1.039e3
+      cppdyn=cpp
+      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 ------------------------------------------------------
+c  Lecture des parametres dans "rcm1d.def" 
+c ------------------------------------------------------
+
+c   Opening parameters file "rcm1d.def"
+c   ---------------------------------------
+      unit =97
+      OPEN(unit,file='rcm1d.def',status='old',form='formatted'
+     .     ,iostat=ierr)
+
+      IF(ierr.ne.0) THEN
+        write(*,*) 'Problem to open "rcm1d.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=REAL(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
+
+c     PRINT *,'nombre d appel au rayonnement par jour ?'
+c     READ(unit,*) nbapp_rad
+c     print*,nbapp_rad
+c LU DANS PHYSIQ.DEF...
+      nbapp_rad = 10.
+
+      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  
+      pa     =  1.e4 
+      preff  = 1.4e5 
+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 ne sert pas ici...
+      albedo(1)=0.3
+c      PRINT *,'Albedo du sol nu ?'
+c      READ(unit,*) albedo(1)
+c      PRINT *,albedo(1)
+
+c   Initialisation speciales "physiq"
+c   ---------------------------------
+
+      CALL init_phys_lmdz(iim,jjm,llm,1,(/1/))
+      call initcomgeomphy
+      call infotrac_init
+      call ini_cpdet
+
+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)
+      CLOSE(unit)
+
+      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)=93.
+      ENDDO
+
+c    Initialisation des traceurs
+c    ---------------------------
+
+      DO iq=1,nqtot
+        DO ilayer=1,nlayer
+           q(ilayer,iq) = 0.
+        ENDDO
+      ENDDO
+
+c    Initialisation des parametres d'oro
+c    -----------------------------------
+
+      zmea(1) = 0.
+      zstd(1) = 0.
+      zsig(1) = 0.
+      zgam(1) = 0.
+      zthe(1) = 0.
+      zpic(1) = 0.
+      zval(1) = 0.
+
+c  Initialisation Ls
+c ------------------
+         zls=0.
+         print*,'Ls=',zls*180./pi
+
+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"...
+
+      solsw(1)    = 0.
+      sollwdown(1)= 0.
+      dlw(1)      = 0.
+      radsol(1)   = 0.
+     
+      radpas      = NINT(1.*day_step/nbapp_rad)
+      soil_model  = .true.
+
+      call phyredem("startphy.nc  ",
+     .              lati,long,
+     .              tsurf,tsoil,albedo,
+     .              solsw,sollwdown,dlw,radsol,
+     .    zmea, zstd, zsig, zgam, zthe, zpic, zval, 
+     .              temp)
+
+c=======================================================================
+c  BOUCLE TEMPORELLE DU MODELE 1D 
+c=======================================================================
+c
+      firstcall=.true.
+      lastcall=.false.
+
+      DO idt=1,ndt
+        IF (idt.eq.ndt) then 
+         lastcall=.true.
+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,
+     ,     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       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       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
+
+      ENDDO   ! fin de la boucle temporelle
+
+c    ========================================================
+c    GESTION DES SORTIE
+c    ========================================================
+
+        print*,"Temperature finale:"
+        print*,temp
+        
+c stabilite
+      DO ilayer=1,nlayer
+        zlay(ilayer) = phi(ilayer)/g/1000.  !en km
+      ENDDO
+      DO ilayer=2,nlayer
+        tmp1(ilayer) =
+     .    (temp(ilayer)-temp(ilayer-1))/(zlay(ilayer)-zlay(ilayer-1)) 
+     .   + 1000.*g/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"
+#include "../dyn3d/cpdet.F"
+
+c***********************************************************************
+      function ssum(n,sx,incx)
+c
+      IMPLICIT NONE
+c
+      integer n,incx,i,ix
+      real ssum,sx((n-1)*incx+1)
+c
+      ssum=0.
+      ix=1
+      do 10 i=1,n
+         ssum=ssum+sx(ix)
+         ix=ix+incx
+10    continue
+c
+      return
+      end
