Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/caldyn0.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/caldyn0.F	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/caldyn0.F	(revision 2366)
@@ -1,1 +1,88 @@
-link ../../dyn3d/caldyn0.F
+      SUBROUTINE caldyn0
+     $ (itau,ucov,vcov,teta,ps,masse,pk,phis ,
+     $  phi,w,pbaru,pbarv,time )
+
+      USE comvert_mod, ONLY: ap,bp
+
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c  Auteur :  P. Le Van
+c
+c   Objet:
+c   ------
+c
+c   Calcul des tendances dynamiques.
+c
+c Modif 04/93 F.Forget
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c   0. Declarations:
+c   ----------------
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comgeom.h"
+
+c   Arguments:
+c   ----------
+
+      INTEGER itau
+      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
+      REAL ps(ip1jmp1),phis(ip1jmp1)
+      REAL pk(iip1,jjp1,llm)
+      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
+      REAL phi(ip1jmp1,llm),masse(ip1jmp1,llm)
+      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
+      REAL time
+
+c   Local:
+c   ------
+
+      REAL ang(ip1jmp1,llm),p(ip1jmp1,llmp1)
+      REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm),psexbarxy(ip1jm)
+      REAL vorpot(ip1jm,llm)
+      REAL w(ip1jmp1,llm),ecin(ip1jmp1,llm),convm(ip1jmp1,llm)
+      REAL bern(ip1jmp1,llm)
+      REAL massebxy(ip1jm,llm), dp(ip1jmp1)
+    
+
+      INTEGER   ij,l
+      EXTERNAL  advect,bernoui,convmas,covcont ,
+     *          enercin,flumass,tourpot,vitvert,sortvarc,traceur,
+     *          pression,psextbar,massdair
+
+c-----------------------------------------------------------------------
+c   Calcul des tendances dynamiques:
+c   --------------------------------
+      CALL covcont  (llm,ucov,vcov,ucont,vcont)
+      CALL pression (ip1jmp1,ap,bp,ps,p)
+      CALL psextbar (ps,psexbarxy)
+      CALL massdair (p,masse)
+      CALL massbar  (masse,massebx,masseby)
+      CALL massbarxy(masse,massebxy)
+      CALL flumass  (massebx,masseby,vcont,ucont,pbaru,pbarv)
+      CALL convmas  (pbaru,pbarv,convm)
+
+      DO ij =1, ip1jmp1
+         dp( ij ) = convm( ij,1 ) / airesurg( ij )
+      ENDDO
+
+      CALL vitvert (convm,w)
+      CALL tourpot (vcov,ucov,massebxy,vorpot)
+      CALL enercin (vcov,ucov,vcont,ucont,ecin)
+      CALL bernoui (ip1jmp1,llm,phi,ecin,bern)
+
+      DO l=1,llm
+         DO ij=1,ip1jmp1
+            ang(ij,l) = ucov(ij,l) + constang(ij)
+         ENDDO
+      ENDDO
+
+        CALL sortvarc0
+     $ (itau,ucov,teta,ps,masse,pk,phis,vorpot,phi,bern,dp,time,vcov)
+
+      RETURN
+      END
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/defrun_new.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/defrun_new.F	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/defrun_new.F	(revision 2366)
@@ -1,1 +1,535 @@
-link ../../dyn3d/defrun_new.F
+      SUBROUTINE defrun_new( tapedef, etatinit )
+#ifndef CPP_PARA
+c
+c-----------------------------------------------------------------------
+c     Auteurs :   L. Fairhead , P. Le Van  .
+c      Modif C. Hourdin F. Forget VERSION MARTIENNE
+c
+c
+c  -------------------------------------------------------------------
+c
+c                    MODIF JUIN 2000 (zoom)
+c       .........     Version  du 29/04/97       ..........
+c
+c   Nouveaux parametres nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,
+c   tetatemp   ajoutes  pour la dissipation   .
+c
+c   Autre parametre ajoute en fin de liste : ** fxyhypb ** 
+c
+c   Si fxyhypb = .TRUE. , choix de la fonction a derivee tangente hyperb.
+c   Sinon , choix de fxynew  , a derivee sinusoidale  ..
+c
+c   ......  etatinit = . TRUE. si defrun_new  est appele dans NEWSTART
+c   ETAT0_LMD  ou  LIMIT_LMD  pour l'initialisation de start.dat (dic) et
+c   de limit.dat (dic)  ...........
+c   Sinon  etatinit = . FALSE .
+c
+c   Donc etatinit = .F.  si on veut comparer les valeurs de  alphax ,
+c   alphay,clon,clat, fxyhypb  lues sur  le fichier  start  avec
+c   celles passees  par run.def ,  au debut du gcm, apres l'appel a 
+c   lectba .  
+c   Ces parametres definissant entre autres la grille et doivent etre
+c   pareils et coherents , sinon il y aura  divergence du gcm .
+c
+c
+c-----------------------------------------------------------------------
+c   Declarations :
+c   --------------
+! to use  'getin'
+      USE ioipsl_getincom
+      use sponge_mod,only: callsponge,nsponge,mode_sponge,tetasponge
+      use control_mod,only: nday, day_step, iperiod, anneeref,
+     &                      iconser, dissip_period, iphysiq
+      USE logic_mod, ONLY: hybrid,purmats,fxyhypb,ysinus,iflag_phys
+      USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy,
+     .			alphax,alphay,taux,tauy
+      IMPLICIT NONE
+
+#include "dimensions.h"
+#include "paramet.h"
+!#include "control.h"
+#include "comdissnew.h"
+#include "iniprint.h"
+c
+c   arguments:
+c   ---------
+      LOGICAL  etatinit ! should be .false. for a call from gcm.F
+                        ! and .true. for a call from newstart.F
+      INTEGER  tapedef  ! unit number to assign to 'run.def' file
+c
+c   local variables:
+c   ---------------
+
+      CHARACTER ch1*72,ch2*72,ch3*72,ch4*8 ! to store various strings
+      integer tapeerr ! unit number for error message
+      parameter (tapeerr=0)
+
+c     REAL clonn,clatt,alphaxx,alphayy
+c     LOGICAL  fxyhypbb
+      INTEGER ierr
+      REAL clonn,clatt,grossismxx,grossismyy
+      REAL dzoomxx,dzoomyy,tauxx,tauyy,temp
+      LOGICAL  fxyhypbb, ysinuss
+
+
+c   initialisations:
+c   ----------------
+ 
+!      lunout=6
+
+c-----------------------------------------------------------------------
+c  Parametres de controle du run:
+c-----------------------------------------------------------------------
+
+
+!Initialisation des parametres "terrestres", qui ne concernent pas
+!le modele martien et ne sont donc plus lues dans "run.def"
+
+        anneeref=0
+        ! Note: anneref is a common in 'control.h'
+
+      OPEN(tapedef,file='run.def',status='old',form='formatted'
+     .     ,iostat=ierr)
+      CLOSE(tapedef) ! first call to getin will open the file
+
+      !lunout: default unit for the text outputs
+      lunout=6
+      CALL getin('lunout', lunout)
+      IF (lunout /= 5 .and. lunout /= 6) THEN
+        OPEN(UNIT=lunout,FILE='lmdz.out',ACTION='write',
+     &       STATUS='unknown',FORM='formatted')
+      ENDIF
+
+      IF(ierr.EQ.0) THEN ! if file run.def is found
+        WRITE(lunout,*) "DEFRUN_NEW: reading file run.def"
+        
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "Number of days to run:"
+        nday=1 ! default value
+        call getin("nday",nday)
+        WRITE(lunout,*)" nday = ",nday
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "Number of dynamical steps per day:",
+     & "(should be a multiple of iperiod)"
+        day_step=960 ! default value
+        call getin("day_step",day_step)
+        WRITE(lunout,*)" day_step = ",day_step
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "periode pour le pas Matsuno (en pas)"
+        iperiod=5 ! default value
+        call getin("iperiod",iperiod)
+        WRITE(lunout,*)" iperiod = ",iperiod
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "periode de sortie des variables de ",
+     &  "controle (en pas)"
+        iconser=120 ! default value
+        call getin("iconser",iconser)
+        WRITE(lunout,*)" iconser = ",iconser
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "periode de la dissipation (en pas)"
+        dissip_period=5 ! default value
+        call getin("idissip",dissip_period)
+        ! if there is a "dissip_period" in run.def, it overrides "idissip"
+        call getin("dissip_period",dissip_period)
+        WRITE(lunout,*)" dissip_period = ",dissip_period
+
+ccc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
+ccc
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "choix de l'operateur de dissipation ",
+     & "(star ou  non star )"
+        lstardis=.true. ! default value
+        call getin("lstardis",lstardis)
+        WRITE(lunout,*)" lstardis = ",lstardis
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "avec ou sans coordonnee hybrides"
+        hybrid=.true. ! default value
+        call getin("hybrid",hybrid)
+        WRITE(lunout,*)" hybrid = ",hybrid
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "nombre d'iterations de l'operateur de ",
+     & "dissipation   gradiv "
+        nitergdiv=1 ! default value
+        call getin("nitergdiv",nitergdiv)
+        WRITE(lunout,*)" nitergdiv = ",nitergdiv
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "nombre d'iterations de l'operateur de ",
+     & "dissipation  nxgradrot"
+        nitergrot=2 ! default value
+        call getin("nitergrot",nitergrot)
+        WRITE(lunout,*)" nitergrot = ",nitergrot
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "nombre d'iterations de l'operateur de ",
+     & "dissipation  divgrad"
+        niterh=2 ! default value
+        call getin("niterh",niterh)
+        WRITE(lunout,*)" niterh = ",niterh
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "temps de dissipation des plus petites ",
+     & "long.d ondes pour u,v (gradiv)"
+        tetagdiv=4000. ! default value
+        call getin("tetagdiv",tetagdiv)
+        WRITE(lunout,*)" tetagdiv = ",tetagdiv
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "temps de dissipation des plus petites ",
+     & "long.d ondes pour u,v(nxgradrot)"
+        tetagrot=5000. ! default value
+        call getin("tetagrot",tetagrot)
+        WRITE(lunout,*)" tetagrot = ",tetagrot
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "temps de dissipation des plus petites ",
+     & "long.d ondes pour  h ( divgrad)"
+        tetatemp=5000. ! default value
+        call getin("tetatemp",tetatemp)
+        WRITE(lunout,*)" tetatemp = ",tetatemp
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "coefficient pour gamdissip"
+        coefdis=0. ! default value
+        call getin("coefdis",coefdis)
+        WRITE(lunout,*)" coefdis = ",coefdis
+c
+c    ...............................................................
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "choix du shema d'integration temporelle ",
+     & "(true==Matsuno ou false==Matsuno-leapfrog)"
+        purmats=.false. ! default value
+        call getin("purmats",purmats)
+        WRITE(lunout,*)" purmats = ",purmats
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "avec ou sans physique"
+!        physic=.true. ! default value
+!        call getin("physic",physic)
+!        WRITE(lunout,*)" physic = ",physic
+        iflag_phys=1 ! default value
+        call getin("iflag_phys",iflag_phys)
+        WRITE(lunout,*)" iflag_phys = ",iflag_phys
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "periode de la physique (en pas)"
+        iphysiq=20 ! default value
+        call getin("iphysiq",iphysiq)
+        WRITE(lunout,*)" iphysiq = ",iphysiq
+
+!        WRITE(lunout,*) ""
+!        WRITE(lunout,*) "choix d'une grille reguliere"
+!        grireg=.true.
+!        call getin("grireg",grireg)
+!        WRITE(lunout,*)" grireg = ",grireg
+
+ccc   .... P.Le Van, ajout le 03/01/96 pour l'ecriture phys ...
+c
+!        WRITE(lunout,*) ""
+!        WRITE(lunout,*) "frequence (en pas) de l'ecriture ",
+!     & "du fichier diagfi.nc"
+!        ecritphy=240
+!        call getin("ecritphy",ecritphy)
+!        WRITE(lunout,*)" ecritphy = ",ecritphy
+
+ccc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
+c     .........   (  modif  le 17/04/96 )   .........
+c
+        if (.not.etatinit ) then 
+
+           clonn=63.
+           call getin("clon",clonn)
+           
+           IF( ABS(clon - clonn).GE. 0.001 )  THEN
+             PRINT *,' La valeur de clon passee par run.def est '
+     *       ,'differente de celle lue sur le fichier start '
+             STOP
+           ENDIF
+c
+c
+           clatt=0.
+           call getin("clat",clatt)
+  
+           IF( ABS(clat - clatt).GE. 0.001 )  THEN
+             PRINT *,' La valeur de clat passee par run.def est '
+     *       ,'differente de celle lue sur le fichier start '
+             STOP
+           ENDIF
+
+           grossismxx=1.0
+           call getin("grossismx",grossismxx)
+
+           if(grossismxx.eq.0) then  
+             write(*,*)
+             write(*,*)'ERREUR : dans run.def, grossismx =0'
+             write(*,*)'Attention a ne pas utiliser une version de'
+             write(*,*)'run.def avant le nouveau zoom LMDZ2.3 (06/2000)'
+             write(*,*)'(Il faut ajouter grossismx,dzoomx,etc... a la'
+             write(*,*)'place de alphax, alphay. cf. dyn3d). '
+             write(*,*)
+             stop
+           end if
+
+           IF( ABS(grossismx - grossismxx).GE. 0.001 )  THEN
+             PRINT *,' La valeur de grossismx passee par run.def est '
+     *       ,'differente de celle lue sur le fichier  start =',
+     *        grossismx
+             if (grossismx.eq.0) then
+                  write(*,*) 'OK,Normal : c est un vieux start'
+     *             , 'd avant le nouveau zoom LMDZ2.3 (06/2000)'
+                 grossismx=grossismxx
+             else
+                   STOP
+             endif
+           ENDIF
+
+           grossismyy=1.0
+           call getin("grossismy",grossismyy)
+
+           IF( ABS(grossismy - grossismyy).GE. 0.001 )  THEN
+             PRINT *,' La valeur de grossismy passee par run.def est '
+     *       ,'differente de celle lue sur le fichier  start =',
+     *        grossismy
+             if (grossismy.eq.0) then
+                  write(*,*) 'OK,Normal : c est un vieux start'
+     *             , 'd avant le nouveau zoom LMDZ2.3 (06/2000)'
+                 grossismy=grossismyy
+             else
+                   STOP
+             endif
+           ENDIF
+
+
+           IF( grossismx.LT.1. )  THEN
+             PRINT *,' ***  ATTENTION !! grossismx < 1 .   *** '
+             STOP
+           ELSE
+             alphax = 1. - 1./ grossismx
+           ENDIF
+
+           IF( grossismy.LT.1. )  THEN
+             PRINT *,' ***  ATTENTION !! grossismy < 1 .   *** '
+             STOP
+           ELSE
+             alphay = 1. - 1./ grossismy
+           ENDIF
+
+           PRINT *,' '
+           PRINT *,' --> In defrun: alphax alphay  ',alphax,alphay
+           PRINT *,' '
+c
+           fxyhypbb=.false.
+           call getin("fxyhypbb",fxyhypbb)
+  
+           IF( .NOT.fxyhypb )  THEN
+             IF( fxyhypbb )     THEN
+                PRINT *,' ********  PBS DANS  DEFRUN  ******** '
+                PRINT *,' *** fxyhypb lu sur le fichier start est F ',
+     *          'alors  qu il est  T  sur  run.def  ***'
+                STOP
+             ENDIF
+           ELSE
+             IF( .NOT.fxyhypbb )   THEN
+                PRINT *,' ********  PBS DANS  DEFRUN  ******** '
+                PRINT *,' ***  fxyhypb lu sur le fichier start est T ',
+     *         'alors  qu il est  F  sur  run.def  ****  '
+                STOP
+             ENDIF
+           ENDIF
+           dzoomxx=0.0
+           call getin("dzoomx",dzoomxx)
+
+           IF( fxyhypb )  THEN
+              IF( ABS(dzoomx - dzoomxx).GE. 0.001 )  THEN
+                PRINT *,' La valeur de dzoomx passee par run.def est '
+     *          ,'differente de celle lue sur le fichier  start '
+                STOP
+              ENDIF
+           ENDIF
+
+           dzoomyy=0.0
+           call getin("dzoomy",dzoomyy)
+
+           IF( fxyhypb )  THEN
+              IF( ABS(dzoomy - dzoomyy).GE. 0.001 )  THEN
+                PRINT *,' La valeur de dzoomy passee par run.def est '
+     *          ,'differente de celle lue sur le fichier  start '
+                STOP
+              ENDIF
+           ENDIF
+
+           tauxx=2.0
+           call getin("taux",tauxx)
+
+           tauyy=2.0
+           call getin("tauy",tauyy)
+
+           IF( fxyhypb )  THEN
+              IF( ABS(taux - tauxx).GE. 0.001 )  THEN
+                WRITE(6,*)' La valeur de taux passee par run.def est', 
+     *             'differente de celle lue sur le fichier  start '
+                CALL ABORT
+              ENDIF
+
+              IF( ABS(tauy - tauyy).GE. 0.001 )  THEN
+                WRITE(6,*)' La valeur de tauy passee par run.def est',
+     *          'differente de celle lue sur le fichier  start '
+                CALL ABORT
+              ENDIF
+           ENDIF
+  
+        ELSE    ! Below, case when etainit=.true.
+
+           WRITE(lunout,*) ""
+           WRITE(lunout,*) "longitude en degres du centre du zoom"
+           clon=63. ! default value
+           call getin("clon",clon)
+           WRITE(lunout,*)" clon = ",clon
+           
+c
+           WRITE(lunout,*) ""
+           WRITE(lunout,*) "latitude en degres du centre du zoom "
+           clat=0. ! default value
+           call getin("clat",clat)
+           WRITE(lunout,*)" clat = ",clat
+
+           WRITE(lunout,*) ""
+           WRITE(lunout,*) "facteur de grossissement du zoom,",
+     & " selon longitude"
+           grossismx=1.0 ! default value
+           call getin("grossismx",grossismx)
+           WRITE(lunout,*)" grossismx = ",grossismx
+
+           WRITE(lunout,*) ""
+           WRITE(lunout,*) "facteur de grossissement du zoom ,",
+     & " selon latitude"
+           grossismy=1.0 ! default value
+           call getin("grossismy",grossismy)
+           WRITE(lunout,*)" grossismy = ",grossismy
+
+           IF( grossismx.LT.1. )  THEN
+            PRINT *,' ***  ATTENTION !! grossismx < 1 .   *** '
+            STOP
+           ELSE
+             alphax = 1. - 1./ grossismx
+           ENDIF
+
+           IF( grossismy.LT.1. )  THEN
+             PRINT *,' ***  ATTENTION !! grossismy < 1 .   *** '
+             STOP
+           ELSE
+             alphay = 1. - 1./ grossismy
+           ENDIF
+
+           PRINT *,' Defrun  alphax alphay  ',alphax,alphay
+c
+           WRITE(lunout,*) ""
+           WRITE(lunout,*) "Fonction  f(y)  hyperbolique  si = .true.",
+     &  ", sinon  sinusoidale"
+           fxyhypb=.false. ! default value
+           call getin("fxyhypb",fxyhypb)
+           WRITE(lunout,*)" fxyhypb = ",fxyhypb
+
+           WRITE(lunout,*) ""
+           WRITE(lunout,*) "extension en longitude de la zone du zoom", 
+     & " (fraction de la zone totale)"
+           dzoomx=0. ! default value
+           call getin("dzoomx",dzoomx)
+           WRITE(lunout,*)" dzoomx = ",dzoomx
+
+           WRITE(lunout,*) ""
+           WRITE(lunout,*) "extension en latitude de la zone du zoom", 
+     & " (fraction de la zone totale)"
+           dzoomy=0. ! default value
+           call getin("dzoomy",dzoomy)
+           WRITE(lunout,*)" dzoomy = ",dzoomy
+
+           WRITE(lunout,*) ""
+           WRITE(lunout,*) "raideur du zoom en  X"
+           taux=2. ! default value
+           call getin("taux",taux)
+           WRITE(lunout,*)" taux = ",taux
+
+           WRITE(lunout,*) ""
+           WRITE(lunout,*) "raideur du zoom en  Y"
+           tauy=2.0 ! default value
+           call getin("tauy",tauy)
+           WRITE(lunout,*)" tauy = ",tauy
+
+        END IF ! of if (.not.etatinit )
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "Use a sponge layer?"
+        callsponge=.true. ! default value
+        call getin("callsponge",callsponge)
+        WRITE(lunout,*)" callsponge = ",callsponge
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "Sponge: number of layers over which",
+     &                    " sponge extends"
+        nsponge=3 ! default value
+        call getin("nsponge",nsponge)
+        WRITE(lunout,*)" nsponge = ",nsponge
+
+        WRITE(lunout,*)""
+        WRITE(lunout,*)"Sponge mode: (forcing is towards ..."
+        WRITE(lunout,*)"  over upper nsponge layers)"
+        WRITE(lunout,*)"  0: (h=hmean,u=v=0)"
+        WRITE(lunout,*)"  1: (h=hmean,u=umean,v=0)"
+        WRITE(lunout,*)"  2: (h=hmean,u=umean,v=vmean)"
+        mode_sponge=2 ! default value
+        call getin("mode_sponge",mode_sponge)
+        WRITE(lunout,*)" mode_sponge = ",mode_sponge
+
+        WRITE(lunout,*) ""
+        WRITE(lunout,*) "Sponge: characteristic time scale tetasponge"
+        WRITE(lunout,*) "(seconds) at topmost layer (time scale then "
+        WRITE(lunout,*) " doubles with decreasing layer index)."
+        tetasponge=50000.0
+        call getin("tetasponge",tetasponge)
+        WRITE(lunout,*)" tetasponge = ",tetasponge
+
+
+      WRITE(lunout,*) '-----------------------------------------------'
+      WRITE(lunout,*) ' '
+      WRITE(lunout,*) ' '
+c
+
+c       Unlike on Earth (cf LMDZ2.2) , always a regular grid on Mars :
+        ysinus = .false. !Mars Mettre a jour
+
+
+      WRITE(lunout,*) '-----------------------------------------------'
+      WRITE(lunout,*) ' '
+      WRITE(lunout,*) ' '
+cc
+      ELSE
+        write(tapeerr,*) ' WHERE IS run.def ? WE NEED IT !!!!!!!!!!!!!!'
+        stop
+      ENDIF ! of IF(ierr.eq.0)
+
+c     Test sur le zoom
+
+      if((grossismx.eq.1).and.(grossismy.eq.1)) then  
+c        Pas de zoom :
+         write(lunout,*) 'No zoom ? -> fxyhypb set to False'
+     &   ,'           (It runs better that way)'
+         fxyhypb = .false.
+      else     
+c        Avec Zoom
+         if (.not.fxyhypb) stop 'With zoom, fxyhypb should be set to T 
+     &in run.def for this version... -> STOP ! '     
+      end if
+#else
+      write(*,*) "defrun_new should not be used in parallel mode!"
+      stop
+#endif
+! of #ifndef CPP_PARA
+      END
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/exner_hyb.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/exner_hyb.F	(revision 2365)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../../dyn3d/exner_hyb.F
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/lect_start_archive.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/lect_start_archive.F	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/lect_start_archive.F	(revision 2366)
@@ -2,13 +2,10 @@
      &     date,tsurf,tsoil,emis,q2,
      &     t,ucov,vcov,ps,h,phisold_newgrid,
