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