Index: /LMDZ.3.3/trunk/libf/phylmd/1DUTILS.h
===================================================================
--- /LMDZ.3.3/trunk/libf/phylmd/1DUTILS.h	(revision 248)
+++ /LMDZ.3.3/trunk/libf/phylmd/1DUTILS.h	(revision 248)
@@ -0,0 +1,920 @@
+      SUBROUTINE initial0(n,x)
+      IMPLICIT NONE
+      INTEGER n,i
+      REAL x(n)
+      DO 10 i=1,n
+         x(i)=0.
+10    CONTINUE
+      RETURN
+      END
+      SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
+      IMPLICIT NONE
+c=======================================================================
+c   passage d'un champ de la grille scalaire a la grille physique
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c   declarations:
+c   -------------
+
+      INTEGER im,jm,ngrid,nfield
+      REAL pdyn(im,jm,nfield)
+      REAL pfi(ngrid,nfield)
+
+      INTEGER i,j,ifield,ig
+
+c-----------------------------------------------------------------------
+c   calcul:
+c   -------
+
+      DO ifield=1,nfield
+c   traitement des poles
+         DO i=1,im
+            pdyn(i,1,ifield)=pfi(1,ifield)
+            pdyn(i,jm,ifield)=pfi(ngrid,ifield)
+         ENDDO
+
+c   traitement des point normaux
+         DO j=2,jm-1
+	    ig=2+(j-2)*(im-1)
+            CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
+	    pdyn(im,j,ifield)=pdyn(1,j,ifield)
+         ENDDO
+      ENDDO
+
+      RETURN
+      END
+
+      SUBROUTINE disvert1d(llm,kappa,sig,dsig,s,ds,dsig1,sdsig)
+      IMPLICIT NONE
+c
+c=======================================================================
+c
+c
+c    s = sigma ** kappa   :  coordonnee  verticale
+c    dsig(l)            : epaisseur de la couche l ds la coord.  s
+c    sig(l)             : sigma a l'interface des couches l et l-1
+c    ds(l)              : distance entre les couches l et l-1 en coord.s
+c
+c=======================================================================
+c
+c   declarations:
+c   -------------
+c
+      integer llm
+      real kappa,pi,x
+      real sig(llm+1),dsig(llm),s(llm),ds(llm),dsig1(llm),sdsig(llm)
+c
+      integer ll,l,lllm,lllmm1,lllmp1
+      real abid,abid2,som,quoi,quand,snorm,sigbid,sbid
+      REAL alpha,beta,h,zd,dz0,dz1
+      REAL gama,delta,deltaz,np
+
+      real nhaut
+      INTEGER ierr,ierr1,ierr2
+      real puiss
+
+      real asig,bsig,csig,esig,zsig,p,zz,sig1
+      REAL z1,z2
+c
+c-----------------------------------------------------------------------
+c
+      lllm=llm
+      lllmm1=lllm-1
+      lllmp1=lllm+1
+      pi=2.*asin(1.)
+
+      OPEN(99,file='sigma.def',status='old',form='formatted',
+     s   iostat=ierr1)
+      if(ierr1.ne.0) then
+         close(99)
+         open(99,file='esasig.def',status='old',form='formatted',
+     s   iostat=ierr2)
+      endif
+
+
+c-----------------------------------------------------------------------
+c   cas 1 on lit les options dans sigma.def:
+c   ----------------------------------------
+
+
+      if (ierr1.eq.0) then
+         PRINT*,'WARNING Lecture de sigma.def'
+         READ(99,*) deltaz
+         READ(99,*) h
+         READ(99,*) beta
+         READ(99,*) gama
+         READ(99,*) delta
+         READ(99,*) np
+         CLOSE(99)
+         alpha=deltaz/(llm*h)
+          do l= 1, llm
+             dsig(l) = (alpha+(1.-alpha)*exp(-beta*(llm-l)))*
+     $          ( (tanh(gama*l)/tanh(gama*llm))**np +
+     $            (1.-l/FLOAT(llm))*delta )
+          enddo
+
+          sig(1)=1.
+          do l=1,llm-1
+             sig(l+1)=sig(l)*(1.-dsig(l))/(1.+dsig(l))
+          enddo
+          sig(llm+1)=0.
+
+          do l = 1, llm
+             dsig(l) = sig(l)-sig(l+1)
+          enddo
+
+      else if(ierr2.eq.0) then
+
+         PRINT*,'WARNING Lecture de esasig.def'
+         READ(99,*) h
+         READ(99,*) dz0
+         READ(99,*) dz1
+         READ(99,*) nhaut
+         CLOSE(99)
+
+         dz0=dz0/h
+         dz1=dz1/h
+
+         sig1=(1.-dz1)/tanh(.5*(llm-1)/nhaut)
+
+         esig=1.
+
+         PRINT*
+         do l=1,20
+            print*,'esig=',esig
+            esig=-log((1./sig1-1.)*exp(-dz0)/esig)/(llm-1.)
+         enddo
+         PRINT*
+         csig=(1./sig1-1.)/(exp(esig)-1.)
+
+         DO L = 2, llm
+           zz=csig*(exp(esig*(l-1.))-1.)
+           sig(l) =1./(1.+zz)
+     &    * tanh(.5*(llm+1-l)/nhaut)
+         ENDDO
+         sig(1)=1.
+         sig(llm+1)=0.
+
+         do  l = 1, llm
+         dsig(l) =sig(l)-sig(l+1)
+         enddo
+
+      else
+
+         print*,'WARNING!!! Ancienne discretisation verticale'
+c        stop
+         h=7.
+         snorm  = 0.
+         do l = 1, llm
+            x = 2.*asin(1.) * (float(l)-0.5) / float(llm+1)
+            dsig(l) = 1.0 + 7.0 * sin(x)**2
+            snorm = snorm + dsig(l)
+         enddo
+         snorm = 1./snorm
+         do l = 1, llm
+            dsig(l) = dsig(l)*snorm
+         enddo
+         sig(llm+1) = 0.
+         do l = llm, 1, -1
+            sig(l) = sig(l+1) + dsig(l)
+         enddo
+
+      endif
+
+c-----------------------------------------------------------------------
+c   calcul de s, ds, sdsig...
+c   -------------------------
+
+       quoi      = 1. + 2.* kappa
+       s( llm )  = 1.
+       s(lllmm1) = quoi
+       IF( llm.gt.2 )  THEN
+          DO  ll = 2, lllmm1
+             l         = lllmp1 - ll
+             quand     = sig(l+1)/ sig(l)
+             s(l-1)    = quoi * (1.-quand) * s(l)  + quand * s(l+1)
+          ENDDO
+       END IF
+c
+       snorm=(1.-.5*sig(2)+kappa*(1.-sig(2)))*s(1)+.5*sig(2)*s(2)
+       DO l = 1, llm
+          s(l)    = s(l)/ snorm
+       ENDDO
+
+       DO l = 2, llm
+          ds(l)   = s(l-1) - s(l)
+       ENDDO
+       ds(1)  = 1. - s(1)
+c
+       DO l  = 1, llm
+          sdsig(l) = s(l) * dsig(l)
+          dsig1(l)= 1./dsig(l)
+       ENDDO
+
+c-----------------------------------------------------------------------
+c
+c     Diagnostique sur la discretisation verticale:
+c     ---------------------------------------------
+c
+      print*,'Diagnostique de la discretisation verticale'
+      print*
+      print*,'comparaison de sig(l) et (s(l)+s(l+1))/2)**(1/K)'
+      do 14 l=1,llm-1
+         sigbid=(0.5*(s(l)+s(l+1)))**(1./kappa)
+         print*,'sig(',l+1,')  = ',sig(l+1),
+     S           '    valeur approchee :',sigbid,'   ',dsig(l)
+14    continue
+      print*
+      print*,'comparaison de s(l) et (sig(l)+sig(l+1))/2)**K'
+      do 15 l=1,llm
+         sbid=(0.5*(sig(l+1)+sig(l)))**kappa
+         print*,'  s(',l,')  = ',s(l),
+     S           '    valeur approchee :',sbid
+15    continue
+c
+      PRINT*,'Altitude approchee z,dz'
+      PRINT*
+      z1=0.
+      print*,'   l       Z      DZ      Ztop   dsig'
+      DO 18 l=1,llm-1
+         z2=-h*log(sig(l+1))
+         write(*,'(i5,3x,4f8.4)') l,-h*log(s(l))/kappa,z2-z1,z2
+     &    ,dsig(l)
+         write(14,'(3x,i5,1f10.4)') l,-h*log(s(l))/kappa
+         z1=z2
+18    CONTINUE
+      write(*,'(i5,3x,3f8.4)') l,-h*log(s(llm))/kappa
+      write(14,'(3x,i5,1f10.4)') l,-h*log(s(llm))/kappa
+
+
+      RETURN
+      END
+
+c---------------------------------------------------------------
+      subroutine adv_vert2(w,q,plev)
+c   Schéma amont pour l'advection verticale de l'eau
+      implicit none
+
+#include "dimensions.h"
+#include "dimphy.h"
+
+      real w(llm+1),wq(llm+1),q(llm),plev(llm+1),dsig(llm+1)
+
+      integer l
+
+      do l=1,llm
+         wq(l)=q(l)*w(l)
+         dsig(l)=(plev(l)-plev(l+1))/plev(1)
+      enddo
+
+      wq(llm+1)=0.
+      do l=1,llm
+         q(l)=(q(l)*dsig(l)+wq(l+1)-wq(l))/(dsig(l)+w(l+1)-w(l))
+      enddo
+      do l=1,llm
+         wq(l)=q(l)*w(l)
+      enddo
+      wq(llm+1)=0.
+      do l=1,llm
+         q(l)=(q(l)*dsig(l)+wq(l+1)-wq(l))/(dsig(l)+w(l+1)-w(l))
+      enddo
+
+      return
+      end
+c---------------------------------------------------------------
+      subroutine adv_vert(w,q)
+c   Schéma amont pour l'advection verticale de l'eau
+      implicit none
+
+#include "dimensions.h"
+#include "dimphy.h"
+#include "comvert1d.h"
+
+      real w(llm+1),wq(llm+1),q(llm)
+
+      integer l
+
+      do l=1,llm
+         wq(l)=q(l)*w(l)
+      enddo
+
+      wq(llm+1)=0.
+      do l=1,llm
+         q(l)=(q(l)*dsig(l)+wq(l+1)-wq(l))/(dsig(l)+w(l+1)-w(l))
+      enddo
+      do l=1,llm
+         wq(l)=q(l)*w(l)
+      enddo
+      wq(llm+1)=0.
+      do l=1,llm
+         q(l)=(q(l)*dsig(l)+wq(l+1)-wq(l))/(dsig(l)+w(l+1)-w(l))
+      enddo
+
+      return
+      end
+      subroutine ug_vg(lm,sig,ug,vg)
+      implicit none
+
+      integer lm
+
+      real u0,uinf,sigu0,alphau
+      real v0,vinf,sigv0,alphav
+
+      real sig(lm),ug(lm),vg(lm)
+      integer l
+
+      data u0,uinf,sigu0/-9.,6.,0.8/
+c     data v0,vinf,sigv0/-1.,4.,0.98/
+      data v0,vinf,sigv0/-1.,7.,0.98/
+
+
+      alphau=-log((uinf-u0)/uinf)/log(sigu0)
+      alphav=-log((vinf-v0)/vinf)/log(sigv0)
+      print*,'Alpha=',alphau,alphav
+
+      do l=1,lm
+         ug(l)=uinf+(u0-uinf)*sig(l)**alphau
+         vg(l)=vinf+(v0-vinf)*sig(l)**alphav
+      enddo
+
+      return
+      end
+      SUBROUTINE abort_gcm(modname, message, ierr)
+     
+      USE IOIPSL
+C
+C Stops the simulation cleanly, closing files and printing various
+C comments
+C
+C  Input: modname = name of calling program
+C         message = stuff to print
+C         ierr    = severity of situation ( = 0 normal )
+
+      character*20 modname
+      integer ierr
+      character*80 message
+
+      write(*,*) 'in abort_gcm'
+      call histclo
+c     call histclo(2)
+c     call histclo(3)
+c     call histclo(4)
+c     call histclo(5)
+      write(*,*) 'out of histclo'
+      write(*,*) 'Stopping in ', modname
+      write(*,*) 'Reason = ',message
+      if (ierr .eq. 0) then
+        write(*,*) 'Everything is cool'
+      else
+        write(*,*) 'Houston, we have a problem ', ierr
+      endif
+      STOP
+      END
+      FUNCTION q_sat(kelvin, millibar)
+c
+      IMPLICIT none
+c======================================================================
+c Autheur(s): Z.X. Li (LMD/CNRS)
+c Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
+c======================================================================
+c Arguments:
+c kelvin---input-R: temperature en Kelvin
+c millibar--input-R: pression en mb
+c
+c q_sat----output-R: vapeur d'eau saturante en kg/kg
+c======================================================================
+c
+      REAL q_sat, kelvin, millibar
+c
+      REAL r2es
+      PARAMETER (r2es=611.14 *18.0153/28.9644)
+c
+      REAL r3les, r3ies, r3es
+      PARAMETER (R3LES=17.269)
+      PARAMETER (R3IES=21.875)
+c
+      REAL r4les, r4ies, r4es
+      PARAMETER (R4LES=35.86)
+      PARAMETER (R4IES=7.66)
+c
+      REAL rtt
+      PARAMETER (rtt=273.16)
+c
+      REAL retv
+      PARAMETER (retv=28.9644/18.0153 - 1.0)
+c
+      REAL zqsat
+      REAL temp, pres
+C     ------------------------------------------------------------------
+c
+c
+      temp = kelvin
+      pres = millibar * 100.0
+c      write(*,*)'kelvin,millibar=',kelvin,millibar
+c      write(*,*)'temp,pres=',temp,pres
+c
+      IF (temp .LE. rtt) THEN
+         r3es = r3ies
+         r4es = r4ies
+      ELSE
+         r3es = r3les
+         r4es = r4les
+      ENDIF
+c
+      zqsat=r2es/pres * EXP ( r3es*(temp-rtt) / (temp-r4es) )
+      zqsat=MIN(0.5,ZQSAT)
+      zqsat=zqsat/(1.-retv  *zqsat)
+c
+      q_sat = zqsat
+c
+      RETURN
+      END
+      subroutine scopy(n,sx,incx,sy,incy)
+c
+      IMPLICIT NONE
+c
+      integer n,incx,incy,ix,iy,i
+      real sx((n-1)*incx+1),sy((n-1)*incy+1)
+c
+      iy=1
+      ix=1
+      do 10 i=1,n
+	 sy(iy)=sx(ix)
+         ix=ix+incx
+         iy=iy+incy
+10    continue
+c
+      return
+      end
+      subroutine wrgradsfi(if,nl,field,name,titlevar)
+      implicit none
+
+c   Declarations
+
+#include "gradsdef.h"
+
+c   arguments
+      integer if,nl
+      real field(imx*jmx*lmx)
+      character*10 name,file
+      character*10 titlevar
+
+c   local
+
+      integer im,jm,lm,i,j,l,lnblnk,iv,iii,iji,iif,ijf
+
+      logical writectl
+
+
+      writectl=.false.
+
+c     print*,if,iid(if),jid(if),ifd(if),jfd(if)
+      iii=iid(if)
+      iji=jid(if)
+      iif=ifd(if)
+      ijf=jfd(if)
+      im=iif-iii+1
+      jm=ijf-iji+1
+      lm=lmd(if)
+
+c     print*,'im,jm,lm,name,firsttime(if)'
+c     print*,im,jm,lm,name,firsttime(if)
+
+      if(firsttime(if)) then
+         if(name.eq.var(1,if)) then
+            firsttime(if)=.false.
+            ivar(if)=1
+         print*,'fin de l initialiation de l ecriture du fichier'
+         print*,file
+           print*,'fichier no: ',if
+           print*,'unit ',unit(if)
+           print*,'nvar  ',nvar(if)
+           print*,'vars ',(var(iv,if),iv=1,nvar(if))
+         else
+            ivar(if)=ivar(if)+1
+            nvar(if)=ivar(if)
+            var(ivar(if),if)=name
+            tvar(ivar(if),if)=titlevar(1:lnblnk(titlevar))
+            nld(ivar(if),if)=nl
+            print*,'initialisation ecriture de ',var(ivar(if),if)
+            print*,'if ivar(if) nld ',if,ivar(if),nld(ivar(if),if)
+         endif
+         writectl=.true.
+         itime(if)=1
+      else
+         ivar(if)=mod(ivar(if),nvar(if))+1
+         if (ivar(if).eq.nvar(if)) then
+            writectl=.true.
+            itime(if)=itime(if)+1
+         endif
+
+         if(var(ivar(if),if).ne.name) then
+           print*,'Il faut stoker la meme succession de champs a chaque'
+           print*,'pas de temps'
+           print*,'fichier no: ',if
+           print*,'unit ',unit(if)
+           print*,'nvar  ',nvar(if)
+           print*,'vars ',(var(iv,if),iv=1,nvar(if))
+
+           stop
+         endif
+      endif
+
+c     print*,'ivar(if),nvar(if),var(ivar(if),if),writectl'
+c     print*,ivar(if),nvar(if),var(ivar(if),if),writectl
+      do l=1,nl
+         irec(if)=irec(if)+1
+c        print*,'Ecrit rec=',irec(if),iii,iif,iji,ijf,
+c    s (l-1)*imd(if)*jmd(if)+(iji-1)*imd(if)+iii
+c    s ,(l-1)*imd(if)*jmd(if)+(ijf-1)*imd(if)+iif
+         write(unit(if)+1,rec=irec(if))
+     s   ((field((l-1)*imd(if)*jmd(if)+(j-1)*imd(if)+i)
+     s   ,i=iii,iif),j=iji,ijf)
+      enddo
+      if (writectl) then
+
+      file=fichier(if)
+c   WARNING! on reecrase le fichier .ctl a chaque ecriture
+      open(unit(if),file=file(1:lnblnk(file))//'.ctl',
+     &         form='formatted',status='unknown')
+      write(unit(if),'(a5,1x,a40)')
+     &       'DSET ','^'//file(1:lnblnk(file))//'.dat'
+
+      write(unit(if),'(a12)') 'UNDEF 1.0E30'
+      write(unit(if),'(a5,1x,a40)') 'TITLE ',title(if)
+      call formcoord(unit(if),im,xd(iii,if),1.,.false.,'XDEF')
+      call formcoord(unit(if),jm,yd(iji,if),1.,.true.,'YDEF')
+      call formcoord(unit(if),lm,zd(1,if),1.,.false.,'ZDEF')
+      write(unit(if),'(a4,i10,a30)')
+     &       'TDEF ',itime(if),' LINEAR 07AUG1998 30MN '
+      write(unit(if),'(a4,2x,i5)') 'VARS',nvar(if)
+      do iv=1,nvar(if)
+c        print*,'if,var(iv,if),nld(iv,if),nld(iv,if)-1/nld(iv,if)'
+c        print*,if,var(iv,if),nld(iv,if),nld(iv,if)-1/nld(iv,if)
+         write(unit(if),1000) var(iv,if),nld(iv,if)-1/nld(iv,if)
+     &     ,99,tvar(iv,if)
+      enddo
+      write(unit(if),'(a7)') 'ENDVARS'
+c
+1000  format(a5,3x,i4,i3,1x,a39)
+
+      close(unit(if))
+
+      endif ! writectl
+
+      return
+
+      END
+
+      subroutine inigrads(if,im
+     s  ,x,fx,xmin,xmax,jm,y,ymin,ymax,fy,lm,z,fz
+     s  ,dt,file,titlel)
+
+
+      implicit none
+
+      integer if,im,jm,lm,i,j,l,lnblnk
+      real x(im),y(jm),z(lm),fx,fy,fz,dt
+      real xmin,xmax,ymin,ymax
+
+      character file*10,titlel*40
+
+#include "gradsdef.h"
+
+      data unit/24,32,34,36,38,40,42,44,46,48/
+      data nf/0/
+
+      if (if.le.nf) stop'verifier les appels a inigrads'
+
+      print*,'Entree dans inigrads'
+
+      nf=if
+      title(if)=titlel
+      ivar(if)=0
+
+      fichier(if)=file(1:lnblnk(file))
+
+      firsttime(if)=.true.
+      dtime(if)=dt
+
+      iid(if)=1
+      ifd(if)=im
+      imd(if)=im
+      do i=1,im
+         xd(i,if)=x(i)*fx
+         if(xd(i,if).lt.xmin) iid(if)=i+1
+         if(xd(i,if).le.xmax) ifd(if)=i
+      enddo
+      print*,'On stoke du point ',iid(if),'  a ',ifd(if),' en x'
+
+      jid(if)=1
+      jfd(if)=jm
+      jmd(if)=jm
+      do j=1,jm
+         yd(j,if)=y(j)*fy
+         if(yd(j,if).gt.ymax) jid(if)=j+1
+         if(yd(j,if).ge.ymin) jfd(if)=j
+      enddo
+      print*,'On stoke du point ',jid(if),'  a ',jfd(if),' en y'
+
+      print*,'Open de dat'
+      print*,'file=',file
+      print*,'fichier(if)=',fichier(if)
+
+      print*,4*(ifd(if)-iid(if))*(jfd(if)-jid(if))
+      print*,file(1:lnblnk(file))//'.dat'
+
+      OPEN (unit(if)+1,FILE=file(1:lnblnk(file))//'.dat',
+     s   FORM='UNFORMATTED',
+     s   ACCESS='DIRECT'
+     s  ,RECL=4*(ifd(if)-iid(if)+1)*(jfd(if)-jid(if)+1))
+
+      print*,'Open de dat ok'
+
+      lmd(if)=lm
+      do l=1,lm
+         zd(l,if)=z(l)*fz
+      enddo
+
+      irec(if)=0
+
+      print*,if,imd(if),jmd(if),lmd(if)
+      print*,'if,imd(if),jmd(if),lmd(if)'
+
+      return
+      end
+      SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
+      IMPLICIT NONE
+c=======================================================================
+c   passage d'un champ de la grille scalaire a la grille physique
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c   declarations:
+c   -------------
+
+      INTEGER im,jm,ngrid,nfield
+      REAL pdyn(im,jm,nfield)
+      REAL pfi(ngrid,nfield)
+
+      INTEGER j,ifield,ig
+
+c-----------------------------------------------------------------------
+c   calcul:
+c   -------
+
+      IF(ngrid.NE.2+(jm-2)*(im-1).AND.ngrid.NE.1)
+     s    STOP 'probleme de dim'
+c   traitement des poles
+      CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
+      CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
+
+c   traitement des point normaux
+      DO ifield=1,nfield
+         DO j=2,jm-1
+	    ig=2+(j-2)*(im-1)
+            CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
+         ENDDO
+      ENDDO
+
+      RETURN
+      END
+
+      subroutine writelim
+     s   (phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice)
+c
+#include "dimensions.h"
+#include "dimphy.h"
+#include "netcdf.inc"
+
+      REAL phy_nat(klon,360)
+      REAL phy_alb(klon,360)
+      REAL phy_sst(klon,360)
+      REAL phy_bil(klon,360)
+      REAL phy_rug(klon,360)
+      REAL phy_ice(klon,360)
+
+      INTEGER ierr
+      INTEGER dimfirst(3)
+      INTEGER dimlast(3)
+c
+      INTEGER nid, ndim, ntim
+      INTEGER dims(2), debut(2), epais(2)
+      INTEGER id_tim
+      INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB
+
+      PRINT*, 'Ecriture du fichier limit'
+c
+      ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
+c
+      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
+     .                       "Fichier conditions aux limites")
+      ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
+      ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
+c
+      dims(1) = ndim
+      dims(2) = ntim
+c
+ccc      ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
+      ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
+      ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
+     .                        "Jour dans l annee")
+ccc      ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
+      ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
+      ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
+     .                        "Nature du sol (0,1,2,3)")
+ccc      ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
+      ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
+      ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
+     .                        "Temperature superficielle de la mer")
+ccc      ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
+      ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
+      ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
+     .                        "Reference flux de chaleur au sol")
+ccc      ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
+      ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
+      ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
+     .                        "Albedo a la surface")
+ccc      ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
+      ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
+      ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
+     .                        "Rugosite")
+c
+      ierr = NF_ENDDEF(nid)
+c
+      DO k = 1, 360
+c
+      debut(1) = 1
+      debut(2) = k
+      epais(1) = klon
+      epais(2) = 1
+c
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
+      ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais,phy_nat(1,k))
+      ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k))
+      ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k))
+      ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k))
+      ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k))
+#else
+      ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
+      ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais,phy_nat(1,k))
+      ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k))
+      ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k))
+      ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k))
+      ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k))
+#endif
+c
+      ENDDO
+c
+      ierr = NF_CLOSE(nid)
+c
+      return
+      end
+
+      SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
+
+c    Auteur :  P. Le Van .
+c
+      IMPLICIT NONE
+
+#include "dimensions.h"
+#include "paramet.h"
+c
+c=======================================================================
+c
+c
+c    s = sigma ** kappa   :  coordonnee  verticale
+c    dsig(l)            : epaisseur de la couche l ds la coord.  s
+c    sig(l)             : sigma a l'interface des couches l et l-1
+c    ds(l)              : distance entre les couches l et l-1 en coord.s
+c
+c=======================================================================
+c
+      REAL pa,preff
+      REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
+      REAL presnivs(llm)
+c
+c   declarations:
+c   -------------
+c
+      REAL sig(llm+1),dsig(llm)
+c
+      INTEGER l
+      REAL snorm
+      REAL alpha,beta,gama,delta,deltaz,h
+      INTEGER np,ierr
+      REAL pi,x
+
+c-----------------------------------------------------------------------
+c
+      pi=2.*ASIN(1.)
+
+      OPEN(99,file='sigma.def',status='old',form='formatted',
+     s   iostat=ierr)
+
+c-----------------------------------------------------------------------
+c   cas 1 on lit les options dans sigma.def:
+c   ----------------------------------------
+
+      IF (ierr.eq.0) THEN
+
+      print*,'WARNING!!! on lit les options dans sigma.def'
+      READ(99,*) deltaz
+      READ(99,*) h
+      READ(99,*) beta
+      READ(99,*) gama
+      READ(99,*) delta
+      READ(99,*) np
+      CLOSE(99)
+      alpha=deltaz/(llm*h)
+c
+
+       DO 1  l = 1, llm
+       dsig(l) = (alpha+(1.-alpha)*exp(-beta*(llm-l)))*
+     $          ( (tanh(gama*l)/tanh(gama*llm))**np +
+     $            (1.-l/FLOAT(llm))*delta )
+   1   CONTINUE
+
+       sig(1)=1.
+       DO 101 l=1,llm-1
+          sig(l+1)=sig(l)*(1.-dsig(l))/(1.+dsig(l))
+101    CONTINUE
+       sig(llm+1)=0.
+
+       DO 2  l = 1, llm
+       dsig(l) = sig(l)-sig(l+1)
+   2   CONTINUE
+c
+
+      ELSE
+c-----------------------------------------------------------------------
+c   cas 2 ancienne discretisation (LMD5...):
+c   ----------------------------------------
+
+      PRINT*,'WARNING!!! Ancienne discretisation verticale'
+
+      h=7.
+      snorm  = 0.
+      DO l = 1, llm
+         x = 2.*asin(1.) * (FLOAT(l)-0.5) / float(llm+1)
+         dsig(l) = 1.0 + 7.0 * SIN(x)**2
+         snorm = snorm + dsig(l)
+      ENDDO
+      snorm = 1./snorm
+      DO l = 1, llm
+         dsig(l) = dsig(l)*snorm
+      ENDDO
+      sig(llm+1) = 0.
+      DO l = llm, 1, -1
+         sig(l) = sig(l+1) + dsig(l)
+      ENDDO
+
+      ENDIF
+
+
+      DO l=1,llm
+        nivsigs(l) = FLOAT(l)
+      ENDDO
+
+      DO l=1,llmp1
+        nivsig(l)= FLOAT(l)
+      ENDDO
+
+c
+c    ....  Calculs  de ap(l) et de bp(l)  ....
+c    .........................................
+c
+c
+c   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
+c
+
+      bp(llmp1) =   0.
+
+      DO l = 1, llm
+cc
+ccc    ap(l) = 0.
+ccc    bp(l) = sig(l)
+
+      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
+      ap(l) = pa * ( sig(l) - bp(l) )
+c
+      ENDDO
+      ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
+
+      PRINT *,' BP '
+      PRINT *,  bp
+      PRINT *,' AP '
+      PRINT *,  ap
+
+      DO l = 1, llm
+       dpres(l) = bp(l) - bp(l+1)
+       presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
+      ENDDO
+
+      PRINT *,' PRESNIVS '
+      PRINT *,presnivs
+
+      RETURN
+      END
Index: /LMDZ.3.3/trunk/libf/phylmd/1Dconv.h
===================================================================
--- /LMDZ.3.3/trunk/libf/phylmd/1Dconv.h	(revision 248)
+++ /LMDZ.3.3/trunk/libf/phylmd/1Dconv.h	(revision 248)
@@ -0,0 +1,1755 @@
+        SUBROUTINE get_uvd(itap,dtime,tsol,qsol,file_fordat
+     s                    ,ht,hq,hw)
+
+        implicit none
+ 
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c cette routine permet d'obtenir u_convg,v_convg,ht,hq et ainsi de
+c pouvoir calculer la convergence et le cisaillement dans la physiq
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+#include "YOMCST.h"
+
+      INTEGER klev
+      REAL play(100)  !pression en Pa au milieu de chaque couche GCM
+      INTEGER JM(100) !pression en Pa au milieu de chaque couche GCM
+      REAL coef1(100) !coefficient d'interpolation
+      REAL coef2(100) !coefficient d'interpolation
+
+      INTEGER nblvlm !nombre de niveau de pression du mesoNH
+      REAL playm(100)  !pression en Pa au milieu de chaque couche Meso-NH
+      REAL hplaym(100) !pression en hPa milieux des couches Meso-NH
+
+      integer i,j,k,ii,ll,in
+      REAL tsol,qsol
+
+      CHARACTER*80 file_forctl,file_fordat,file_start
+
+      COMMON/physiq1/klev,play,JM,coef1,coef2
+      COMMON/physiq2/nblvlm,playm,hplaym
+
+c======================================================================
+c methode: on va chercher les donnees du mesoNH de meteo france, on y
+c          a acces a tout pas detemps grace a la routine rdgrads qui
+c          est une boucle lisant dans ces fichiers.
+c          Puis on interpole ces donnes sur les 11 niveaux du gcm et
+c          et sur les pas de temps de ce meme gcm
+c======================================================================
+c input:
+c       pasmax     :nombre de pas de temps maximum du mesoNH
+c       dt         :pas de temps du meso_NH (en secondes)
+c----------------------------------------------------------------------
+      integer pasmax,dt
+      save pasmax,dt
+c----------------------------------------------------------------------
+c arguments:
+c           itap   :compteur de la physique(le nombre de ces pas est
+c                   fixe dans la subroutine calcul_ini_gcm de interpo
+c                   -lation
+c           dtime  :pas detemps du gcm (en secondes)
+c           ht     :convergence horizontale de temperature(K/s)
+c           hq     :    "         "       d'humidite (kg/kg/s)
+c           hw     :vitesse verticale moyenne (m/s**2)
+c----------------------------------------------------------------------
+        integer itap
+        real dtime
+        real ht(100)
+        real hq(100)
+        real hw(100)
+c----------------------------------------------------------------------
+c Variables internes de get_uvd (note : l'interpolation temporelle
+c est faite entre les pas de temps before et after, sur les variables
+c definies sur la grille du SCM)
+c     time0     :date initiale en secondes
+c     time      :temps associe a chaque pas
+c     pas       :numero du pas du meso_NH
+c     pasprev   :numero du pas precedent
+c     htaft     :advection horizontale de temp. au pas de temps after
+c     hqaft     :    "         "      d'humidite        "
+c     hwaft     :vitesse verticalle moyenne  au pas de temps after
+c     htbef     :idem htaft, mais pour le pas de temps before
+c     hqbef     :voir hqaft
+c     hwbef     :voir hwaft
+c----------------------------------------------------------------------
+        integer time0,pas,pasprev
+        save time0,pas,pasprev
+        real time
+        real htaft(100),hqaft(100),hwaft(100)
+        save htaft,hqaft,hwaft
+        real htbef(100),hqbef(100),hwbef(100)
+        save htbef,hqbef,hwbef
+        integer timeaft,timebef
+        save timeaft,timebef
+        integer temps
+        character*4 string
+c----------------------------------------------------------------------
+c variables arguments de la subroutine rdgrads
+c---------------------------------------------------------------------
+        integer icompt      !compteur de rdgrads
+        real z(100)         ! altitude (grille Meso)
+        real ht_mes(100)    !convergence horizontale de temperature
+                            !-(grille Meso)
+        real hq_mes(100)    !convergence horizontale d'humidite
+                            !(grille Meso)
+        real hw_mes(100)    !vitesse verticale moyenne
+                            !(grille Meso)
+c
+c---------------------------------------------------------------------
+c variable argument de la subroutine copie
+c---------------------------------------------------------------------
+c SB        real pplay(100)    !pression en milieu de couche du gcm
+c SB                            !argument de la physique
+c---------------------------------------------------------------------
+c variables destinees a la lecture du pas de temps du fichier de donnees
+c---------------------------------------------------------------------
+       character*80 aaa,atemps,spaces,apasmax
+       integer nch,imn,ipa
+c---------------------------------------------------------------------
+c  procedures appelees
+        external rdgrads    !lire en iterant dans forcing.dat
+c---------------------------------------------------------------------
+               print*,'le pas itap est:',itap
+c*** on determine le pas du meso_NH correspondant au nouvel itap ***
+c*** pour aller chercher les champs dans rdgrads                 ***
+        time=time0+itap*dtime
+        temps=int(time/dt+1)
+        pas=min(temps,pasmax)
+             print*,'le pas Meso est:',pas
+c
+c
+c===================================================================
+c
+c*** on remplit les champs before avec les champs after du pas   ***
+c*** precedent en format gcm                                     ***
+        if(pas.gt.pasprev)then
+               do i=1,klev
+                  htbef(i)=htaft(i)
+                  hqbef(i)=hqaft(i)
+                  hwbef(i)=hwaft(i)
+               enddo
+               timebef=pasprev*dt
+               timeaft=timebef+dt
+               icompt=(pas-1)*(nblvlm*4)
+                       print*,'le pas pas est:',pas
+c*** on va chercher les nouveaux champs after dans toga.dat     ***
+c*** champs en format meso_NH                                   ***
+c          open(99,FILE='forcing.dat',FORM='UNFORMATTED',
+
+          write(*,'(a)') 'OPEN dans get_uvd de '//file_fordat
+          open(99,FILE=file_fordat,FORM='UNFORMATTED',
+     .             ACCESS='DIRECT',RECL=4)
+          call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes)
+          do i = 1,nblvlm
+            ht_mes(i) = ht_mes(i)*(hplaym(i)/1000.)**rkappa
+          enddo
+c
+               print*,'ht_mes ',(ht_mes(i),i=1,nblvlm)
+               print*,'hq_mes ',(hq_mes(i),i=1,nblvlm)
+               print*,'hw_mes ',(hw_mes(i),i=1,nblvlm)
+c*** on interpole les champs meso_NH sur les niveaux de pression***
+c*** gcm . on obtient le nouveau champ after                    ***
+            do k=1,klev
+             if (JM(k) .eq. 0) then
+         htaft(k)=coef1(k)*tsol+coef2(k)*ht_mes(jm(k)+1)
+         hqaft(k)=coef1(k)*qsol+coef2(k)*hq_mes(jm(k)+1)
+         hwaft(k)=              coef2(k)*hw_mes(jm(k)+1)
+             else
+         htaft(k)=coef1(k)*ht_mes(jm(k))+coef2(k)*ht_mes(jm(k)+1)
+         hqaft(k)=coef1(k)*hq_mes(jm(k))+coef2(k)*hq_mes(jm(k)+1)
+         hwaft(k)=coef1(k)*hw_mes(jm(k))+coef2(k)*hw_mes(jm(k)+1)
+             endif
+            enddo
+          pasprev=pas
+         else
+                      print*,'timebef est:',timebef
+         endif      !fin du bloc relatif au passage au pas
+                    !de temps (meso) suivant
+c*** si on atteint le pas max des donnees experimentales ,on     ***
+c*** on conserve les derniers champs calcules                    ***
+      if(pas.ge.pasmax)then
+          do ll=1,klev
+               ht(ll)=htaft(ll)
+               hq(ll)=hqaft(ll)
+               hw(ll)=hwaft(ll)
+          enddo
+      else
+c*** on interpole sur les pas de temps de 10mn du gcm a partir   ***
+c** des pas de temps de 1h du meso_NH                            ***
+         do j=1,klev
+         ht(j)=((timeaft-time)*htbef(j)+(time-timebef)*htaft(j))/dt
+         hq(j)=((timeaft-time)*hqbef(j)+(time-timebef)*hqaft(j))/dt
+         hw(j)=((timeaft-time)*hwbef(j)+(time-timebef)*hwaft(j))/dt
+         enddo
+       endif
+c
+c-------------------------------------------------------------------
+c
+         return
+c
+c-----------------------------------------------------------------------
+c on sort les champs de "convergence" pour l'instant initial 'in'
+c ceci se passe au pas temps itap=0 de la physique
+c-----------------------------------------------------------------------
+        entry get_uvd2(itap,file_forctl,file_fordat,file_start
+     s                ,ht,hq,hw)
+             print*,'le pas itap est:',itap
+c
+c===================================================================
+c
+      write(*,*) '   '
+      write(*,*) 'FICHIERS A LIRE DANS GET_UVD2:   '
+      write(*,'(a)') 'fichier forcing.ctl: '//file_forctl
+      write(*,'(a)') 'fichier forcing.dat: '//file_fordat
+      write(*,'(a)') 'fichier start18.data: '//file_start
+      write(*,*) '   '
+
+c!! en attendant de pouvoir compiler les fns CERN, en prescrit
+c!! les variables imn et pasmax a la main...
+c!!
+       write(*,'(a)') 'OPEN '//file_forctl
+       open(97,FILE=file_forctl,FORM='FORMATTED')
+c
+c------------------
+      do i=1,1000
+      read(97,1000,end=999) string
+ 1000 format (a4)
+      if (string .eq. 'TDEF') go to 50
+      enddo
+ 50   backspace(97)
+c-------------------------------------------------------------------
+c   *** on lit le pas de temps dans le fichier de donnees ***
+c   *** "forcing.ctl" et pasmax                           ***
+c-------------------------------------------------------------------
+      read(97,2000) aaa
+2000  format (a80)
+         print*,'aaa est',aaa
+      aaa=spaces(aaa,1)
+         print*,'aaa',aaa
+      call getsch(aaa,' ',' ',5,atemps,nch)
+         print*,'atemps est',atemps
+        atemps=atemps(1:nch-2)
+         print*,'atemps',atemps
+        read(atemps,*) imn
+        dt=imn*60
+         print*,'le pas de temps dt',dt
+      call getsch(aaa,' ',' ',2,apasmax,nch)
+        apasmax=apasmax(1:nch)
+        read(apasmax,*) ipa
+        pasmax=ipa
+         print*,'pasmax est',pasmax
+      CLOSE(97)
+
+c CASE_E:
+c!!         imn = 60
+c!!         ipa = 8
+c TOGA:
+c!!         imn = 360
+c!!         ipa = 480
+
+         dt=imn*60
+         pasmax=ipa
+         print*,'le pas de temps dt',dt
+         print*,'pasmax est',pasmax
+
+
+c------------------------------------------------------------------
+c *** onlit le pas de temps initial de la simulation dans ***
+c *** "start.data"                                        ***
+c------------------------------------------------------------------
+c      open(98,file='start18.data',form='formatted')
+      write(*,'(a)') 'OPEN '//file_start
+      open(98,FILE=file_start,FORM='FORMATTED')
+          read(98,*)in
+                  pasprev=in
+                    print*,'le pas in ini est:',pasprev
+C
+Cjyg     Correction de la date du demarrage.
+CC                  time0=dt*pasprev
+                  time0=dt*(pasprev-1)
+C
+          close(98)
+c
+c      open(99,FILE='forcing.dat',FORM='UNFORMATTED',
+      write(*,'(a)') 'OPEN '//file_fordat
+      open(99,FILE=file_fordat,FORM='UNFORMATTED',
+     .          ACCESS='DIRECT',RECL=4)
+                  icompt=(in-1)*(nblvlm*4)
+        call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes)
+          do i = 1,nblvlm
+            ht_mes(i) = ht_mes(i)*(hplaym(i)/1000.)**rkappa
+          enddo
+c
+               print*,'ht_mes ',(ht_mes(i),i=1,nblvlm)
+               print*,'hq_mes ',(hq_mes(i),i=1,nblvlm)
+               print*,'hw_mes ',(hw_mes(i),i=1,nblvlm)
+c----------------------------------------------------------------------
+c on a obtenu des champs initiaux sur les niveaux du meso_NH
+c on interpole sur les niveaux du gcm(niveau pression bien sur!)
+c-----------------------------------------------------------------------
+         do ii=1,klev
+         htaft(ii)=coef1(ii)*ht_mes(JM(ii))+coef2(ii)*ht_mes(JM(ii)+1)
+         hqaft(ii)=coef1(ii)*hq_mes(JM(ii))+coef2(ii)*hq_mes(JM(ii)+1)
+         hwaft(ii)=coef1(ii)*hw_mes(JM(ii))+coef2(ii)*hw_mes(JM(ii)+1)
+             enddo
+c valeurs initiales des champs de convergence
+          do k=1,klev
+             ht(k)=htaft(k)
+             hq(k)=hqaft(k)
+             hw(k)=hwaft(k)
+          enddo
+        close(99)
+        close(98)
+c
+c-------------------------------------------------------------------
+c
+ 100    return
+c
+999     continue
+        stop 'erreur lecture, file forcing.ctl'
+        end
+ 
+ 
+      SUBROUTINE cool_pool(istep
+     e                    ,n_cooling,dt_cooling,dq_cooling
+     s                    ,dt_cool,dq_cool)
+       implicit none
+C***************************************************************
+C*                                                             *
+C* COOL_POOL                                                   *
+C*                                                             *
+C*                                                             *
+C* written by   : Gilles Foret RAMSES, 15/09/97, 22.00.2       *
+C* modified by :  Sandrine Bony 10/09/98                       *
+C***************************************************************
+c Arguments
+c =========
+c Input
+c -----
+c   istep : numero du pas de temps
+c   n_cooling: nbre de pas de temps ou la pertubation nominale
+c              est appliquee (ensuite, la pertubation decroit
+c              exponentiellement).
+c   dt_cooling : pertubation nominale en temperature
+c   dq_cooling : pertubation nominale en humidite
+c Output
+c ------
+c    dt_cool : pertubation en temperature
+c    dq_cool : pertubation en humidite
+c
+c Variables internes
+c ==================
+c    scale :  facteur applique a la pertubation nominale
+c
+#include "dimensions.h"
+#include "dimphy.h"
+c
+      integer n_cooling,k,istep
+      real dt_cooling(klev),dq_cooling(klev),scale
+      real dt_cool(klev),dq_cool(klev)
+c
+      if (istep .le. n_cooling ) then
+        scale = 1.
+      else
+        scale = 4**(min(15,istep-n_cooling))
+      endif
+c
+        do k = 1,klev
+          dt_cool(k) = dt_cooling(k)/scale
+          dq_cool(k) = dq_cooling(k)/scale
+        enddo
+c
+      return
+      end
+      SUBROUTINE advect_tvl(dtime,t,q,vu_f,vv_f,t_f,q_f
+     :                     ,d_t_adv,d_q_adv)
+      implicit none
+
+#include "dimensions.h"
+#include "dimphy.h"
+
+      integer k
+      real dtime, fact, du, dv, cx, cy, alx, aly
+      real t(klev), q(klev,3)
+     :   , vu_f(klev), vv_f(klev), t_f(klev), q_f(klev,3)
+
+      real d_t_adv(klev), d_q_adv(klev,3)
+
+c Velocity of moving cell
+      data cx,cy /12., -2./
+
+c Dimensions of moving cell
+      data alx,aly /100 000.,150 000./
+
+      do k = 1, klev
+            du = abs(vu_f(k)-cx)/alx
+            dv = abs(vv_f(k)-cy)/aly
+            fact = dtime *(du+dv-du*dv*dtime)
+            d_t_adv(k) = fact * (t_f(k)-t(k))
+            d_q_adv(k,1) = fact * (q_f(k,1)-q(k,1))
+            d_q_adv(k,2) = fact * (q_f(k,2)-q(k,2))
+            d_q_adv(k,3) = fact * (q_f(k,3)-q(k,3))
+      enddo
+
+      return
+      end
+      SUBROUTINE copie(klevgcm,playgcm,psolgcm,file_forctl)
+      implicit none
+
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c cette routine remplit les COMMON physiq1 et physiq2.h
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      INTEGER JM
+      INTEGER klev !nombre de niveau de pression du GCM
+      INTEGER nblvlm !nombre de niveau de pression du mesoNH
+
+      REAL playm(100) !pression en Pa au milieu de chaque couche Meso-NH
+      REAL hplaym(100)!pression en hecto-Pa des milieux de couche Meso-NH
+      REAL play(100)  !pression en Pa au milieu de chaque couche GCM
+      REAL coef1(100)   !coefficient d'interpolation
+      REAL coef2(100)   !coefficient d'interpolation
+
+      COMMON/physiq1/klev,play,JM,coef1,coef2
+      COMMON/physiq2/nblvlm,playm,hplaym
+
+      integer i,k,klevgcm
+      real playgcm(klevgcm) ! pression en milieu de couche du gcm
+      real psolgcm
+      character*80 file_forctl
+
+      klev = klevgcm
+
+c---------------------------------------------------------------------
+c pression au milieu des couches du gcm dans la physiq
+c (SB: remplace le call conv_lipress_gcm(playgcm) )
+c---------------------------------------------------------------------
+
+       do k = 1, klev
+        play(k) = playgcm(k)
+        print*,'la pression gcm est:',play(k)
+       enddo
+
+c----------------------------------------------------------------------
+c lecture du descripteur des donnees Meso-NH (forcing.ctl):
+c  -> nb niveaux du meso.NH (nblvlm) + pressions meso.NH
+c (on remplit le COMMON physiq2)
+c----------------------------------------------------------------------
+
+      call mesolupbis(file_forctl)
+
+      print*,'la valeur de nblvlm est:',nblvlm
+
+c----------------------------------------------------------------------
+c etude de la correspondance entre les niveaux meso.NH et GCM;
+c calcul des coefficients d'interpolation coef1 et coef2
+c (on remplit le COMMON physiq1)
+c----------------------------------------------------------------------
+
+      call corresbis(psolgcm)
+
+c---------------------------------------------------------
+c TEST sur le remplissage de physiq1 et physiq2:
+c---------------------------------------------------------
+ 
+      write(*,*) ' '
+      write(*,*) 'TESTS physiq1 et physiq2 dans copie.F '
+      write(*,*) '--------------------------------------'
+      write(*,*) 'GCM: nb niveaux:',klev,' et pression, coeffs:'
+      do k = 1, klev
+      write(*,*) play(k), coef1(k), coef2(k)
+      enddo
+      write(*,*) 'MESO-NH: nb niveaux:',nblvlm,' et pression:'
+      do k = 1, nblvlm
+      write(*,*) playm(k), hplaym(k)
+      enddo
+      write(*,*) ' '
+
+      end
+      SUBROUTINE writeg1d(ngrid,nx,x,nom,titre)
+      IMPLICIT NONE
+c.......................................................................
+c
+c  ecriture de x pour GRADS-1D
+c
+c  in :
+c         * ngrid      ---> pour controler que l'on est bien en 1D
+c         * nx         ---> taille du vecteur a stocker
+c                             "1" pour une variable de surface
+c                             "nlayer" pour une variable de centre de couche
+c                             "nlayer+1" pour une variable d'interface
+c         * x          ---> variable a stocker
+c         * nom        ---> nom "pour grads"
+c         * titre      ---> titre "pour grads"
+c
+c.......................................................................
+c
+#include "comg1d.h"
+c
+c.......................................................................
+c  declaration des arguments 
+c
+      INTEGER ngrid,nx
+      REAL x(nx)
+      CHARACTER*(*) nom
+      CHARACTER*(*) titre
+c
+c  declaration des arguments 
+c....................................................................... 
+c  declaration des variables locales
+c
+      INTEGER ilayer,ivar
+      LOGICAL test 
+c
+c  declaration des variables locales
+c.......................................................................
+c  contole 1D
+c
+c     print*,'ngrid=',ngrid
+      IF (ngrid.NE.1) return
+c
+c  contole 1D
+c.......................................................................
+c  ouverture du fichier au premier appel
+c
+      IF (g1d_premier) THEN
+        OPEN (g1d_unitfich,FILE=g1d_nomfich
+     &       ,FORM='unformatted',ACCESS='direct',RECL=4)
+        g1d_irec=0
+        g1d_nvar=0
+        g1d_premier=.false.
+      ENDIF
+c
+c  ouverture du fichier au premier appel
+c.......................................................................
+c  pour l'ecriture du fichier ctl
+c
+      test=.true.
+      DO ivar=1,g1d_nvar
+        IF (nom.EQ.g1d_nomvar(ivar)) test=.false.
+      ENDDO
+      IF (test) THEN
+        g1d_nvar=g1d_nvar+1
+        g1d_nomvar(g1d_nvar)=nom
+        g1d_titrevar(g1d_nvar)=titre
+        IF (nx.EQ.1) THEN
+           g1d_dimvar(g1d_nvar)=0
+        ELSEIF (nx.EQ.g1d_nlayer) THEN
+           g1d_dimvar(g1d_nvar)=g1d_nlayer
+        ELSEIF (nx.EQ.g1d_nlayer+1) THEN
+           g1d_dimvar(g1d_nvar)=g1d_nlayer
+        ELSE
+           PRINT *,'._. probleme de dimension dans GRADS-1D ._.'
+        ENDIF
+      ENDIF
+c
+c  pour l'ecriture du fichier ctl
+c.......................................................................
+c  ecriture
+c
+      IF (nx.EQ.1) THEN
+        g1d_irec=g1d_irec+1
+        WRITE(g1d_unitfich,REC=g1d_irec) x(1)
+      ELSE
+        DO ilayer=1,g1d_nlayer
+          g1d_irec=g1d_irec+1
+          WRITE(g1d_unitfich,REC=g1d_irec) x(ilayer)
+        ENDDO
+      ENDIF
+c
+c  ecriture
+c.......................................................................
+c
+10001 CONTINUE
+c
+c.......................................................................
+c
+      RETURN
+      END
+
+
+
+
+
+
+c SB      SUBROUTINE endg1d(ngrid,nlayer,zlayer,ndt)
+      SUBROUTINE endg1d(ngrid,nlayer,player,ndt,dt)
+      IMPLICIT NONE
+c.......................................................................
+c
+c  ecriture du fichier de controle pour GRADS-1D
+c
+c  in :
+c         * ngrid      ---> pour controler que l'on est bien en 1D
+c         * nlayer     ---> nombre de couches
+c         * zlayer     ---> altitude au centre de chaque couche (km)
+c         * player     ---> pression au centre de chaque couche (hPa)
+c         * ndt        ---> nombre de pas de temps
+c         * dt         ---> valeur du pas de temps (s)
+c
+c.......................................................................
+c
+#include "comg1d.h"
+c
+c.......................................................................
+c  declaration des arguments 
+c
+      INTEGER ngrid,nlayer
+c SB      REAL zlayer(nlayer)
+      REAL player(nlayer)
+      INTEGER ndt
+      REAL dt,dtm
+c
+c  declaration des arguments 
+c....................................................................... 
+c  declaration des variables locales
+c
+      INTEGER ivar,ilayer
+c
+c  declaration des variables locales
+c.......................................................................
+c  contole 1D
+c
+      IF (ngrid.NE.1) GOTO 10001
+c
+c  contole 1D
+c.......................................................................
+c
+      IF (nlayer.ne.g1d_nlayer) 
+     &   PRINT *,'._. probleme de dimension dans GRADS-1D ._.'
+c
+c.......................................................................
+c
+      CLOSE (g1d_unitfich)
+c
+c.......................................................................
+c
+      dtm = dt/60.
+
+      OPEN (g1d_unitctl,FILE=g1d_nomctl,FORM='formatted'
+     s     ,status='new')
+      WRITE (g1d_unitctl,'(a4,2x,a20)') 'DSET',g1d_nomfich
+      WRITE (g1d_unitctl,'(a5,2x,a20)') 'UNDEF ','1.E+30'
+      WRITE (g1d_unitctl,'(a11)') 'FORMAT YREV'
+      WRITE (g1d_unitctl,'(a5,2x,a30)') 'TITLE ','champs 1D'
+      WRITE (g1d_unitctl,'(a5,i4,a20)') 'XDEF ',1,' LINEAR 0 1'
+      WRITE (g1d_unitctl,'(a5,i4,a20)') 'YDEF ',1,' LINEAR 0 1'
+      WRITE (g1d_unitctl,'(a5,i4,a20)') 'ZDEF ',g1d_nlayer,' LEVELS'
+      WRITE (g1d_unitctl,'(5(1x,f13.5))')
+c SB     &      (zlayer(ilayer),ilayer=1,g1d_nlayer)
+     &      (player(ilayer)/100.,ilayer=1,g1d_nlayer)
+c SB      WRITE (g1d_unitctl,'(a4,2x,i10,a25)')
+c SB     &      'TDEF ',ndt,' LINEAR 02JAN1987 1HR '
+      WRITE (g1d_unitctl,'(a4,2x,i10,a20,i3,a3)')
+     &      'TDEF ',ndt,' LINEAR 02JAN1987 ',INT(dtm),'MN '
+      WRITE (g1d_unitctl,'(a5,i5)') 'VARS ',g1d_nvar
+      DO ivar=1,g1d_nvar
+      WRITE (g1d_unitctl,'(a5,3x,i4,i3,1x,a39)') 
+     &       g1d_nomvar(ivar),g1d_dimvar(ivar),99,g1d_titrevar(ivar)
+      ENDDO
+      WRITE (g1d_unitctl,'(a7)') 'ENDVARS'
+      CLOSE (g1d_unitctl)
+c
+c.......................................................................
+c
+10001 CONTINUE
+c
+c.......................................................................
+c
+      RETURN
+      END
+      SUBROUTINE mesolupbis(file_forctl)
+      implicit none 
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c Lecture descripteur des donnees MESO-NH (forcing.ctl):
+c -------------------------------------------------------
+c
+c     Cette subroutine lit dans le fichier de controle "essai.ctl"
+c     et affiche le nombre de niveaux du Meso-NH ainsi que les valeurs
+c     des pressions en milieu de couche du Meso-NH (en Pa puis en hPa).
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      INTEGER nblvlm !nombre de niveau de pression du mesoNH
+      REAL playm(100)  !pression en Pa milieu de chaque couche Meso-NH
+      REAL hplaym(100) !pression en hPa des milieux de couche Meso-NH
+      COMMON/physiq2/nblvlm,playm,hplaym
+
+      INTEGER i,lu,k,mlz,mlzh,j
+ 
+      character*80 file_forctl
+
+      character*4 a
+      character*80 aaa,anblvl,spaces
+      integer nch
+
+      lu=9
+c      open (lu,file='forcing.ctl')
+      open(lu,file=file_forctl,form='formatted')
+c
+      do i=1,1000
+      read(lu,1000,end=999) a
+      if (a .eq. 'ZDEF') go to 100
+      enddo
+c
+ 100  backspace(lu)
+      print*,'  DESCRIPTION DES 2 MODELES : '
+      print*,' '
+c
+      read(lu,2000) aaa
+ 2000  format (a80)
+       aaa=spaces(aaa,1)
+       call getsch(aaa,' ',' ',2,anblvl,nch)
+         read(anblvl,*) nblvlm
+
+c      write(*,*) 'ATTENTION! dans mesolupbis on rentre
+c     : nblvlm a la main car pas de bibliotheque CERN..:'
+c CASE_e:
+c!          nblvlm = 43
+c TOGA:
+c!!          nblvlm = 40
+c
+      print*,'nbre de niveaux de pression Meso-NH :',nblvlm
+      print*,' '
+      print*,'pression en Pa de chaque couche du meso-NH :'
+c
+      read(lu,*) (playm(mlz),mlz=1,nblvlm)
+c      Si la pression est en HPa, la multiplier par 100
+      if (playm(1) .lt. 10000.) then
+        do mlz = 1,nblvlm
+         playm(mlz) = playm(mlz)*100.
+        enddo
+      endif
+      print*,(playm(mlz),mlz=1,nblvlm)
+c
+ 1000 format (a4)
+ 1001 format(5x,i2)
+c
+      print*,' '
+      do mlzh=1,nblvlm
+      hplaym(mlzh)=playm(mlzh)/100.
+      enddo
+c
+      print*,'pression en hPa de chaque couche du meso-NH: '
+      print*,(hplaym(mlzh),mlzh=1,nblvlm)
+c
+      close (lu)
+      return
+c
+ 999  stop 'erreur lecture des niveaux pression des donnees'
+      end
+      SUBROUTINE GETSCH(STR,DEL,TRM,NTH,SST,NCH)
+C***************************************************************
+C*                                                             *
+C*                                                             *
+C* GETSCH                                                      *
+C*                                                             *
+C*                                                             *
+C* modified by :                                               *
+C***************************************************************
+C*   Return in SST the character string found between the NTH-1 and NTH
+C*   occurence of the delimiter 'DEL' but before the terminator 'TRM' in
+C*   the input string 'STR'. If TRM=DEL then STR is considered unlimited.
+C*   NCH=Length of the string returned in SST or =-1 if NTH is <1 or if
+C*   NTH is greater than the number of delimiters in STR.
+      IMPLICIT INTEGER (A-Z)
+      CHARACTER STR*(*),DEL*1,TRM*1,SST*(*)
+      NCH=-1
+      SST=' '
+      IF(NTH.GT.0) THEN
+        IF(TRM.EQ.DEL) THEN
+          LENGTH=LEN(STR)
+        ELSE
+          LENGTH=INDEX(STR,TRM)-1
+          IF(LENGTH.LT.0) LENGTH=LEN(STR)
+        ENDIF
+C*     Find beginning and end of the NTH DEL-limited substring in STR
+        END=-1
+        DO 1,N=1,NTH
+        IF(END.EQ.LENGTH) RETURN
+        BEG=END+2
+        END=BEG+INDEX(STR(BEG:LENGTH),DEL)-2
+        IF(END.EQ.BEG-2) END=LENGTH
+C*        PRINT *,'NTH,LENGTH,N,BEG,END=',NTH,LENGTH,N,BEG,END
+    1   CONTINUE
+        NCH=END-BEG+1
+        IF(NCH.GT.0) SST=STR(BEG:END)
+      ENDIF
+      END
+      SUBROUTINE rdgrads(itape,icount,nl,z,ht,hq,hw)
+      IMPLICIT none
+
+      INTEGER itape,icount,icomp, nl
+      real z(nl),ht(nl),hq(nl),hw(nl)
+c
+      INTEGER i, k
+c
+      icomp = icount
+c
+c
+         do k=1,nl
+            icomp=icomp+1
+            read(itape,rec=icomp)z(k)
+         enddo
+         do k=1,nl
+            icomp=icomp+1
+            read(itape,rec=icomp)hT(k)
+         enddo
+         do k=1,nl
+            icomp=icomp+1
+            read(itape,rec=icomp)hQ(k)
+         enddo
+         do k=1,nl
+            icomp=icomp+1
+            read(itape,rec=icomp)hw(k)
+         enddo
+c
+c
+      RETURN
+      END
+      SUBROUTINE corresbis(psol)
+      implicit none
+
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c Cette subroutine calcule et affiche les valeurs des coefficients
+c d'interpolation qui serviront dans la formule d'interpolation elle-
+c meme.
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      INTEGER klev    !nombre de niveau de pression du GCM
+      REAL play(100)  !pression en Pa au milieu de chaque couche GCM
+      INTEGER JM(100)
+      REAL coef1(100) !coefficient d'interpolation
+      REAL coef2(100) !coefficient d'interpolation
+
+      INTEGER nblvlm !nombre de niveau de pression du mesoNH
+      REAL playm(100) !pression en Pa milieu de chaque couche Meso-NH
+      REAL hplaym(100)!pression en hPa des milieux de couche Meso-NH
+
+      COMMON/physiq1/klev,play,JM,coef1,coef2
+      COMMON/physiq2/nblvlm,playm,hplaym
+
+      REAL psol
+      REAL val
+      INTEGER k, mlz, mlzh
+
+
+      do k=1,klev
+         val=play(k)
+       if (val .gt. playm(1)) then
+          mlz = 0
+          JM(1) = mlz
+          coef1(1)=(playm(mlz+1)-val)
+     *             /(playm(mlz+1)-psol)
+          coef2(1)=(val-psol)
+     *             /(playm(mlz+1)-psol)
+       else
+         do mlz=1,nblvlm
+          if (     val .le. playm(mlz)
+     *       .and. val .gt. playm(mlz+1))then
+           JM(k)=mlz
+           coef1(k)=(playm(mlz+1)-val)
+     *              /(playm(mlz+1)-playm(mlz))
+           coef2(k)=(val-playm(mlz))
+     *              /(playm(mlz+1)-playm(mlz))
+          endif
+c
+         enddo
+       endif
+      enddo
+c
+      if (play(klev) .le. playm(nblvlm)) then
+         mlz=nblvlm-1
+         JM(klev)=mlz
+         coef1(klev)=(playm(mlz+1)-val)
+     *            /(playm(mlz+1)-playm(mlz))
+         coef2(klev)=(val-playm(mlz))
+     *            /(playm(mlz+1)-playm(mlz))
+      endif
+c
+      print*,' '
+      print*,'         INTERPOLATION  : '
+      print*,' '
+      print*,'correspondance de 9 niveaux du GCM sur les 53 du meso-NH:'
+      print*,(JM(k),k=1,klev)
+      print*,'correspondance de 9 niveaux du GCM sur les 53 du meso-NH:'
+      print*,(JM(k),k=1,klev)
+      print*,' '
+      print*,'valeurs du premier coef d"interpolation pour les 9 niveaux
+     *: '
+      print*,(coef1(k),k=1,klev)
+      print*,' '
+      print*,'valeurs du deuxieme coef d"interpolation pour les 9 niveau
+     *x: '
+      print*,(coef2(k),k=1,klev)
+c
+      return
+      end
+      SUBROUTINE phyredem (fichnom,dtime,radpas,co2_ppm,solaire,
+     .           rlat,rlon,tsol,tsoil,deltat,qsol,snow,
+     .           radsol,rugmer,agesno,
+     .           zmea,zstd,zsig,zgam,zthe,zpic,zval,rugsrel,
+     .           t_ancien, q_ancien)
+      IMPLICIT none
+c======================================================================
+c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
+c Objet: Ecriture de l'etat de redemarrage pour la physique
+c======================================================================
+#include "dimensions.h"
+#include "dimphy.h"
+#include "netcdf.inc"
+#include "indicesol.h"
+#include "dimsoil.h"
+#include "clesphys.h"
+#include "control.h"
+#include "temps.h"
+c======================================================================
+      CHARACTER*(*) fichnom
+      REAL dtime
+      INTEGER radpas
+      REAL rlat(klon), rlon(klon)
+      REAL co2_ppm
+      REAL solaire
+      REAL tsol(klon,nbsrf)
+      REAL tsoil(klon,nsoilmx,nbsrf)
+      REAL deltat(klon)
+      REAL qsol(klon,nbsrf)
+      REAL snow(klon,nbsrf)
+      REAL radsol(klon)
+      REAL rugmer(klon)
+      REAL agesno(klon)
+      REAL zmea(klon)
+      REAL zstd(klon)
+      REAL zsig(klon)
+      REAL zgam(klon)
+      REAL zthe(klon)
+      REAL zpic(klon)
+      REAL zval(klon)
+      REAL rugsrel(klon)
+      REAL t_ancien(klon,klev), q_ancien(klon,klev)
+c
+      INTEGER nid, nvarid, idim1, idim2, idim3
+      INTEGER ierr
+      INTEGER length
+      PARAMETER (length=100)
+      REAL tab_cntrl(length)
+c
+      INTEGER isoil, nsrf
+      CHARACTER*7 str7
+      CHARACTER*2 str2
+c
+      ierr = NF_CREATE(fichnom, NF_CLOBBER, nid)
+      IF (ierr.NE.NF_NOERR) THEN
+        write(6,*)' Pb d''ouverture du fichier '//fichnom
+        write(6,*)' ierr = ', ierr
+        CALL ABORT
+      ENDIF
+c
+      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 28,
+     .                       "Fichier redemmarage physique")
+c
+      ierr = NF_DEF_DIM (nid, "index", length, idim1)
+      ierr = NF_DEF_DIM (nid, "points_physiques", klon, idim2)
+      ierr = NF_DEF_DIM (nid, "horizon_vertical", klon*klev, idim3)
+c
+      ierr = NF_ENDDEF(nid)
+c
+      DO ierr = 1, length
+         tab_cntrl(ierr) = 0.0
+      ENDDO
+      tab_cntrl(1) = dtime
+      tab_cntrl(2) = radpas
+      tab_cntrl(3) = co2_ppm
+      tab_cntrl(4) = solaire
+      tab_cntrl(5) = iflag_con
+      tab_cntrl(6) = nbapp_rad
+
+      IF( cycle_diurne ) tab_cntrl( 7 ) = 1.
+      IF(   soil_model ) tab_cntrl( 8 ) = 1.
+      IF(     new_oliq ) tab_cntrl( 9 ) = 1.
+      IF(     ok_orodr ) tab_cntrl(10 ) = 1.
+      IF(     ok_orolf ) tab_cntrl(11 ) = 1.
+
+      tab_cntrl(13) = dayref
+      tab_cntrl(14) = anneeref
+      tab_cntrl(13) = day_end
+      tab_cntrl(14) = anne_ini
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "controle", NF_DOUBLE, 1, idim1,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "controle", NF_FLOAT, 1, idim1,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 22,
+     .                        "Parametres de controle")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "longitude", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "longitude", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 32,
+     .                        "Longitudes de la grille physique")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlon)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlon)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "latitude", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "latitude", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 31,
+     .                        "Latitudes de la grille physique")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlat)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlat)
+#endif
+c
+c
+      DO nsrf = 1, nbsrf
+        IF (nsrf.LE.99) THEN
+        WRITE(str2,'(i2.2)') nsrf
+        ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+        ierr = NF_DEF_VAR (nid, "TS"//str2, NF_DOUBLE, 1, idim2,nvarid)
+#else
+        ierr = NF_DEF_VAR (nid, "TS"//str2, NF_FLOAT, 1, idim2,nvarid)
+#endif
+        ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
+     .                        "Temperature de surface No."//str2)
+        ierr = NF_ENDDEF(nid)
+        ELSE
+        PRINT*, "Trop de sous-mailles"
+        CALL abort
+        ENDIF
+#ifdef NC_DOUBLE
+        ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tsol(1,nsrf))
+#else
+        ierr = NF_PUT_VAR_REAL (nid,nvarid,tsol(1,nsrf))
+#endif
+      ENDDO
+c
+      DO nsrf = 1, nbsrf
+      DO isoil=1, nsoilmx
+        IF (isoil.LE.99 .AND. nsrf.LE.99) THEN
+        WRITE(str7,'(i2.2,"srf",i2.2)') isoil,nsrf
+        ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+        ierr = NF_DEF_VAR (nid, "Tsoil"//str7,NF_DOUBLE,1,idim2,nvarid)
+#else
+        ierr = NF_DEF_VAR (nid, "Tsoil"//str7,NF_FLOAT,1,idim2,nvarid)
+#endif
+        ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 29,
+     .                        "Temperature du sol No."//str7)
+        ierr = NF_ENDDEF(nid)
+        ELSE
+        PRINT*, "Trop de couches"
+        CALL abort
+        ENDIF
+#ifdef NC_DOUBLE
+        ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tsoil(1,isoil,nsrf))
+#else
+        ierr = NF_PUT_VAR_REAL (nid,nvarid,tsoil(1,isoil,nsrf))
+#endif
+      ENDDO
+      ENDDO
+c
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "DELTAT", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "DELTAT", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 33,
+     .                        "Ecart de la SST (pour slab-ocean)")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,deltat)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,deltat)
+#endif
+c
+      DO nsrf = 1, nbsrf
+        IF (nsrf.LE.99) THEN
+        WRITE(str2,'(i2.2)') nsrf
+        ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+        ierr = NF_DEF_VAR (nid,"QS"//str2,NF_DOUBLE,1,idim2,nvarid)
+#else
+        ierr = NF_DEF_VAR (nid,"QS"//str2,NF_FLOAT,1,idim2,nvarid)
+#endif
+        ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 25,
+     .                        "Humidite de surface No."//str2)
+        ierr = NF_ENDDEF(nid)
+        ELSE
+        PRINT*, "Trop de sous-mailles"
+        CALL abort
+        ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,qsol(1,nsrf))
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,qsol(1,nsrf))
+#endif
+      ENDDO
+c
+      DO nsrf = 1, nbsrf
+        IF (nsrf.LE.99) THEN
+        WRITE(str2,'(i2.2)') nsrf
+        ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+        ierr = NF_DEF_VAR (nid,"SNOW"//str2,NF_DOUBLE,1,idim2,nvarid)
+#else
+        ierr = NF_DEF_VAR (nid,"SNOW"//str2,NF_FLOAT,1,idim2,nvarid)
+#endif
+        ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 22,
+     .                        "Neige de surface No."//str2)
+        ierr = NF_ENDDEF(nid)
+        ELSE
+        PRINT*, "Trop de sous-mailles"
+        CALL abort
+        ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,snow(1,nsrf))
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,snow(1,nsrf))
+#endif
+      ENDDO
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "RADS", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "RADS", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
+     .                        "Rayonnement net a la surface")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,radsol)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,radsol)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "RUGMER", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "RUGMER", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
+     .                        "Longueur de rugosite sur mer")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rugmer)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rugmer)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "AGESNO", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "AGESNO", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 15,
+     .                        "Age de la neige")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,agesno)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,agesno)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "ZMEA", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "ZMEA", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zmea)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,zmea)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "ZSTD", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "ZSTD", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zstd)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,zstd)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "ZSIG", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "ZSIG", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zsig)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,zsig)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "ZGAM", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "ZGAM", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zgam)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,zgam)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "ZTHE", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "ZTHE", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zthe)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,zthe)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "ZPIC", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "ZPIC", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zpic)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,zpic)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "ZVAL", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "ZVAL", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zval)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,zval)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "RUGSREL", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "RUGSREL", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rugsrel)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rugsrel)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "TANCIEN", NF_DOUBLE, 1, idim3,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "TANCIEN", NF_FLOAT, 1, idim3,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,t_ancien)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,t_ancien)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "QANCIEN", NF_DOUBLE, 1, idim3,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "QANCIEN", NF_FLOAT, 1, idim3,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,q_ancien)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,q_ancien)
+#endif
+c
+      ierr = NF_CLOSE(nid)
+c
+      RETURN
+      END
+      subroutine physdem(lonfi, latfi,phystep,radpas,co2_ppm,
+     .                   solaire, ts, ws, 
+     .                   sn, radsol, deltat, rugmer,
+     .                   agesno, zmea, zstd, zsig,
+     .                   zgam, zthe, zpic, zval,
+     .                   rugsrel)
+
+      IMPLICIT none
+c-------------------------------------------------------------
+C Author : L. Fairhead
+C Date   : 01/10/1999
+C Objet  : Ecriture des etats initiaux physiques
+c-------------------------------------------------------------
+c
+c
+c
+      INTEGER ivap
+      PARAMETER (ivap=1)
+c
+      REAL qsolmax
+      PARAMETER ( qsolmax = 150.0 )
+c
+#include "dimensions.h"
+#include "paramet.h"
+#include "dimphy.h"
+#include "control.h"
+#include "netcdf.inc"
+c
+      INTEGER nid
+
+c Ajout de quelques parametres orographiques (F. LOTT janvier 1995)
+
+      REAL zmea(iip1,jjp1),zstd(iip1,jjp1)
+      REAL zsig(iip1,jjp1),zgam(iip1,jjp1),zthe(iip1,jjp1)
+      REAL zpic(iip1,jjp1),zval(iip1,jjp1)
+      REAL rugsrel(iip1,jjp1)
+      INTEGER idayref,anneeref
+
+
+      integer ierr, idim1, idim2, nvarid
+
+c
+      REAL phystep
+      INTEGER radpas
+      REAL co2_ppm
+      REAL solaire
+      REAL latfi(klon), lonfi(klon)
+      REAL champhys(klon)
+      REAL ts(klon)
+      REAL deltat(klon)
+      REAL ws(klon)
+      REAL sn(klon)
+      REAL radsol(klon)
+      REAL rugmer(klon)
+      REAL agesno(klon)
+      INTEGER length
+      PARAMETER (length=100)
+      REAL tab_cntrl(length)
+      real pi
+
+c
+
+#include "serre.h"
+#include "clesphys.h"
+#include "fxyprim.h"
+c-----------------------------------------------------------------------
+c
+c  stockage sur le fichier Physique:
+c
+      pi=2.*asin(1.)
+      ierr = NF_CREATE("startphy.nc", NF_CLOBBER, nid)
+      IF (ierr.NE.NF_NOERR) THEN
+        WRITE(6,*)' Pb d''ouverture du fichier startphy.nc'
+        WRITE(6,*)' ierr = ', ierr
+        CALL ABORT
+      ENDIF
+c
+      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 28,
+     .                       "Fichier demmarage physique")
+c
+      ierr = NF_DEF_DIM (nid, "index", length, idim1)
+      ierr = NF_DEF_DIM (nid, "points_physiques", klon, idim2)
+c
+      ierr = NF_ENDDEF(nid)
+c
+      DO ierr = 1, length
+         tab_cntrl(ierr) = 0.0
+      ENDDO
+      tab_cntrl(1)  = phystep
+      tab_cntrl(2)  = radpas
+      tab_cntrl(3)  = co2_ppm
+      tab_cntrl(4)  = solaire
+      tab_cntrl(5)  = iflag_con
+      tab_cntrl(6)  = nbapp_rad
+c
+cc     Modif ( P. Le Van )
+c
+       tab_cntrl( 7 ) = 0. 
+       tab_cntrl( 8 ) = 0. 
+       tab_cntrl( 9 ) = 0. 
+       tab_cntrl(10 ) = 0. 
+       tab_cntrl(11 ) = 0. 
+       tab_cntrl(12 ) = 0. 
+
+      IF(  cycle_diurne )  tab_cntrl( 7 ) = 1. 
+      IF(   soil_model  )  tab_cntrl( 8 ) = 1. 
+      IF(    new_oliq   )  tab_cntrl( 9 ) = 1. 
+      IF(    ok_orodr   )  tab_cntrl(10 ) = 1. 
+      IF(    ok_orolf   )  tab_cntrl(11 ) = 1. 
+      IF(  ok_limitvrai )  tab_cntrl(12 ) = 1. 
+
+      tab_cntrl(13)  = dayref
+      tab_cntrl(14)  = anneeref
+
+
+cc   ***    new_oliq   (  commentaires de L. LI dans routine physique ) 
+cc   ***  ok_orodr  et ok_orolf   si on appelle l'orographie      ****
+
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "controle", NF_DOUBLE, 1, idim1,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "controle", NF_FLOAT, 1, idim1,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 22,
+     .                        "Parametres de controle")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "longitude", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "longitude", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 32,
+     .                        "Longitudes de la grille physique")
+      ierr = NF_ENDDEF(nid)
+
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lonfi)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lonfi)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "latitude", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "latitude", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 31,
+     .                        "Latitudes de la grille physique")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,latfi)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,latfi)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "TS", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "TS", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 25,
+     .                        "Temperature de la surface")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ts)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,ts)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "QS", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "QS", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 15,
+     .                        "Humidite du sol")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ws)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,ws)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "SNOW", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "SNOW", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 5,
+     .                        "Neige")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,sn)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,sn)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "RADS", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "RADS", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
+     .                        "Rayonnement net a la surface")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,radsol)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,radsol)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "DELTAT", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "DELTAT", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 33,
+     .                        "Ecart de la SST (pour slab-ocean)")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,deltat)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,deltat)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "RUGMER", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "RUGMER", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
+     .                        "Longueur de rugosite sur mer")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rugmer)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rugmer)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "AGESNO", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "AGESNO", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 15,
+     .                        "Age de la neige")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,agesno)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,agesno)
+#endif
+c
+      CALL gr_dyn_fi(1, iip1, jjp1, klon, zmea, champhys)
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "ZMEA", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "ZMEA", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,champhys)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,champhys)
+#endif
+c
+      CALL gr_dyn_fi(1, iip1, jjp1, klon, zstd, champhys)
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "ZSTD", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "ZSTD", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,champhys)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,champhys)
+#endif
+
+      CALL gr_dyn_fi(1, iip1, jjp1, klon, zsig, champhys)
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "ZSIG", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "ZSIG", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,champhys)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,champhys)
+#endif
+
+      CALL gr_dyn_fi(1, iip1, jjp1, klon, zgam, champhys)
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "ZGAM", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "ZGAM", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,champhys)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,champhys)
+#endif
+
+      CALL gr_dyn_fi(1, iip1, jjp1, klon, zthe, champhys)
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "ZTHE", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "ZTHE", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,champhys)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,champhys)
+#endif
+
+      CALL gr_dyn_fi(1, iip1, jjp1, klon, zpic, champhys)
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "ZPIC", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "ZPIC", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,champhys)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,champhys)
+#endif
+
+      CALL gr_dyn_fi(1, iip1, jjp1, klon, zval, champhys)
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "ZVAL", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "ZVAL", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,champhys)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,champhys)
+#endif
+
+      CALL gr_dyn_fi(1, iip1, jjp1, klon, rugsrel, champhys)
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "RUGSREL", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "RUGSREL", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,champhys)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,champhys)
+#endif
+c
+      ierr = NF_CLOSE(nid)
+
+      RETURN
+
+      END
+*CMZ :          28/02/95  17.58.56  by  Unknown
+*-- Author :
+      CHARACTER*(*) FUNCTION SPACES(STR,NSPACE)
+C
+C CERN PROGLIB# M433    SPACES          .VERSION KERNFOR  4.14  860211
+C ORIG.  6/05/86 M.GOOSSENS/DD
+C
+C-    The function value SPACES returns the character string STR with
+C-    leading blanks removed and each occurence of one or more blanks
+C-    replaced by NSPACE blanks inside the string STR
+C
+      CHARACTER*(*) STR
+C
+      LENSPA = LEN(SPACES)
+      SPACES = ' '
+      IF (NSPACE.LT.0) NSPACE = 0
+      IBLANK = 1
+      ISPACE = 1
+  100 INONBL = INDEXC(STR(IBLANK:),' ')
+      IF (INONBL.EQ.0) THEN
+          SPACES(ISPACE:) = STR(IBLANK:)
+                                                    GO TO 999
+      ENDIF
+      INONBL = INONBL + IBLANK - 1
+      IBLANK = INDEX(STR(INONBL:),' ')
+      IF (IBLANK.EQ.0) THEN
+          SPACES(ISPACE:) = STR(INONBL:)
+                                                    GO TO 999
+      ENDIF
+      IBLANK = IBLANK + INONBL - 1
+      SPACES(ISPACE:) = STR(INONBL:IBLANK-1)
+      ISPACE = ISPACE + IBLANK - INONBL + NSPACE
+      IF (ISPACE.LE.LENSPA)                         GO TO 100
+  999 END
+      FUNCTION INDEXC(STR,SSTR)
+C
+C CERN PROGLIB# M433    INDEXC          .VERSION KERNFOR  4.14  860211
+C ORIG. 26/03/86 M.GOOSSENS/DD
+C
+C-    Find the leftmost position where substring SSTR does not match
+C-    string STR scanning forward
+C
+      CHARACTER*(*) STR,SSTR
+C
+      LENS   = LEN(STR)
+      LENSS  = LEN(SSTR)
+C
+      DO 10 I=1,LENS-LENSS+1
+          IF (STR(I:I+LENSS-1).NE.SSTR) THEN
+              INDEXC = I
+                                         GO TO 999
+          ENDIF
+   10 CONTINUE
+      INDEXC = 0
+C
+  999 END
Index: /LMDZ.3.3/trunk/libf/phylmd/comg1d.h
===================================================================
--- /LMDZ.3.3/trunk/libf/phylmd/comg1d.h	(revision 248)
+++ /LMDZ.3.3/trunk/libf/phylmd/comg1d.h	(revision 248)
@@ -0,0 +1,52 @@
+c.......................................................................
+c  le COMMON pour GRADS-1D
+c
+c  on peut se dire : "on ne sauvera pas plus de 1000 variables ... hein ?"
+c
+      INTEGER g1d_nvarmx
+      PARAMETER(g1d_nvarmx=1000)
+c
+c         * g1d_nlayer     ---> nombre de couches verticales
+c         * g1d_nomfich    ---> nom du fichier grads
+c         * g1d_unitfich   ---> code du fichier grads
+c         * g1d_nomctl     ---> nom du fichier ctl
+c         * g1d_unitctl    ---> code du fichier ctl
+c         * g1d_premier    ---> variable logique pour dire si le fichier
+c                               est deja ouvert
+c         * g1d_irec       ---> indice de derniere ecriture
+c         * g1d_nvar       ---> nombre de variables deja definies a la
+c                               derniere ecriture
+c         * g1d_nomvar     ---> noms des vecteurs existants
+c         * g1d_dimvar     ---> taille des vecteurs
+c         * g1d_titrevar   ---> titres des vecteurs
+c         * g1d_tmp1       ---> caractere 
+c         * g1d_tmp2       ---> caractere 
+c
+      INTEGER g1d_nlayer
+      CHARACTER*100 g1d_nomfich
+      INTEGER g1d_unitfich
+      CHARACTER*100 g1d_nomctl
+      INTEGER g1d_unitctl
+      LOGICAL g1d_premier
+      INTEGER g1d_irec
+      INTEGER g1d_nvar
+      CHARACTER*100 g1d_nomvar
+      INTEGER g1d_dimvar
+      CHARACTER*100 g1d_titrevar
+      CHARACTER*100 g1d_tmp1,g1d_tmp2
+c
+      COMMON/COMG1DI/g1d_nlayer
+     &             ,g1d_unitfich
+     &             ,g1d_unitctl
+     &             ,g1d_irec
+     &             ,g1d_nvar
+      COMMON/COMG1DC/g1d_dimvar(0:g1d_nvarmx)
+     &             ,g1d_nomfich
+     &             ,g1d_nomctl
+     &             ,g1d_nomvar(0:g1d_nvarmx)
+     &             ,g1d_titrevar(0:g1d_nvarmx)
+     &             ,g1d_tmp1
+     &             ,g1d_tmp2
+      COMMON/COMG1DL/g1d_premier
+c
+c.......................................................................
Index: /LMDZ.3.3/trunk/libf/phylmd/comvert1d.h
===================================================================
--- /LMDZ.3.3/trunk/libf/phylmd/comvert1d.h	(revision 248)
+++ /LMDZ.3.3/trunk/libf/phylmd/comvert1d.h	(revision 248)
@@ -0,0 +1,10 @@
+c-----------------------------------------------------------------------
+c   INCLUDE 'comvert.h'
+
+      COMMON/comvert/
+     . s(klev),sig(klev+1),ds(klev),dsig(klev),dsig1(klev),sdsig(klev)
+     . ,sig_s(klev)
+
+      REAL s,sig,ds,dsig,dsig1,sdsig,sig_s
+
+c-----------------------------------------------------------------------
Index: /LMDZ.3.3/trunk/libf/phylmd/testphys1d.F
===================================================================
--- /LMDZ.3.3/trunk/libf/phylmd/testphys1d.F	(revision 248)
+++ /LMDZ.3.3/trunk/libf/phylmd/testphys1d.F	(revision 248)
@@ -0,0 +1,673 @@
+      PROGRAM gcm1d
+      IMPLICIT NONE
+
+#include "dimensions.h"
+#include "dimphy.h"
+#include "YOMCST.h"
+#include "comg1d.h"
+#include "clesphys.h"
+#include "control.h"
+
+c    Arguments :
+c    -----------
+
+      integer ngrid,nlayer,longcles,nqmax
+      parameter (longcles=20,nqmax=3)
+cfleur on  ne se sert pas de nqmax mais de nqmx defini dnas dimensions.h depend de la 
+c      compilation
+      integer radpas
+      integer nday,it,iflag_con,unit,i,itap,n_cooling
+      integer iphys_ver,iadv_tvl,i_cvg,i_hum
+
+      real fnday
+      real h,kappa,z1,z2,sbid,sigbid
+      real plev(klev+1),play(klev), psol
+      real temp(klev),u(klev),v(klev),tsurf,q(klev,nqmx),w(klev+1)
+      real ug(klev),vg(klev)
+      real du(klev),dv(klev),dt(klev),dpsrf,dq(klev,nqmx)
+      real du_dyn(klev),dv_dyn(klev)
+     :    ,dt_dyn(klev),dq_dyn(klev,nqmax)
+      real phis,presnivs(klev),clesphy0(longcles)
+      real time,timestep,ecritphy,day,tho
+      real co2_ppm,solaire
+      real rlat,rlon,tsol,radsol,psol_f,tsol_f,qsol_f,qsol
+      real rugmer,rugsrel,snow,agesno,deltat,zmea,zstd,zsig,sn
+      real zgam,zthe,zpic,zval
+      real ema_sig(klev),ema_w(klev),ncst_cbmf,ema_cbmfz,ema_pcb
+     : ,zz_f(klev),vu_f(klev),vv_f(klev),t_f(klev),q_f(klev,nqmax)
+
+      real temp0(klev),q0(klev,3)
+
+      real dt_cooling(klev),dq_cooling(klev)
+      real d_t_cool(klev),d_q_cool(klev)
+      real d_t_adv(klev),d_q_adv(klev,nqmax)
+      real d_t_cvg(klev),d_q_cvg(klev)
+      real ht(100),hq(100),hw(100)
+
+      real phy_nat(360)
+      real phy_alb(360)
+      real phy_sst(360)
+      real phy_bil(360)
+      real phy_rug(360)
+      real phy_ice(360)
+
+      logical cycle_diurne,soil_model,new_oliq,ok_orodr
+     :       ,ok_orolf,ok_limitvrai
+
+      logical firstcall,lastcall,flag_cool
+      logical itsourcecont
+c      itsourcecont permet de choisir entre les sources pour plus de 5 traceurs
+      character*80 ans,file_forctl, file_fordat, file_start,file
+      character*2 cnbl
+
+c-----------------------------------------------------------------------
+
+      COMMON/comvert/
+     * s(llm),sig(llm+1),ds(llm),dsig(llm),dsig1(llm),sdsig(llm)
+     . ,sig_s(llm)
+
+      REAL s,sig,ds,dsig,dsig1,sdsig,sig_s
+
+c-----------------------------------------------------------------------
+c INCLUDE 'temps.h':
+
+      COMMON/temps/itaufin,dtd,
+     s  day_ini,day_end,anne_ini
+
+      INTEGER  itaufin
+      INTEGER*4 day_ini,day_end,anne_ini
+      REAL     dtd
+
+c-----------------------------------------------------------------------
+c   dynamical tendencies
+
+      INTEGER l,ierr,aslun,nlevel,iq,ll
+
+      REAL longitude,latitude
+      REAL zlay(klev),phi(klev)
+      REAL paire 
+
+      DATA latitude,longitude/0.,0./
+
+c-----------------------------------------------------------------------
+c    Initialisations  des constantes
+c    -------------------------------
+
+c constantes
+c ----------
+
+      time=0.
+      tho=3600. 
+      it=0
+
+      call suphec
+
+c parametres lus dans execution_1D:
+c ---------------------------------
+
+      read(*,*) fnday
+      print*, 'fnday',fnday
+      read(*,*) ecritphy
+      print*, 'ecritphy',ecritphy
+      read(*,*) timestep
+      print*, 'timestep',timestep
+      read(*,*) iflag_con 
+      print*, 'iflag_con',iflag_con
+      read(*,*) flag_cool 
+      print*, 'flag_cool',flag_cool
+      read(*,'(a)') cnbl
+      print*, 'cnbl',cnbl
+      read(*,*) radpas   
+      print*, 'radpas',radpas
+      read(*,*) anne_ini ! sans importance mais il faut qqchose
+      print*, 'anne_ini',anne_ini
+      
+c si ecritphy=0: on ecrit la physique a chaque pas de temps:
+
+      if (ecritphy.eq.0.) then
+       ecritphy = timestep/rday
+      endif
+      write(*,*) 'ECRITPHY = ',ecritphy
+
+      firstcall=.true.
+      lastcall=.false.
+      day=1.
+      itsourcecont=.true.
+
+      open(78,file='transil.dat',form='unformatted',
+     s    access='direct',recl=4*llm*llm)
+
+      ngrid=1
+      IF (ngrid.NE.klon) THEN
+         PRINT*,'STOP in inifis'
+         PRINT*,'Probleme de dimenesions :'
+         PRINT*,'ngrid     = ',ngrid
+         PRINT*,'klon   = ',klon
+         STOP
+      ENDIF
+      nlayer=klev
+      nlevel=nlayer+1
+
+c   -------------------------------------------------------------------
+c   Initialisation de la discretisation verticale
+c   ---------------------------------------------
+c sb      CALL disvert(llm,rkappa,sig,dsig,s,ds,dsig1,sdsig)
+
+c   -------------------------------------------------------------------
+c   Profils initiaux.
+c   -----------------
+
+      unit = 80
+      open(unit,file='start'//cnbl//'.data',form='formatted')
+      read(unit,*) 
+      read(unit,*) 
+      read(unit,*) (play(l),l=1,nlayer)
+      read(unit,*) (plev(l),l=1,nlevel)
+      read(unit,*) (temp(l),l=1,nlayer)
+      read(unit,*) (q(l,1),l=1,nlayer)
+      read(unit,*) (q(l,2),l=1,nlayer)
+      read(unit,*) (ema_sig(l),l=1,nlayer)
+      read(unit,*) (ema_w(l),l=1,nlayer) ! reinitialise + loin???
+      read(unit,*) ncst_cbmf, ema_cbmfz, ema_pcb
+      read(unit,*) (zz_f(l),l=1,nlayer)
+      read(unit,*) (vu_f(l),l=1,nlayer)
+      read(unit,*) (vv_f(l),l=1,nlayer)
+      close(unit)
+      write(*,*) 'Lecture fichier start.data ok'
+
+      DO l = 1, nlayer
+         u(l) = vu_f(l)
+         v(l) = vv_f(l)
+c##################
+cfleur ATTENTION q(3) ce nest plus la glace mais le radon si on active phytrac
+c##################3
+         q(l,3) = 0. ! on initialise la glace a zero
+         PRINT*,plev(l),play(l)
+c on recopie le profil initial dans les champs de forcing:
+         t_f(l) = temp(l)
+         q_f(l,1) = q(l,1)
+         q_f(l,2) = q(l,2)
+         q_f(l,3) = q(l,3)
+      ENDDO
+
+
+c Lecture/creation des conditions aux limites:
+c --------------------------------------------
+
+      unit = 81
+      open(unit,file='startphy.data',form='formatted')
+      read(unit,*)
+      read(unit,*) co2_ppm
+      read(unit,*) solaire
+      read(unit,*) rlat
+      read(unit,*) tsol
+      read(unit,*) radsol
+      close(unit)
+      write(*,*) 'Lecture fichier startphy.data ok'
+
+      unit = 82
+      open(unit,file='condsurf.data',form='formatted')
+      read(unit,*) psol_f
+      read(unit,*) tsol_f
+      read(unit,*) qsol_f
+      close(unit)
+      write(*,*) 'Lecture fichier condsurf.data ok'
+
+      if (play(1).lt.10000.) then
+        do l = 1, nlayer
+         play(l) = play(l)*100.
+        enddo
+      endif
+      if (plev(1).lt.10000.) then
+        do l = 1, nlevel
+         plev(l) = plev(l)*100.
+        enddo
+      endif
+
+      if (abs(psol_f-plev(1)) .gt. 1.) then
+         print *,' Incompatibilite entre psol et profil'
+     :   ,' de pression'
+         print *,' psol = ',psol_f,' plev(1) = ',plev(1)
+         stop
+      endif
+
+      tsol    = tsol_f ! temp au sol prise dans condsurf.data
+      psol    = psol_f ! pression au sol prise dans condsurf.data
+      qsol    = qsol_f
+      rugmer  = 0.0001 ! valeur de cchlim.data
+      rugsrel = 0.0    ! (rugsrel = rugoro)
+      snow    = 0.0
+      sn=0.
+      agesno  = 50.0
+      rlon    = 0.0
+      deltat  = 0.0    ! ne sert que pour les slab_ocean
+      phis = 0.
+      zmea = 0.
+      zstd = 0.
+      zsig = 0.
+      zgam = 0.
+      zthe = 0.
+      zpic = 0.
+      zval = 0.
+
+      do i=1,360
+      phy_sst(i) = tsol
+      phy_nat(i) = 0.0  ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
+      phy_alb(i) = 0.15 ! albedo land only (old value condsurf_jyg=0.3)
+      phy_bil(i) = 1.0  ! ne sert que pour les slab_ocean
+      phy_rug(i) = 0.1  ! longueur rugosite utilisee sur land only 
+      phy_ice(i) = 0.0  ! fraction de glace (?)
+      enddo
+
+      call writelim
+     s   (phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice)
+
+
+
+c controles du run (en 3D: lus dans run.def) :
+
+      cycle_diurne = .FALSE.
+      soil_model   = .FALSE.
+      new_oliq     = .FALSE.
+      ok_orodr     = .FALSE.
+      ok_orolf     = .FALSE.
+      ok_limitvrai = .FALSE.
+
+      do i = 1, longcles
+       clesphy0(i) = 0.
+      enddo
+
+      nbapp_rad = NINT(86400./radpas/timestep)
+
+      clesphy0(1) = FLOAT( iflag_con )
+      clesphy0(2) = FLOAT( nbapp_rad )
+      if (cycle_diurne ) clesphy0(3) = 1.
+      if (  soil_model ) clesphy0(4) = 1.
+      if (    new_oliq ) clesphy0(5) = 1.
+      if (    ok_orodr ) clesphy0(6) = 1.
+      if (    ok_orolf ) clesphy0(7) = 1.
+      if ( ok_limitvrai) clesphy0(8) = 1.
+
+
+c -------------------------------------------------------------------
+c Discretisation verticale:
+c -------------------------
+
+c calcul sig,dsig,s,ds,dsig1,sdsig a partir des play, plev lus:
+
+      do l = 1, nlevel
+       sig(l)= plev(l)/plev(1)
+      enddo
+
+      do l = 1, nlayer
+       s(l)= sig(l)**RKAPPA
+      enddo
+
+      do l = 2, nlayer
+       ds(l) = s(l-1) - s(l)
+      enddo
+      ds(1)  = 1. - s(1)
+
+      do l = 1, nlayer
+       dsig(l)  = sig(l)-sig(l+1)
+       sdsig(l) = s(l) * dsig(l)
+       dsig1(l) = 1./dsig(l)
+       presnivs(l) = play(l)/100./100.
+      enddo
+
+      print*,'Diagnostique de la discretisation verticale'
+
+      print*
+      h=7.
+      kappa = RKAPPA
+      print*,'comparaison de sig(l) et (s(l)+s(l+1))/2)**(1/K)'
+      do 14 l=1,llm-1
+         sigbid=(0.5*(s(l)+s(l+1)))**(1./kappa)
+         print*,'sig(',l+1,')  = ',sig(l+1),
+     :           '    valeur approchee :',sigbid,'   ',dsig(l)
+14    continue
+      print*
+      print*,'comparaison de s(l) et (sig(l)+sig(l+1))/2)**K'
+      do 15 l=1,llm
+         sbid=(0.5*(sig(l+1)+sig(l)))**kappa
+         print*,'  s(',l,')  = ',s(l),
+     :           '    valeur approchee :',sbid
+15    continue
+
+      print*,'Altitude approchee z,dz'
+      print*
+      z1=0.
+      print*,'   l       Z      DZ      Ztop   dsig'
+      DO 18 l=1,llm-1
+         z2=-h*log(sig(l+1))
+         write(*,'(i5,3x,4f8.4)') l,-h*log(s(l))/kappa,z2-z1,z2
+     :    ,dsig(l)
+         write(14,'(3x,i5,1f10.4)') l,-h*log(s(l))/kappa
+         z1=z2
+18    CONTINUE
+      write(*,'(i5,3x,3f8.4)') l,-h*log(s(llm))/kappa
+      write(14,'(3x,i5,1f10.4)') l,-h*log(s(llm))/kappa
+
+c      print*,'Correspondance approx pression-altitude:'
+c      z1 = 0.
+c      do l = 1, llm
+c       z1 = z1 + 287.*temp(l)/9.81/play(l)*(plev(l)-plev(l+1))
+c       write(*,*) l,play(l),z1
+c      enddo
+
+c -------------------------------------------------------------------
+c Ecriture de l'etat initial:
+c ---------------------------
+
+c     call ini_fis(timestep,radpas,iflag_con,anne_ini
+c    :            ,co2_ppm,solaire,rlat,rlon,tsol,deltat
+c    :            ,qsol,snow,radsol,rugmer,agesno,zmea,zstd,zsig,zgam
+c    :            ,zthe,zpic,zval,rugsrel
+c    :            ,phy_sst,phy_nat,phy_alb,phy_bil,phy_rug,phy_ice)
+
+      call physdem(rlon, rlat, timestep,radpas,co2_ppm,
+     .                   solaire,tsol, qsol,
+     .                   sn, radsol, deltat, rugmer,
+     .                   agesno, zmea, zstd, zsig,
+     .                   zgam, zthe, zpic, zval,
+     .                   rugsrel)
+
+
+
+c -------------------------------------------------------------------
+c Options de la simulation:
+c -------------------------
+
+      unit = 83
+      open(unit,file='version.data',form='formatted')
+      read(unit,*) iphys_ver,iadv_tvl,i_cvg,i_hum
+      close(unit)
+      write(*,*) 'Lecture fichier version.data ok'
+
+      write(*,*) 'physiq # ',iphys_ver,' (0=std; 1=forced glob cis;
+     : 2=forced loc cis)'
+      write(*,*) 'advection switch : ',iadv_tvl,' (0=no adv ;
+     : 1=forced adv)'
+      write(*,*) 'convergence switch : ',i_cvg,' (0=no conv;
+     : 1=forced conv)'
+
+
+
+c Eventuellement, preparation de la "bulle froide":
+c -------------------------------------------------
+
+      IF (flag_cool) THEN
+      unit = 84
+      open(unit,file='cool_buble.data',form='formatted')
+      read(unit,*) n_cooling
+      do l = 1, nlayer
+      read(unit,*) dt_cooling(l),dq_cooling(l)
+      enddo
+      close(unit)
+      write(*,*) 'Lecture fichier cool_buble.data ok'
+      ENDIF
+
+
+c Eventuellement, preparation du forcage par la convergence:
+c ----------------------------------------------------------
+
+      IF (i_cvg .EQ. 1) THEN
+      file_forctl = 'forcing.ctl'
+      file_fordat = 'forcing.dat'
+      file_start  = 'start'//cnbl//'.data'
+
+      call copie(klev,play,psol,file_forctl)
+      call get_uvd2(itap,file_forctl,file_fordat,file_start
+     :             ,ht,hq,hw)
+      ENDIF
+
+
+c-----------------------------------------------------------------------
+c  initialisation pour GRADS-1D
+c  ----------------------------
+
+      g1d_nlayer=nlayer
+      g1d_nomfich='grads1d.dat'
+      g1d_unitfich=30
+      g1d_nomctl='grads1d.ctl'
+      g1d_unitctl=31
+      g1d_premier=.true.
+
+      file='sort'
+      call inigrads(1,1
+     s  ,0.,1.,-2.,2.,1,0.,-2.,2.,1.
+     s  ,llm,presnivs,1000.
+     s  ,1800.,file,'Diagconvect')
+      print*,'Fin de Initialisation de wrgras'
+
+
+c=======================================================================
+c   DEBUT DE L'INTEGRATION TEMPORELLE:
+c   ==================================
+
+      itap = 1
+
+1     continue
+
+c  calcul du geopotentiel:
+c  -----------------------
+
+      phi(1)=RD*temp(1)*
+     :(plev(1)-play(1))/(.5*(plev(1)+play(1)))
+      do l = 1, nlayer-1
+         phi(l+1)=phi(l)+RD*(temp(l)+temp(l+1))*
+     :   (play(l)-play(l+1))/(play(l)+play(l+1))
+      enddo
+
+c        PRINT*,'altitude km  et T (K)'
+c        WRITE(6,'(2f10.2)') (.001*phi(l)/g,temp(l),l=1,nlayer)
+
+
+c pour etre coherent avec LMDZ.3, on passe 2 arguments 
+c en plus a physiq: phis et aire (pas utilises en 1D):
+      paire = 1. ! aire de la maille
+
+
+c  eventuellement, induction de la convection par chauffage du
+c  bas et refroidissement du haut de la couche limite:
+c  ------------------------------------------------------------
+
+      IF (flag_cool) THEN
+
+      call cool_pool(itap
+     e           ,n_cooling,dt_cooling,dq_cooling
+     s           ,d_t_cool,d_q_cool)
+
+      do l = 1, nlayer
+       temp(l) = temp(l) + d_t_cool(l)
+       q(l,1)  = q(l,1)  + d_q_cool(l)
+      enddo
+
+      ELSE
+
+      do l = 1, nlayer 
+       d_t_cool(l) = 0.
+       d_q_cool(l) = 0.
+      enddo
+
+      ENDIF
+
+
+c  eventuellement, ajouter une tendance relative a une
+c  translation du domaine (ex: pour suivre une ligne de grains):
+c  -------------------------------------------------------------
+
+      IF (iadv_tvl .EQ. 1) THEN
+
+      call advect_tvl(timestep,temp,q,vu_f,vv_f,t_f,q_f
+     :               ,d_t_adv,d_q_adv)
+
+      do l = 1, nlayer
+       temp(l) = temp(l) + d_t_adv(l)
+       q(l,1)  = q(l,1)  + d_q_adv(l,1)
+       q(l,2)  = q(l,2)  + d_q_adv(l,2)
+       q(l,3)  = q(l,3)  + d_q_adv(l,3)
+      enddo
+
+      ELSE
+
+      do l = 1, nlayer 
+       d_t_adv(l) = 0.
+       d_q_adv(l,1) = 0.
+       d_q_adv(l,2) = 0.
+       d_q_adv(l,3) = 0.
+      enddo
+
+      ENDIF
+
+
+c  eventuellement, ajouter une tendance relative a la
+c  convergence de grande echelle:
+c  -------------------------------------------------------------
+
+      IF (i_cvg .EQ. 1) THEN
+
+      call get_uvd(itap,timestep,tsol,qsol,file_fordat,ht,hq,hw)
+
+      do l = 1, nlayer
+
+       temp0(l) = temp(l) ! memoire de "l'avant dynamique"
+       q0(l,1) = q(l,1)
+       q0(l,2) = q(l,2)
+       q0(l,3) = q(l,3)
+
+       d_t_cvg(l) = ht(l) * timestep
+       d_q_cvg(l) = hq(l) * timestep
+       temp(l) = temp(l) + d_t_cvg(l)
+       q(l,1)  = q(l,1)  + d_q_cvg(l)
+c       write(*,*)'d_t_cvg,d_q_cvg:',d_t_cvg(l),d_q_cvg(l) 
+       if (q(l,1).lt.0.) then ! sb
+       print*,'OVAP negative dans main!'
+       print*,'itap,l,q,d_q_cvg:',itap,l,q(l,1),d_q_cvg(l),hq(l)
+       q(l,1)  = MAX(q(l,1),1.e-10) ! evite les humidites negatives
+       endif
+      enddo
+
+      if (itap.ge.30) then
+        print*,'hq(l),ht(l),dq,dt,q:'
+     :  ,hq(16),ht(16),d_q_cvg(16),d_t_cvg(16),q(16,1)
+      endif
+
+      ELSE
+
+      do l = 1, nlayer
+       d_t_cvg(l) = 0.
+       d_q_cvg(l) = 0.
+      enddo
+
+      ENDIF
+
+      do l = 1, nlayer
+       dt_dyn(l) = d_t_cvg(l) / timestep
+       dq_dyn(l,1) = d_q_cvg(l) / timestep
+       dq_dyn(l,2) = 0.0
+       dq_dyn(l,3) = 0.0
+      enddo
+
+
+c   calcul des tendances physiques:
+c   -------------------------------
+
+      if (itsourcecont) then
+        do iq=5,nqmx
+          ll=min(iq-4,18)
+          if(firstcall) then 
+c          q(ll,iq)=5/((plev(ll)-plev(ll+1))/rg)
+c           q(ll,iq)=5 
+c          q(ll,iq)=5*play(ll)/((plev(ll)-plev(ll+1))/rg)/(RD*temp(ll))
+            q(ll,iq)=5*play(ll)/(RD*temp(ll))
+          endif
+        print*,'q ll iq',q(ll,iq),ll,iq
+        enddo
+      else
+        do iq=5,nqmx
+          ll=min(iq-4,18)
+          q(ll,iq)=q(ll,iq)+1
+             do l=1,llm
+               q(l,iq)=q(l,iq)*(1-timestep/tho)
+        print*,'q ll iq',q(ll,iq),ll,iq
+             enddo   
+        enddo
+ 
+      endif
+
+c      CALL physiq(ecritphy,ngrid,nlayer,nqmx,
+c     :            firstcall,lastcall,
+c     :            day,day,time,timestep,
+c     :            plev,play,phi,phis,paire,
+c     :            u,v,temp,q,
+c     :            du_dyn, dv_dyn,  dt_dyn, dq_dyn,
+c     :            w,
+c     :            du,dv,dt,dq,dpsrf)
+
+c     CALL physiq (ecritphy,ngrid,nlayer,nqmx,firstcall,lastcall,
+c    ,   day,day,time,timestep,plev,play,phi,phis,paire,
+c     ,             presnivs,clesphy0,u,v,temp,q,  
+c pour calculer proprement la tendance de cape liee a la
+c dynamique, il faut entrer aussi temp0 et q0:
+c    ,             presnivs,clesphy0,u,v,temp,q,temp0,q0,  
+c    ,             du_dyn,dv_dyn,dt_dyn,dq_dyn,w,
+c - sorties
+c    s             du,dv,dt,dq,dpsrf)
+
+      print*,'PAS DE TEMPS ',timestep
+      print*,'ATTENTION!!! Il faudra passer temp0 et q0'
+      CALL physiq(ngrid,nlayer,nqmx,
+     :            firstcall,lastcall,
+     :            day,day,time,timestep,
+     :            plev,play,phi,phis,paire,presnivs,clesphy0,
+     :            u,v,temp,q,
+     :            w,
+     :            du,dv,dt,dq,dpsrf)
+
+
+
+
+      firstcall=.false.
+
+
+c   Ajout des tendances
+c   -------------------
+      DO l=1,nlayer
+         u(l)=u(l)+timestep*du(l)
+         v(l)=v(l)+timestep*dv(l)
+
+         temp(l)=temp(l)+timestep*dt(l)
+      ENDDO
+
+      do iq=1,nqmx
+         do l=1,nlayer
+            q(l,iq)=q(l,iq)+timestep*dq(l,iq)
+         enddo
+      enddo
+c      write(78,rec=itap) ((q(l,iq),l=1,llm),iq=5,22)
+      write(78,rec=itap) ((q(l,iq),iq=5,22),l=1,llm)
+      itap = itap + 1
+
+
+         time=time+timestep/rday
+         it=it+1
+         if(mod(it,2000).eq.0) print*,'TIME=',time
+         if(time.gt.fnday) then
+            CALL endg1d(1,nlayer,play,int(time/ecritphy),timestep)
+            stop
+         endif
+
+         GOTO 1
+
+      END
+
+c=======================================================================
+c=======================================================================
+c    FIN DU PROGRAMMES
+c    CI-DESSOUS, QUELQUES SOUS-PROGRAMMES UTILS
+c=======================================================================
+c=======================================================================
+
+#include "1DUTILS.h"
+#include "1Dconv.h"