-     &     q,qsurf,surfith,nid)
-
-!      USE surfdat_h
+     &     q,qsurf,tankCH4,surfith,nid)
+
       USE comchem_h, only : cnames, nkim
       USE comchem_newstart_h
       USE comsoil_h, ONLY: nsoilmx, layer, mlayer, volcapa, inertiedat
       USE infotrac, ONLY: tname, nqtot
-!      USE control_mod
-! to use  'getin'
       USE callkeys_mod, only: callchim
       USE comvert_mod, ONLY: ap,bp,aps,bps,preff
@@ -17,14 +14,5 @@
 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    Routine to load variables from the "start_archive.nc" file
 c
 c=======================================================================
@@ -32,12 +20,9 @@
       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"
+      include "dimensions.h"
+      include "paramet.h"
+      include "comgeom2.h"
+      include "netcdf.inc"
+
 c=======================================================================
 c   Declarations
@@ -49,4 +34,5 @@
 c------------------------------------
       INTEGER   imold,jmold,lmold,nsoilold,nqold
+
 
 c Variables pour les lectures des fichiers "ini" 
@@ -60,20 +46,5 @@
       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
+      REAL,INTENT(OUT) :: date
       INTEGER   memo
 !      character (len=50) :: tmpname
@@ -81,14 +52,7 @@
 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)
+      REAL,INTENT(OUT) :: vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
+      REAL,INTENT(OUT) :: h(iip1,jjp1,llm),ps(iip1,jjp1)
+      REAL,INTENT(OUT) :: q(iip1,jjp1,llm,nqtot)
 
 c Physique sur grille scalaire 
@@ -97,12 +61,14 @@
 c variable physique
 c------------------
-      REAL tsurf(ngrid) ! surface temperature
-      REAL tsoil(ngrid,nsoilmx) ! soil temperature
-      REAL emis(ngrid)
-      REAL q2(ngrid,llm+1),qsurf(ngrid,nqtot)
+      REAL,INTENT(OUT) :: tsurf(ngrid) ! surface temperature
+      REAL,INTENT(OUT) :: tsoil(ngrid,nsoilmx) ! soil temperature
+      REAL,INTENT(OUT) :: emis(ngrid)
+      REAL,INTENT(OUT) :: q2(ngrid,llm+1),qsurf(ngrid,nqtot)
+      REAL,INTENT(OUT) :: tankCH4(ngrid)
 c     REAL phisfi(ngrid)
 
       INTEGER i,j,l
-      INTEGER nid,nvarid
+      INTEGER,INTENT(IN) :: nid
+      INTEGER :: nvarid
 c     REAL year_day,periheli,aphelie,peri_day
 c     REAL obliquit,z0,emin_turb,lmixmin
@@ -116,11 +82,12 @@
 c------------------------------------------------------
       real us(iip1,jjp1,llm),vs(iip1,jjp1,llm)
-      REAL phisold_newgrid(iip1,jjp1)
-      REAL t(iip1,jjp1,llm)
+      REAL,INTENT(OUT) :: phisold_newgrid(iip1,jjp1)
+      REAL,INTENT(OUT) :: t(iip1,jjp1,llm)
       real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx)
       real inertiedatS(iip1,jjp1,nsoilmx)
       real emisS(iip1,jjp1)
       REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot)
-      
+      REAL tankCH4S(iip1,jjp1)      
+ 
       real ptotal
 
@@ -148,4 +115,5 @@
       real, dimension(:,:), allocatable :: tsurfold
       real, dimension(:,:), allocatable :: emisold
+      real, dimension(:,:), allocatable :: tankCH4old
       real, dimension(:,:,:,:), allocatable :: qold
       
@@ -167,5 +135,5 @@
       real, dimension(:), allocatable :: newval
 
-      real surfith(iip1,jjp1) ! surface thermal inertia
+      real,intent(out) :: surfith(iip1,jjp1) ! surface thermal inertia
       ! surface thermal inertia at old horizontal grid resolution
       real, dimension(:,:), allocatable :: surfithold 
@@ -293,4 +261,5 @@
       allocate(tsurfold(imold+1,jmold+1))
       allocate(emisold(imold+1,jmold+1))
+      allocate(tankCH4old(imold+1,jmold+1))
       allocate(q2old(imold+1,jmold+1,lmold+1))
 !      allocate(tsoilold(imold+1,jmold+1,nsoilmx))
@@ -312,7 +281,7 @@
       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
+      write(*,*) 'lect_start_archive: q2',ngrid,llm+1
+      write(*,*) 'lect_start_archive: q2S',iip1,jjp1,llm+1
+      write(*,*) 'lect_start_archive: q2old',imold+1,jmold+1,lmold+1
 
 !-----------------------------------------------------------------------
@@ -327,5 +296,5 @@
       ierr = NF_INQ_VARID (nid, "controle", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "Lect_start_archive: champ <controle> est absent"
+         PRINT*, "Lect_start_archive: Field <controle> not found"
          CALL abort
       ENDIF
@@ -348,5 +317,5 @@
       ierr = NF_INQ_VARID (nid, "rlonv", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <rlonv> est absent"
+         PRINT*, "lect_start_archive: Field <rlonv> not found"
          CALL abort
       ENDIF
@@ -357,5 +326,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <rlonv>"
+         PRINT*, "lect_start_archive: Failed loading <rlonv>"
          CALL abort
       ENDIF
@@ -363,5 +332,5 @@
       ierr = NF_INQ_VARID (nid, "rlatu", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <rlatu> est absent"
+         PRINT*, "lect_start_archive: Field <rlatu> not found"
          CALL abort
       ENDIF 
@@ -372,5 +341,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <rlatu>"
+         PRINT*, "lect_start_archive: Failed loading <rlatu>"
          CALL abort
       ENDIF
@@ -378,5 +347,5 @@
       ierr = NF_INQ_VARID (nid, "rlonu", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <rlonu> est absent"
+         PRINT*, "lect_start_archive: Field <rlonu> not found"
          CALL abort
       ENDIF
@@ -387,5 +356,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <rlonu>"
+         PRINT*, "lect_start_archive: Failed loading <rlonu>"
          CALL abort
       ENDIF
@@ -393,5 +362,5 @@
       ierr = NF_INQ_VARID (nid, "rlatv", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <rlatv> est absent"
+         PRINT*, "lect_start_archive: Field <rlatv> not found"
          CALL abort
       ENDIF
@@ -402,5 +371,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <rlatv>"
+         PRINT*, "lect_start_archive: Failed loading <rlatv>"
          CALL abort
       ENDIF
@@ -413,5 +382,5 @@
       ierr = NF_INQ_VARID (nid, "aps", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <aps> est absent"
+         PRINT*, "lect_start_archive: Field <aps> not found"
          apsold=0
          PRINT*, "<aps> set to 0"
@@ -423,5 +392,5 @@
 #endif
          IF (ierr .NE. NF_NOERR) THEN
-            PRINT*, "lect_start_archive: Lecture echouee pour <aps>"
+            PRINT*, "lect_start_archive: Failed loading <aps>"
          ENDIF
       ENDIF
@@ -429,5 +398,5 @@
       ierr = NF_INQ_VARID (nid, "bps", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <bps> est absent"
+         PRINT*, "lect_start_archive: Field <bps> not found"
          PRINT*, "It must be an old start_archive, lets look for sig_s"
          ierr = NF_INQ_VARID (nid, "sig_s", nvarid)
@@ -443,5 +412,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <bps>"
+         PRINT*, "lect_start_archive: Failed loading <bps>"
          CALL abort
       END IF
@@ -532,5 +501,5 @@
       ierr=NF_INQ_VARID(nid,"preskim",nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <preskim> est absent"
+         PRINT*, "lect_start_archive: Field <preskim> not found"
          CALL abort
       ENDIF
@@ -541,5 +510,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <preskim>"
+         PRINT*, "lect_start_archive: Failed reading <preskim>"
          CALL abort
       ENDIF
@@ -551,5 +520,5 @@
       ierr = NF_INQ_VARID (nid, "phisinit", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <phisinit> est absent"
+         PRINT*, "lect_start_archive: Field <phisinit> not found"
          CALL abort
       ENDIF
@@ -560,5 +529,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <phisinit>"
+         PRINT*, "lect_start_archive: Failed loading <phisinit>"
          CALL abort
       ENDIF
@@ -579,5 +548,5 @@
          ierr = NF_INQ_DIMID (nid, "temps", nvarid)
          IF (ierr .NE. NF_NOERR) THEN
-            PRINT*, "lect_start_archive: Le champ <Time> est absent"
+            PRINT*, "lect_start_archive: Field <Time> not found"
             CALL abort
          endif
@@ -594,5 +563,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <Time>"
+         PRINT*, "lect_start_archive: Failed loading <Time>"
          CALL abort
       ENDIF
@@ -600,11 +569,11 @@
       write(*,*)
       write(*,*)
-      write(*,*) 'Differentes dates des etats initiaux stockes:'
-      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+      write(*,*) 'Available dates for the stored initial conditions:'
+      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))
+        write(*,*) 'initial state for day ' ,int(timelist(i))
 c       write(*,6) nint(timelist(i)),nint(mod(timelist(i),669)),
 c    .    sollong(i)
@@ -614,5 +583,5 @@
  
       write(*,*)
-      write(*,*) 'Choix de la date'
+      write(*,*) 'Choice for the date'
  123  read(*,*,iostat=ierr) date
       if(ierr.ne.0) goto 123
@@ -627,10 +596,10 @@
         write(*,*)
         write(*,*)
-        write(*,*) 'He alors... Y sait pas lire !?!'
+        write(*,*) "Wrong value... can't you read !?!"
         write(*,*)
-        write(*,*) 'Differentes dates des etats initiaux stockes:'
-        write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+        write(*,*) 'Available dates for the stored initial conditions:'
+        write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
         do i=1,timelen
-          write(*,*) 'etat initial au jour martien' ,nint(timelist(i))
+          write(*,*) 'initial state for day ' ,nint(timelist(i))
 c         write(*,6) nint(timelist(i)),nint(mod(timelist(i),669))
         end do
@@ -644,5 +613,5 @@
 
 c-----------------------------------------------------------------------
-c 5.1 Lecture des champs 2D (emis,ps,tsurf,Tg[10], qsurf)
+c 5.1 Lecture des champs 2D (emis,ps,tsurf,Tg[10], qsurf,tankCH4)
 c-----------------------------------------------------------------------
  
@@ -650,8 +619,7 @@
       count=(/imold+1,jmold+1,1,0/)
        
-
       ierr = NF_INQ_VARID (nid, "emis", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <emis> est absent"
+         PRINT*, "lect_start_archive: Field <emis> not found"
          CALL abort
       ENDIF
@@ -662,5 +630,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <emis>"
+         PRINT*, "lect_start_archive: Failed loading <emis>"
          CALL abort
       ENDIF
@@ -668,5 +636,5 @@
       ierr = NF_INQ_VARID (nid, "ps", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <ps> est absent"
+         PRINT*, "lect_start_archive: Field <ps> not found"
          CALL abort
       ENDIF
@@ -677,5 +645,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <ps>"
+         PRINT*, "lect_start_archive: Failed loading <ps>"
          CALL abort
       ENDIF
@@ -683,5 +651,5 @@
       ierr = NF_INQ_VARID (nid, "tsurf", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <tsurf> est absent"
+         PRINT*, "lect_start_archive: Field <tsurf> not found"
          CALL abort
       ENDIF
@@ -692,5 +660,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <tsurf>"
+         PRINT*, "lect_start_archive: Failed loading <tsurf>"
          CALL abort
       ENDIF
@@ -698,5 +666,5 @@
       ierr = NF_INQ_VARID (nid, "q2surf", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <q2surf> est absent"
+         PRINT*, "lect_start_archive: Field <q2surf> not found"
          CALL abort
       ENDIF
@@ -707,5 +675,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <q2surf>"
+         PRINT*, "lect_start_archive: Failed loading <q2surf>"
          CALL abort
       ENDIF
@@ -721,10 +689,4 @@
         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
@@ -751,9 +713,22 @@
      &             " 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
-
+c
+      ierr = NF_INQ_VARID (nid, "tankCH4", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Field <tankCH4> not found"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,tankCH4old)
+#else
+      ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,tqnkCH4old)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Failed loading <tankCH4>"
+         CALL abort
+      ENDIF
 
 !-----------------------------------------------------------------------
@@ -813,22 +788,4 @@
       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 Read 3D upper chemistry fields, if needed
@@ -847,8 +804,8 @@
       IF (ierr .NE. NF_NOERR) THEN ! H_up not found
       
-         PRINT*, "lect_start_archive: Le champ <H_up> est absent..."        
+         PRINT*, "lect_start_archive: Field <H_up> not found..."        
          IF (callchim) THEN
-           PRINT*, "... mais callchim=.TRUE. dans callphys.def !"
-           PRINT*, "Verifiez start_archive.nc ou desactivez callchim !"
+           PRINT*, "... but callchim=.TRUE. in callphys.def !"
+           PRINT*, "Check start_archive.nc or deactivate callchim !"
            CALL abort
          ELSE
@@ -861,7 +818,7 @@
       
         IF (.not.callchim) THEN
-          PRINT*, "lect_start_archive: Le champ <H_up> est present..."
-          PRINT*, "... mais callchim=.FALSE. dans callphys.def !"
-          PRINT*, "Si vous voulez gerer la chimie activez callchim !"
+          PRINT*, "lect_start_archive: Field <H_up> found..."
+          PRINT*, "... but callchim=.FALSE. in callphys.def !"
+          PRINT*, "If you want to have chemistry, activate callchim !"
           ! CALL abort ! This is too violent to abort here we can start from an archive with chemistry and don't want to use it - JVO !
         ELSE
@@ -874,5 +831,5 @@
 #endif
           IF (ierr .NE. NF_NOERR) THEN
-             PRINT*, "lect_start_archive: Lecture echouee pour <H_up>"
+             PRINT*, "lect_start_archive: Failed reading <H_up>"
              CALL abort
           ENDIF
@@ -882,6 +839,6 @@
             ierr=NF_INQ_VARID(nid,trim(cnames(iq))//"_up",nvarid)
             IF (ierr .NE. NF_NOERR) THEN
-              PRINT*, "lect_start_archive: Le champ <"
-     &//trim(cnames(iq))//"_up> est absent..."
+              PRINT*, "lect_start_archive: Field <"
+     &//trim(cnames(iq))//"_up> not found..."
               CALL abort
             ENDIF
@@ -894,5 +851,5 @@
 #endif
             IF (ierr .NE. NF_NOERR) THEN
-               PRINT*, "lect_start_archive: Lecture echouee pour <"
+               PRINT*, "lect_start_archive: Failed reading <"
      &//trim(cnames(iq))//"_up>"
                CALL abort
@@ -916,5 +873,5 @@
       ierr = NF_INQ_VARID (nid,"temp", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <temp> est absent"
+         PRINT*, "lect_start_archive: Field <temp> not found"
          CALL abort
       ENDIF
@@ -925,5 +882,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <temp>"
+         PRINT*, "lect_start_archive: Failed loading <temp>"
          CALL abort
       ENDIF
@@ -931,5 +888,5 @@
       ierr = NF_INQ_VARID (nid,"u", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <u> est absent"
+         PRINT*, "lect_start_archive: Field <u> not found"
          CALL abort
       ENDIF
@@ -940,5 +897,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <u>"
+         PRINT*, "lect_start_archive: Failed loading <u>"
          CALL abort
       ENDIF
@@ -946,5 +903,5 @@
       ierr = NF_INQ_VARID (nid,"v", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <v> est absent"
+         PRINT*, "lect_start_archive: Field <v> not found"
          CALL abort
       ENDIF
@@ -955,5 +912,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <v>"
+         PRINT*, "lect_start_archive: Failed loading <v>"
          CALL abort
       ENDIF
@@ -961,5 +918,5 @@
       ierr = NF_INQ_VARID (nid,"q2atm", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Le champ <q2atm> est absent"
+         PRINT*, "lect_start_archive: Field <q2atm> not found"
          CALL abort
       ENDIF
@@ -970,5 +927,5 @@
 #endif
       IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "lect_start_archive: Lecture echouee pour <q2atm>"
+         PRINT*, "lect_start_archive: Failed loading <q2atm>"
          CALL abort
       ENDIF
@@ -1071,6 +1028,6 @@
 
       write(*,*)
-      write(*,*)'Ancienne grille: masse de l atm :',ptotalold
-      write(*,*)'Nouvelle grille: masse de l atm :',ptotal
+      write(*,*)'Old grid: atmospheric mass :',ptotalold
+      write(*,*)'New grid: atmospheric mass :',ptotal
       write (*,*) 'Ratio new atm./ old atm =', ptotal/ptotalold 
       write(*,*)
@@ -1429,5 +1386,3 @@
       deallocate(ykim_upoldS)
 
-!      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 2365)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/newstart.F	(revision 2366)
@@ -39,4 +39,5 @@
       use iniphysiq_mod, only: iniphysiq
       use phyetat0_mod, only: phyetat0
+      use exner_hyb_m, only: exner_hyb
       use tracer_h
       implicit none
@@ -125,5 +126,4 @@
       REAL :: xpn,xps,xppn(iim),xpps(iim)
       REAL :: p3d(iip1, jjp1, llm+1)
-      REAL :: beta(iip1,jjp1,llm)
 !      REAL dteta(ip1jmp1,llm)
 
@@ -542,5 +542,5 @@
      &   date,tsurf,tsoil,emis,q2,
      &   t,ucov,vcov,ps,teta,phisold_newgrid,
-     &   q,qsurf,surfith,nid)
+     &   q,qsurf,tankCH4,surfith,nid)
         write(*,*) "OK, read start_archive file"
         ! copy soil thermal inertia
@@ -1120,5 +1120,5 @@
 c-----------------------------------------------------------------------
 
-      CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
+      CALL exner_hyb(ip1jmp1, ps, p3d, pks, pk, pkf)
 ! Calcul de la temperature potentielle teta
 
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/start2archive.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/start2archive.F	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/start2archive.F	(revision 2366)
@@ -22,9 +22,5 @@
       USE comsoil_h
       USE comchem_h, only : cnames, nkim, nlaykim_up, preskim, ykim_up
-!      USE comgeomfi_h, ONLY: lati, long, area
-!      use control_mod
-!      use comgeomphy, only: initcomgeomphy
-! to use  'getin'
-      USE ioipsl_getincom
+      USE ioipsl_getincom, only: getin
       USE planete_mod, only: year_day
       USE mod_const_mpi, ONLY: COMM_LMDZ
@@ -35,6 +31,8 @@
       USE temps_mod, ONLY: day_ini
       USE iniphysiq_mod, ONLY: iniphysiq
+      use phys_state_var_mod, only: phys_state_var_init
       use phyetat0_mod, only: phyetat0
       use tracer_h
+      use exner_hyb_m, only: exner_hyb
       implicit none
 
@@ -44,9 +42,5 @@
       include "comdissip.h"
       include "comgeom.h"
-!#include "control.h"
-
-!#include "dimphys.h"
-!#include "planete.h"
-!#include"advtrac.h"
+
       include "netcdf.inc"
 c-----------------------------------------------------------------------
@@ -62,5 +56,4 @@
       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
@@ -262,4 +255,7 @@
       Lmodif=0
 
+! Allocate saved arrays (as in firstcall of physiq)
+      call phys_state_var_init(nqtot)
+      
 ! Initialize tracer names, indexes and properties
       CALL initracer2(nqtot,tname)
@@ -312,5 +308,5 @@
 
       CALL pression(ip1jmp1, ap, bp, ps, p3d)
-      call exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
+      call exner_hyb(ip1jmp1, ps, p3d, pks, pk, pkf)
 
 c=======================================================================
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/write_archive.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/write_archive.F	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/write_archive.F	(revision 2366)
@@ -37,8 +37,8 @@
       implicit none
 
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-#include "netcdf.inc"
+      include "dimensions.h"
+      include "paramet.h"
+      include "comgeom.h"
+      include "netcdf.inc"
 
 c-----------------------------------------------------------------------
@@ -83,5 +83,5 @@
 
               write (*,*) "====================="
-              write (*,*) "creation de ",nom
+              write (*,*) "defining ",nom
               call def_var(nid,nom,titre,unite,4,id,varid,ierr)
 
@@ -113,6 +113,6 @@
 
            if (ierr.ne.NF_NOERR) then
-              write(*,*) "***** PUT_VAR matter in write_archive"
-              write(*,*) "***** with ",nom," ",nf_STRERROR(ierr)
+              write(*,*) "***** PUT_VAR problem in write_archive"
+              write(*,*) "***** with ",trim(nom)," ",nf_STRERROR(ierr)
               call abort
            endif
@@ -140,5 +140,5 @@
 	  ! define the variable
 	  write(*,*)"====================="
-	  write(*,*)"defining ",nom
+	  write(*,*)"defining variable ",trim(nom)
 	  call def_var(nid,nom,titre,unite,4,id,varid,ierr)
 	  
@@ -183,5 +183,5 @@
           ! define the variable
           write(*,*)"====================="
-          write(*,*)"defining ",nom
+          write(*,*)"defining variable ",trim(nom)
           call def_var(nid,nom,titre,unite,4,id,varid,ierr)
 
@@ -203,4 +203,5 @@
            ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,px)
 #endif
+
 
 
@@ -222,5 +223,5 @@
 
               write (*,*) "====================="
-              write (*,*) "creation de ",nom
+              write (*,*) "defining variable ",trim(nom)
 
               call def_var(nid,nom,titre,unite,3,id,varid,ierr)
@@ -243,5 +244,5 @@
 
            if (ierr.ne.NF_NOERR) then
-              write(*,*) "***** PUT_VAR matter in write_archive"
+              write(*,*) "***** PUT_VAR problem in write_archive"
               write(*,*) "***** with ",nom,nf_STRERROR(ierr)
               call abort
@@ -264,5 +265,5 @@
 
               write (*,*) "====================="
-              write (*,*) "creation de ",nom
+              write (*,*) "defining variable ",trim(nom)
 
               call def_var(nid,nom,titre,unite,1,id,varid,ierr)
@@ -279,5 +280,5 @@
 #endif
            if (ierr.ne.NF_NOERR) then
-              write(*,*) "***** PUT_VAR matter in write_archive"
+              write(*,*) "***** PUT_VAR problem in write_archive"
               write(*,*) "***** with ",nom,nf_STRERROR(ierr)
               call abort
@@ -289,5 +290,4 @@
         endif ! of if (dim.eq.3) else if (dim.eq.-3) ....
 
-      return
       end
 
Index: trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90	(revision 2366)
@@ -17,5 +17,6 @@
       use comcstfi_mod, only: pi, mugaz, cpp
       use callkeys_mod, only: diurnal,tracer,seashaze,corrk_recombin,   &
-                              strictboundcorrk,specOLR,diagdtau
+                              strictboundcorrk,specOLR,diagdtau,        &
+                              tplanckmin,tplanckmax
       use geometry_mod, only: latitude
 
@@ -115,4 +116,8 @@
       real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
       real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
+
+      ! Miscellaneous :
+      character(len=100) :: message
+      character(len=10),parameter :: subname="callcorrk"
 
       logical OLRz
@@ -377,5 +382,23 @@
             endif
          endif
+
+         if (tlevrad(k).lt.tplanckmin) then
+            print*,'Minimum temperature is outside the boundaries for'
+            print*,'Planck function integration set in callphys.def, aborting.'
+            print*,"k=",k," tlevrad(k)=",tlevrad(k)
+            print*,"tplanckmin=",tplanckmin
+            message="Minimum temperature outside Planck function bounds - Change tplanckmin in callphys.def"
+            call abort_physic(subname,message,1)
+          else if (tlevrad(k).gt.tplanckmax) then
+            print*,'Maximum temperature is outside the boundaries for'
+            print*,'Planck function integration set in callphys.def, aborting.'
+            print*,"k=",k," tlevrad(k)=",tlevrad(k)
+            print*,"tplanckmax=",tplanckmax
+            message="Maximum temperature outside Planck function bounds - Change tplanckmax in callphys.def"
+            call abort_physic(subname,message,1)
+          endif
+
       enddo
+
       do k=1,L_NLAYRAD+1
          if(tmid(k).lt.tgasmin)then
Index: trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 2366)
@@ -60,4 +60,9 @@
       real,save :: szangle
 !$OMP THREADPRIVATE(szangle)
+
+      real,save :: tplanckmin
+      real,save :: tplanckmax
+      real,save :: dtplanck
+!$OMP THREADPRIVATE(tplanckmin,tplanckmax,dtplanck)
       real,save :: Fat1AU
       real,save :: stelTbb
@@ -76,5 +81,6 @@
       real,save :: surfemis
 !$OMP THREADPRIVATE(surfalbedo,surfemis)
-
+      real,save :: noseason_day
+!$OMP THREADPRIVATE(noseason_day)
       logical,save :: iscallphys=.false.!existence of callphys.def
 !$OMP THREADPRIVATE(iscallphys)
Index: trunk/LMDZ.TITAN/libf/phytitan/comsaison_h.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/comsaison_h.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/comsaison_h.F90	(revision 2366)
@@ -4,10 +4,13 @@
        implicit none 
 
-       integer isaison
-       logical callsais
-       real dist_star,declin,right_ascen
+!       integer,save :: isaison
+!       logical,save :: callsais
+!!$OMP THREADPRIVATE(isaison,callsais)
+
+       real,save :: dist_star,declin,right_ascen
+!$OMP THREADPRIVATE(dist_star,declin,right_ascen)
 
        real, allocatable, dimension(:) :: mu0,fract
-!$OMP THREADPRIVATE(isaison,callsais,dist_star,declin,mu0,fract)
+!$OMP THREADPRIVATE(mu0,fract)
 
        end module comsaison_h
Index: trunk/LMDZ.TITAN/libf/phytitan/dyn1d/abort_gcm.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/dyn1d/abort_gcm.F	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/dyn1d/abort_gcm.F	(revision 2366)
@@ -1,1 +1,1 @@
-link ../../dyn3d/abort_gcm.F
+link ../../../../LMDZ.COMMON/libf/dyn3d/abort_gcm.F
Index: trunk/LMDZ.TITAN/libf/phytitan/dyn1d/comconst_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/dyn1d/comconst_mod.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/dyn1d/comconst_mod.F90	(revision 2366)
@@ -1,1 +1,1 @@
-link ../../dyn3d/comconst_mod.F90
+link ../../../../LMDZ.COMMON/libf/dyn3d_common/comconst_mod.F90
Index: trunk/LMDZ.TITAN/libf/phytitan/dyn1d/comgeom.h
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/dyn1d/comgeom.h	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/dyn1d/comgeom.h	(revision 2366)
@@ -1,1 +1,1 @@
-link ../../dyn3d/comgeom.h
+link ../../../../LMDZ.COMMON/libf/dyn3d_common/comgeom.h
Index: trunk/LMDZ.TITAN/libf/phytitan/dyn1d/comvert_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/dyn1d/comvert_mod.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/dyn1d/comvert_mod.F90	(revision 2366)
@@ -1,1 +1,1 @@
-link ../../dyn3d/comvert_mod.F90
+link ../../../../LMDZ.COMMON/libf/dyn3d_common/comvert_mod.F90
Index: trunk/LMDZ.TITAN/libf/phytitan/dyn1d/control_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/dyn1d/control_mod.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/dyn1d/control_mod.F90	(revision 2366)
@@ -1,1 +1,1 @@
-link ../../dyn3d/control_mod.F90
+link ../../../../LMDZ.COMMON/libf/dyn3d_common/control_mod.F90
Index: trunk/LMDZ.TITAN/libf/phytitan/dyn1d/disvert.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/dyn1d/disvert.F	(revision 2365)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../../dyn3d/disvert.F
Index: trunk/LMDZ.TITAN/libf/phytitan/dyn1d/disvert_noterre.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/dyn1d/disvert_noterre.F	(revision 2366)
+++ trunk/LMDZ.TITAN/libf/phytitan/dyn1d/disvert_noterre.F	(revision 2366)
@@ -0,0 +1,1 @@
+link ../../../../LMDZ.COMMON/libf/dyn3d_common/disvert_noterre.F
Index: trunk/LMDZ.TITAN/libf/phytitan/dyn1d/ener_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/dyn1d/ener_mod.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/dyn1d/ener_mod.F90	(revision 2366)
@@ -1,1 +1,1 @@
-link ../../dyn3d/ener_mod.F90
+link ../../../../LMDZ.COMMON/libf/dyn3d_common/ener_mod.F90
Index: trunk/LMDZ.TITAN/libf/phytitan/dyn1d/infotrac.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/dyn1d/infotrac.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/dyn1d/infotrac.F90	(revision 2366)
@@ -1,1 +1,1 @@
-link ../../dyn3d/infotrac.F90
+link ../../../../LMDZ.COMMON/libf/dyn3d_common/infotrac.F90
Index: trunk/LMDZ.TITAN/libf/phytitan/dyn1d/logic_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/dyn1d/logic_mod.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/dyn1d/logic_mod.F90	(revision 2366)
@@ -1,1 +1,1 @@
-link ../../dyn3d/logic_mod.F90
+link ../../../../LMDZ.COMMON/libf/dyn3d/logic_mod.F90
Index: trunk/LMDZ.TITAN/libf/phytitan/dyn1d/mod_const_mpi.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/dyn1d/mod_const_mpi.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/dyn1d/mod_const_mpi.F90	(revision 2366)
@@ -1,1 +1,1 @@
-link ../../dyn3d/mod_const_mpi.F90
+link ../../../../LMDZ.COMMON/libf/dyn3d/mod_const_mpi.F90
Index: trunk/LMDZ.TITAN/libf/phytitan/dyn1d/paramet.h
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/dyn1d/paramet.h	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/dyn1d/paramet.h	(revision 2366)
@@ -1,1 +1,1 @@
-link ../../dyn3d/paramet.h
+link ../../../../LMDZ.COMMON/libf/dyn3d_common/paramet.h
Index: trunk/LMDZ.TITAN/libf/phytitan/dyn1d/rcm1d.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/dyn1d/rcm1d.F	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/dyn1d/rcm1d.F	(revision 2366)
@@ -28,5 +28,5 @@
      &                       presnivs,pseudoalt,scaleheight
       USE vertical_layers_mod, ONLY: init_vertical_layers
-      USE logic_mod, ONLY: hybrid,autozlevs
+      USE logic_mod, ONLY: hybrid
       use regular_lonlat_mod, only: init_regular_lonlat
       use planete_mod, only: ini_planete_mod
@@ -35,4 +35,5 @@
       use mod_interface_dyn_phys, only: init_interface_dyn_phys
       use inifis_mod, only: inifis
+      use phys_state_var_mod, only: phys_state_var_init
       use physiq_mod, only: physiq
       implicit none
@@ -118,4 +119,5 @@
 
 !     added by RW for autozlevs computation
+      logical autozlevs
       real nu, xx, pMIN, zlev, Htop
       real logplevs(llm)
@@ -132,23 +134,20 @@
 c INITIALISATION
 c=======================================================================
-! initialize "serial/parallel" related stuff
-!      CALL init_phys_lmdz(iim,jjm+1,llm,1,(/(jjm-1)*iim+2/))
-!      CALL init_phys_lmdz(1,1,llm,1,(/1/))
-!      call initcomgeomphy
-
-      !! those are defined in surfdat_h.F90
-      IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(1))
-      IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(1))
-      IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(1))
-      IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(1))
-      IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(1))
-      IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(1))
-      IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(1))
-      !! those are defined in comdiurn_h.F90
-      IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(1))
-      IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(1))
-      IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(1))
-      IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(1))
-
+
+      ! read nq from traceur.def
+      open(90,file='traceur.def',status='old',form='formatted',
+     &       iostat=ierr)
+      if (ierr.eq.0) then
+        read(90,*,iostat=ierr) nq
+      else
+        nq=0
+      endif
+      close(90)
+      
+      ! Initialize dimphy module
+      call init_dimphy(1,llm) 
+      ! now initialize arrays using phys_state_var_init
+      call phys_state_var_init(nq)
+      
       saveprofile=.false.
       saveprofile=.true.
@@ -480,5 +479,5 @@
 !      call init_vertical_layers(nlayer,preff,scaleheight,
 !     &                      ap,bp,aps,bps,presnivs,pseudoalt)
-      call init_dimphy(1,nlayer) ! Initialize dimphy module
+!      call init_dimphy(1,nlayer) ! Initialize dimphy module
       call ini_planete_mod(nlayer,preff,ap,bp)
 
@@ -643,5 +642,5 @@
       endif
 
-      call disvert
+      call disvert_noterre
       ! now that disvert has been called, initialize module vertical_layers_mod
       call init_vertical_layers(nlayer,preff,scaleheight,
@@ -669,4 +668,6 @@
          ENDDO
          
+
+
          DO ilayer=1,nlayer
 !     zlay(ilayer)=-300.E+0 *r*log(play(ilayer)/plev(1))
Index: trunk/LMDZ.TITAN/libf/phytitan/dyn1d/serre_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/dyn1d/serre_mod.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/dyn1d/serre_mod.F90	(revision 2366)
@@ -1,1 +1,1 @@
-link ../../dyn3d/serre_mod.F90
+link ../../../../LMDZ.COMMON/libf/dyn3d_common/serre_mod.F90
Index: trunk/LMDZ.TITAN/libf/phytitan/gfluxi.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/gfluxi.F	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/gfluxi.F	(revision 2366)
@@ -85,6 +85,6 @@
          LAMDA(L) = ALPHA(L)*(1.0D0-W0(L)*COSBAR(L))/UBARI
          
-         NT    = int(TLEV(2*L)*NTfac)   - NTstar+1
-         NT2   = int(TLEV(2*L+2)*NTfac) - NTstar+1
+         NT    = int(TLEV(2*L)*NTfac)   - NTstart+1
+         NT2   = int(TLEV(2*L+2)*NTfac) - NTstart+1
          
 ! AB : PLANCKIR(NW,NT) is replaced by P1, the linear interpolation result for a temperature NT
@@ -112,13 +112,13 @@
       ! -- same results for most thin atmospheres
       ! -- and stabilizes integrations
-      NT    = int(TLEV(2*L+1)*NTfac) - NTstar+1
+      NT    = int(TLEV(2*L+1)*NTfac) - NTstart+1
       !! For deep, opaque, thick first layers (e.g. Saturn)
       !! what is below works much better, not unstable, ...
       !! ... and actually fully accurate because 1st layer temp (JL) 
-      !NT    = int(TLEV(2*L)*NTfac) - NTstar+1
+      !NT    = int(TLEV(2*L)*NTfac) - NTstart+1
       !! (or this one yields same results
-      !NT    = int( (TLEV(2*L)+TLEV(2*L+1))*0.5*NTfac ) - NTstar+1
-      
-      NT2   = int(TLEV(2*L)*NTfac)   - NTstar+1
+      !NT    = int( (TLEV(2*L)+TLEV(2*L+1))*0.5*NTfac ) - NTstart+1
+      
+      NT2   = int(TLEV(2*L)*NTfac)   - NTstart+1
       
 ! AB : PLANCKIR(NW,NT) is replaced by P1, the linear interpolation result for a temperature NT
Index: trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 2366)
@@ -11,5 +11,4 @@
   use init_print_control_mod, only: init_print_control
   use radinc_h, only: ini_radinc_h
-  use datafile_mod
   use comdiurn_h, only: sinlat, coslat, sinlon, coslon
   use comgeomfi_h, only: totarea, totarea_planet
@@ -59,4 +58,5 @@
 !   declarations:
 !   -------------
+  use datafile_mod
   use ioipsl_getin_p_mod, only: getin_p
   IMPLICIT NONE
@@ -90,5 +90,7 @@
 
   ! Initialize some "temporal and calendar" related variables
+#ifndef MESOSCALE
   CALL init_time(day_ini,pdaysec,nday,ptimestep)
+#endif
 
   ! read in some parameters from "run.def" for physics,
@@ -147,4 +149,9 @@
      call getin_p("season",season)
      write(*,*) " season = ",season
+     
+     write(*,*) "No seasonal cycle: initial day to lock the run during restart"
+     noseason_day=0.0 ! default value
+     call getin_p("noseason_day",noseason_day)
+     write(*,*) "noseason_day=", noseason_day
 
      write(*,*) "Tidally resonant rotation ?"
@@ -276,4 +283,19 @@
      write(*,*) "strictboundcorrk = ",strictboundcorrk
 
+     write(*,*) "Minimum atmospheric temperature for Planck function integration ?"
+     tplanckmin=30.0 ! default value
+     call getin_p("tplanckmin",tplanckmin)
+     write(*,*) " tplanckmin = ",tplanckmin
+ 
+     write(*,*) "Maximum atmospheric temperature for Planck function integration ?"
+     tplanckmax=1500.0 ! default value
+     call getin_p("tplanckmax",tplanckmax)
+     write(*,*) " tplanckmax = ",tplanckmax
+ 
+     write(*,*) "Temperature step for Planck function integration ?"
+     dtplanck=0.1 ! default value
+     call getin_p("dtplanck",dtplanck)
+     write(*,*) " dtplanck = ",dtplanck
+ 
      write(*,*) "call gaseous absorption in the visible bands?", &
                     "(matters only if callrad=T)"
@@ -548,9 +570,7 @@
      endif ! of if (force_cpp)
      
-     
      call su_gases(nlayer,tracer)     
-     call calc_cpp_mugaz
-     
-     
+     call calc_cpp_mugaz 
+
      PRINT*,'--------------------------------------------'
      PRINT*
@@ -579,4 +599,5 @@
 
   ! Initializations for comgeomfi_h
+#ifndef MESOSCALE
   totarea=SSUM(ngrid,parea,1)
   call planetwide_sumval(parea,totarea_planet)
@@ -594,11 +615,11 @@
      coslon(ig)=cos(plon(ig))
   ENDDO
-
+#endif
   ! initialize variables in radinc_h
-  call ini_radinc_h(nlayer)
+  call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
  
   ! allocate "comsoil_h" arrays
   call ini_comsoil_h(ngrid)
-      
+    
   END SUBROUTINE inifis
 
Index: trunk/LMDZ.TITAN/libf/phytitan/optcv.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optcv.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/optcv.F90	(revision 2366)
@@ -11,5 +11,4 @@
                           callclouds,callmufi,seashaze,uncoupl_optic_haze
   use tracer_h, only: nmicro,nice
-  use MMP_OPTICS
 
   implicit none
Index: trunk/LMDZ.TITAN/libf/phytitan/phys_state_var_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/phys_state_var_mod.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/phys_state_var_mod.F90	(revision 2366)
@@ -10,12 +10,18 @@
 ! Declaration des variables
       USE dimphy, only : klon,klev
+      use comchem_h, only: nkim
+      USE callkeys_mod, only: callchim
       USE comsoil_h, only : nsoilmx
       use comsaison_h, only: mu0, fract
       use radcommon_h, only: gzlat, gzlat_ig, Cmk 
       USE radinc_h, only : L_NSPECTI, L_NSPECTV,naerkind
+      use surfdat_h, only: phisfi, albedodat,  &
+                        zmea, zstd, zsig, zgam, zthe
       use turb_mod, only: q2,sensibFlux,wstar,ustar,tstar,hfmax_th,zmax_th
 
-      real,allocatable,dimension(:,:),save :: ztprevious ! Previous loop Atmospheric Temperature (K) ! Useful for Dynamical Heating calculation.
-      real,allocatable,dimension(:,:),save :: zuprevious ! Previous loop Zonal Wind (m.s-1)          ! Useful for Zonal Wind tendency calculation.
+      real,allocatable,dimension(:,:),save :: ztprevious ! Previous loop Atmospheric Temperature (K)
+! Useful for Dynamical Heating calculation.
+      real,allocatable,dimension(:,:),save :: zuprevious
+!$OMP THREADPRIVATE(ztprevious,zuprevious)
 
       real, dimension(:),allocatable,save ::  tsurf                ! Surface temperature (K).
@@ -40,5 +46,4 @@
       real,dimension(:),allocatable,save :: fluxgrd     ! Surface conduction flux (W.m-2).
       real,dimension(:,:),allocatable,save :: qsurf     ! Tracer on surface (e.g. kg.m-2).
-
 !$OMP THREADPRIVATE(emis,dtrad,fluxrad_sky,fluxrad,capcal,fluxgrd,qsurf,q2)
 
@@ -48,4 +53,5 @@
       real,dimension(:),allocatable,save :: fluxsurf_sw     ! Incident Short Wave (stellar) surface flux (W.m-2).
       real,dimension(:),allocatable,save :: fluxsurfabs_sw  ! Absorbed Short Wave (stellar) flux by the surface (W.m-2).
+!$OMP THREADPRIVATE(fluxsurf_lw,fluxsurf_sw,fluxsurfabs_sw)
 
       real,dimension(:),allocatable,save :: fluxtop_lw      ! Outgoing LW (IR) flux to space (W.m-2).
@@ -53,4 +59,5 @@
       real,dimension(:),allocatable,save :: fluxtop_dn      ! Incoming SW (stellar) radiation at the top of the atmosphere (W.m-2).
       real,dimension(:),allocatable,save :: fluxdyn         ! Horizontal heat transport by dynamics (W.m-2).
+!$OMP THREADPRIVATE(fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn)
 
       real,dimension(:,:),allocatable,save :: OLR_nu        ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1)).
@@ -59,12 +66,16 @@
       real,dimension(:,:),allocatable,save :: zdtsw         ! SW heating tendencies (K/s).
       !real,dimension(:),allocatable,save :: sensibFlux      ! Turbulent flux given by the atmosphere to the surface (W.m-2).
+!$OMP THREADPRIVATE(OLR_nu,OSR_nu,zdtlw,zdtsw)
+
       real,dimension(:,:,:),allocatable,save :: int_dtauv   ! VI optical thickness of layers within narrowbands for diags ().
       real,dimension(:,:,:),allocatable,save :: int_dtaui   ! IR optical thickness of layers within narrowbands for diags ().
-
-!$OMPTHREADPRIVATE(fluxsurf_lw,fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu,&        
-        !$OMP zdtlw,zdtsw,sensibFlux,int_dtauv,int_dtaui))
+!$OMP THREADPRIVATE(int_dtaui,int_dtauv) 
 
       real,allocatable,dimension(:,:),save :: qsurf_hist
 !$OMP THREADPRIVATE(qsurf_hist)
+
+      ! For chemistry
+      real, dimension(:,:,:), allocatable, save   :: dycchi  ! NB : Only for chem tracers. Saved since chemistry is not called every step.
+!$OMP THREADPRIVATE(dycchi)
  
 CONTAINS
@@ -119,4 +130,5 @@
         ALLOCATE(int_dtaui(klon,klev,L_NSPECTI))
         ALLOCATE(int_dtauv(klon,klev,L_NSPECTV))
+        allocate(dycchi(klon,klev,nkim))
         ! This is defined in comsaison_h
         ALLOCATE(mu0(klon))
@@ -126,4 +138,12 @@
         ALLOCATE(gzlat_ig(klev))
         ALLOCATE(Cmk(klev))
+        ! This is defined in surfdat_h
+        ALLOCATE(albedodat(klon))
+        ALLOCATE(phisfi(klon))
+        ALLOCATE(zmea(klon))
+        ALLOCATE(zstd(klon))
+        ALLOCATE(zsig(klon))
+        ALLOCATE(zgam(klon))
+        ALLOCATE(zthe(klon))
         ! This is defined in turb_mod
         allocate(wstar(klon))
@@ -171,4 +191,5 @@
         DEALLOCATE(int_dtaui)
         DEALLOCATE(int_dtauv)
+        DEALLOCATE(dycchi)
         DEALLOCATE(mu0)
         DEALLOCATE(fract)
@@ -176,4 +197,11 @@
         DEALLOCATE(gzlat_ig)
         DEALLOCATE(Cmk)
+        DEALLOCATE(phisfi)
+        DEALLOCATE(albedodat)
+        DEALLOCATE(zmea)
+        DEALLOCATE(zstd)
+        DEALLOCATE(zsig)
+        DEALLOCATE(zgam)
+        DEALLOCATE(zthe)
         deallocate(wstar)
         deallocate(ustar)
Index: trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 2366)
@@ -16,5 +16,4 @@
       use radinc_h, only : L_NSPECTI,L_NSPECTV
       use radcommon_h, only: sigma, gzlat, grav, BWNV
-      use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe
       use comchem_h, only: nkim, cnames, nlaykim_up, ykim_up, ykim_tot, botCH4
       use comdiurn_h, only: coslat, sinlat, coslon, sinlon
@@ -34,7 +33,9 @@
       use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
       use time_phylmdz_mod, only: daysec
+#ifndef MESOSCALE
       use logic_mod, only: moyzon_ch
       use moyzon_mod, only:  zphibar, zphisbar, zplevbar, zplaybar, &
                              zzlevbar, zzlaybar, ztfibar, zqfibar
+#endif
       use callkeys_mod
       use phys_state_var_mod
@@ -42,5 +43,4 @@
 #ifndef MESOSCALE
       use vertical_layers_mod, only: presnivs, pseudoalt
-      use ioipsl_getin_p_mod, only: getin_p
       use mod_phys_lmdz_omp_data, ONLY: is_omp_master
 #else
@@ -56,5 +56,4 @@
       use wxios, only: wxios_context_init, xios_context_finalize
 #endif
-      use MMP_OPTICS
       use muphy_diag
       implicit none
@@ -177,5 +176,4 @@
 !   -------
 
-
       integer,intent(in) :: ngrid             ! Number of atmospheric columns.
       integer,intent(in) :: nlayer            ! Number of atmospheric layers.
@@ -211,8 +209,7 @@
 ! Local saved variables:
 ! ----------------------
-
       integer,save :: day_ini                                      ! Initial date of the run (sol since Ls=0).
       integer,save :: icount                                       ! Counter of calls to physiq during the run.
-
+!$OMP THREADPRIVATE(day_ini,icount)
 
 ! Local variables :
@@ -289,5 +286,4 @@
 
 
-
 ! local variables for DIAGNOSTICS : (diagfi & stat)
 ! -------------------------------------------------
@@ -299,5 +295,4 @@
       real zdtdyn(ngrid,nlayer)                          ! Dynamical Heating (K/s).
       real zdudyn(ngrid,nlayer)                          ! Dynamical Zonal Wind tendency (m.s-2).
-!$OMP THREADPRIVATE(ztprevious,zuprevious)
 
       real zhorizwind(ngrid,nlayer) ! Horizontal Wind ( sqrt(u**+v*v))
@@ -329,6 +324,4 @@
       real tf, ntf    
 
-!$OMP THREADPRIVATE(qsurf_hist)
-   
       ! Miscellaneous :
       character(len=10) :: tmp1
@@ -336,4 +329,6 @@
       
       character*2 :: str2
+
+#ifndef MESOSCALE
 
 ! Local variables for Titan chemistry and microphysics
@@ -348,6 +343,4 @@
 
       ! Molar fraction tendencies ( chemistry, condensation and evaporation ) for tracers (mol/mol/s)
-      real, dimension(:,:,:), allocatable, save   :: dycchi         ! NB : Only for chem tracers. Saved since chemistry is not called every step.
-!$OMP THREADPRIVATE(dycchi)
       real, dimension(ngrid,nlayer,nq)            :: dyccond        ! Condensation rate. NB : for all tracers, as we want to use indx on it.
       real, dimension(ngrid,nlayer,nq)            :: dyccondbar     ! For 2D chemistry
@@ -360,7 +353,9 @@
       real :: i2e(ngrid,nlayer)      ! int 2 ext factor ( X.kg-1 -> X.m-3 for diags )
 
+#ifdef USE_QTEST
       real,save,dimension(:,:,:), allocatable :: tpq ! Tracers for decoupled microphysical tests ( temporary in 01/18 )
 !$OMP THREADPRIVATE(tpq)
       real,dimension(ngrid,nlayer,nq) :: dtpq ! (temporary in 01/18)
+#endif
 
       logical file_ok
@@ -385,4 +380,6 @@
       END SUBROUTINE calmufi
     END INTERFACE
+
+#endif
       
 !==================================================================================================
@@ -396,8 +393,9 @@
 ! --------------------------------
       if (firstcall) then
-#ifndef MESOSCALE
+
+#ifdef USE_QTEST
         allocate(tpq(ngrid,nlayer,nq))
         tpq(:,:,:) = pq(:,:,:)
-
+#endif
         ! Initialisation of nmicro as well as tracers names, indexes ...
         if (ngrid.ne.1) then ! Already done in rcm1d
@@ -405,6 +403,7 @@
         endif
 
-        call phys_state_var_init(nq)
-
+        ! Allocate saved arrays (except for 1D model, where this has already been done)
+#ifndef MESOSCALE
+        if (ngrid>1) call phys_state_var_init(nq)
 #endif
 
@@ -453,4 +452,5 @@
          endif
 
+#ifndef MESOSCALE
 !        Initialize names and timestep for chemistry
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -463,7 +463,5 @@
                stop
             endif
-
-            allocate(dycchi(ngrid,nlayer,nkim)) ! only for chemical tracers
-            
+           
             ! Chemistry timestep
             ctimestep = ptimestep*REAL(ichim)
@@ -482,4 +480,5 @@
 
          ENDIF
+#endif
 
 #ifdef CPP_XIOS
@@ -492,5 +491,6 @@
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 #ifndef MESOSCALE
-         call phyetat0(startphy_file,ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
+         call phyetat0(startphy_file,                                 &
+                       ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
                        day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,tankCH4)                       
 #else
@@ -498,4 +498,5 @@
          q2(:,:)=0.0
          qsurf(:,:)=0.0
+         tankCH4(:)=0.0
          day_ini = pday
 #endif
@@ -534,4 +535,5 @@
          call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
 
+
          if(tlocked)then
             print*,'Planet is tidally locked at resonance n=',nres
@@ -576,11 +578,10 @@
 
 #ifndef MESOSCALE
-         if (ngrid.ne.1) then
-            ! Note : no need to create a restart file in 1d.
+         if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
             call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
-                          ptimestep,pday+nday,time_phys,cell_area,          &
-                          albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
-         endif
-#endif
+                         ptimestep,pday+nday,time_phys,cell_area,          &
+                         albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
+         endif
+#endif         
          
          ! XIOS outputs
@@ -624,5 +625,5 @@
          call stellarlong(zday,zls)
       else
-         call stellarlong(float(day_ini),zls)
+         call stellarlong(noseason_day,zls)
       end if
 
@@ -686,4 +687,5 @@
       ! JVO 19 : We shall always have correct altitudes in chemistry no matter what's in physics
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#ifndef MESOSCALE
       if (moyzon_ch) then ! Zonal averages
          
@@ -716,4 +718,5 @@
 
       else !  if not moyzon
+#endif
       
         DO ig=1,ngrid
@@ -730,5 +733,7 @@
         ENDDO
 
+#ifndef MESOSCALE
       endif  ! moyzon
+#endif
 
       ! -------------------------------------------------------------------------------------
@@ -1028,4 +1033,8 @@
       if (tracer) then
 
+
+#ifndef MESOSCALE
+!! JVO 20 : For now, no chemistry or microphysics in MESOSCALE, but why not in the future ?
+
   ! -------------------
   !   V.1 Microphysics
@@ -1195,4 +1204,6 @@
             
          endif ! end of 'callchim'
+
+#endif
 
   ! ---------------
Index: trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90	(revision 2366)
@@ -1,4 +1,4 @@
 module radcommon_h
-      use radinc_h, only: L_NSPECTI, L_NSPECTV, NTstar, NTstop
+      use radinc_h, only: L_NSPECTI, L_NSPECTV, NTstart, NTstop
       implicit none
 
@@ -78,5 +78,5 @@
       REAL,SAVE :: tstellar ! Stellar brightness temperature (SW)
 
-      real*8,save :: planckir(L_NSPECTI,NTstop-NTstar+1)
+      REAL*8, DIMENSION(:,:), ALLOCATABLE, SAVE :: planckir
 
       real*8,save :: PTOP
Index: trunk/LMDZ.TITAN/libf/phytitan/radinc_h.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/radinc_h.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/radinc_h.F90	(revision 2366)
@@ -72,26 +72,33 @@
 
       ! For Planck function integration: 
-      ! equivalent temperatures are 1/NTfac of these values
-      integer, parameter :: NTstar = 500
-      integer, parameter :: NTstop = 15000 ! new default for all non hot Jupiter runs
-      real*8, parameter :: NTfac = 1.0D+1  
-      !integer, parameter :: NTstar = 1000
-      !integer, parameter :: NTstop = 25000
-      !real*8,parameter :: NTfac = 5.0D+1    
-      !integer, parameter :: NTstar = 2000
-      !integer, parameter :: NTstop = 50000
-      !real*8,parameter :: NTfac = 1.0D+2    
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      ! Integration boundary temperatures are NTstart/NTfac and Ntstop/NTfac
+      ! -- JVO 20 : Now read boundary T and integration dT as inputs in callphys.def
+      !             NTstart, Nstop and NTfac then set by ini_radinc_h
+      !             Smart user can adjust values depending he's running hot or cold atm
+      !             Default is wide range : 30K-1500K, with 0.1K step
+      !             ->  NTstart=300, Nstop=15000, NTfac=10
+      integer :: NTstart, NTstop
+      real*8  :: NTfac
+
 
 contains
 
-  subroutine ini_radinc_h(nbp_lev)
+  subroutine ini_radinc_h(nbp_lev,tplanckmin,tplanckmax,dtplanck)
   ! Initialize module variables
   implicit none
   integer,intent(in) :: nbp_lev
+  real*8, intent(in) :: tplanckmin
+  real*8, intent(in) :: tplanckmax
+  real*8, intent(in) :: dtplanck
   
   L_NLAYRAD = nbp_lev
-  L_LEVELS = 2*(nbp_lev-1)+3
+  L_LEVELS  = 2*(nbp_lev-1)+3
   L_NLEVRAD = nbp_lev+1
-  
+
+  NTfac   = 1.D0 / dtplanck
+  NTstart = int(tplanckmin * NTfac)
+  NTstop  = int(tplanckmax * NTfac)
+ 
   end subroutine
 
Index: trunk/LMDZ.TITAN/libf/phytitan/setspi.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/setspi.F90	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/setspi.F90	(revision 2366)
@@ -22,5 +22,5 @@
 !==================================================================
 
-      use radinc_h,    only: L_NSPECTI,NTstar,NTstop,NTfac
+      use radinc_h,    only: L_NSPECTI,NTstart,NTstop,NTfac
       use radcommon_h, only: BWNI,WNOI,DWNI,WAVEI,planckir,sigma
       use datafile_mod, only: datadir, corrkdir, banddir
@@ -151,5 +151,7 @@
       print*,''
       print*,'setspi: Current Planck integration range:'
-      print*,'T = ',dble(NTstar)/NTfac, ' to ',dble(NTstop)/NTfac,' K.'
+      print*,'T = ',dble(NTstart)/NTfac, ' to ',dble(NTstop)/NTfac,' K.'
+
+      IF(.NOT.ALLOCATED(planckir)) ALLOCATE(planckir(L_NSPECTI,NTstop-NTstart+1))
 
       do NW=1,L_NSPECTI
@@ -158,5 +160,5 @@
          bpa = (b+a)/2.0D0
          bma = (b-a)/2.0D0
-         do nt=NTstar,NTstop
+         do nt=NTstart,NTstop
             T   = dble(NT)/NTfac
             ans = 0.0D0
@@ -167,5 +169,5 @@
             end do
 
-            planckir(NW,nt-NTstar+1) = ans*bma/(PI*DWNI(NW))
+            planckir(NW,nt-NTstart+1) = ans*bma/(PI*DWNI(NW))
          end do
       end do
@@ -174,5 +176,5 @@
       if(forceEC)then
          print*,'setspi: Force F=sigma*eps*T^4 for all values of T!'
-         do nt=NTstar,NTstop
+         do nt=NTstart,NTstop
             plancksum=0.0D0
             T=dble(NT)/NTfac
@@ -180,10 +182,10 @@
             do NW=1,L_NSPECTI
                plancksum=plancksum+  &
-                  planckir(NW,nt-NTstar+1)*DWNI(NW)*pi
+                  planckir(NW,nt-NTstart+1)*DWNI(NW)*pi
             end do
 
             do NW=1,L_NSPECTI
-               planckir(NW,nt-NTstar+1)=     &
-                  planckir(NW,nt-NTstar+1)*  &
+               planckir(NW,nt-NTstart+1)=     &
+                  planckir(NW,nt-NTstart+1)*  &
                           sigma*(dble(nt)/NTfac)**4/plancksum
             end do
@@ -194,7 +196,7 @@
          ! check energy conservation at lower temperature boundary
          plancksum=0.0D0
-         nt=NTstar
+         nt=NTstart
          do NW=1,L_NSPECTI
-            plancksum=plancksum+planckir(NW,nt-NTstar+1)*DWNI(NW)*pi
+            plancksum=plancksum+planckir(NW,nt-NTstart+1)*DWNI(NW)*pi
          end do
          print*,'setspi: At lower limit:'
@@ -206,5 +208,5 @@
          nt=NTstop
          do NW=1,L_NSPECTI
-            plancksum=plancksum+planckir(NW,nt-NTstar+1)*DWNI(NW)*pi
+            plancksum=plancksum+planckir(NW,nt-NTstart+1)*DWNI(NW)*pi
          end do
          print*,'setspi: At upper limit:'
Index: trunk/LMDZ.TITAN/libf/phytitan/sfluxi.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/sfluxi.F	(revision 2365)
+++ trunk/LMDZ.TITAN/libf/phytitan/sfluxi.F	(revision 2366)
@@ -67,6 +67,6 @@
       TSURF = TLEV(L_LEVELS)
 
-      NTS   = int(TSURF*NTfac)-NTstar+1
-      NTT   = int(TTOP *NTfac)-NTstar+1
+      NTS   = int(TSURF*NTfac)-NTstart+1
+      NTT   = int(TTOP *NTfac)-NTstart+1
 
 !JL12 corrects the surface planck function so that its integral is equal to sigma Tsurf^4
