program create_netcdf

      IMPLICIT NONE

include "/d3/crilmd/LMDZ20091210.trunk/netcdf-4.0.1/include/netcdf.inc"
                 
! lm=nb niveaux atm., ls= nb niveaux dans sol
! nbt= nb de tendances
! nbstep= nb de tendances
! nbstep2= nb de flux
! nbv= nb variables dans amma.nc
      INTEGER lm,lm2,ls,nbstep,nbstep2,nbt
      PARAMETER (lm=29,lm2=29,ls=1,nbstep=7,nbstep2=48,nbt=9)  
      
      INTEGER k,l,jk,ls2,psurf2,nb_forc,lev_forc
!

!temps
! nbstep = nb de forcages
! nbstep2= nb de flux
      REAL*8 time(nbstep),heure(nbstep)
      REAL*8 time2(nbstep2),hflux(nbstep2)
      REAL*8 time_out(nbstep2)
      INTEGER an,mois,jour
      

!niveaux verticaux
      REAL*8 z(lm)
!     REAL*8 z1(lm),zz(lm2,nbstep),zz2(lm)
      REAL*8 z1(lm),zz2(lm)
      REAL*8 psurf
!     PARAMETER(psurf=98800.)
      REAL*8 pres(lm)

!champs initiaux
      REAL*8 theta(lm),qv(lm),u(lm),v(lm)
      REAL*8 pp(lm),rv(lm),Temp(lm)
      REAL*8 z0,z0t,lati,long,pclay,psand,tg(ls+1),sm(ls+1)
      
!flux
      REAL*8 sens(nbstep2),flat(nbstep2)

!forcages grande echelle
!     REAL*8 pp1(lm2,nbstep),vitu(lm2,nbstep),vitv(lm2,nbstep)
!     REAL*8 vitw(lm2,nbstep),dth(lm2,nbstep),drv(lm2,nbstep)
!     REAL*8 dt(lm2,nbstep),dq(lm2,nbstep)
!
!     REAL*8 pp2(lm2,nbstep),vitu2(lm2,nbstep),vitv2(lm2,nbstep)
!     REAL*8 vitw2(lm2,nbstep),dth2(lm2,nbstep),drv2(lm2,nbstep)
!     REAL*8 dt2(lm2,nbstep),dq2(lm2,nbstep)
!
! tend1(niveau,temps,nb param)
! les parametres sont: z,pres,du,dv,dw,dth,drv,dt,dq
! tend2 est identiques sf qu il est a la meme grille temporelle que les flux
! tend_final a la meme grille vertical que les profils initiaux
      REAL*8 tend1(lm2,nbstep,nbt)
      REAL*8 tend2(lm2,nbstep2,nbt)
      REAL*8 tend_final(lm2,nbstep2,nbt)

      INTEGER tend_u,tend_v,tend_w,tend_t,tend_q
      INTEGER nudg_u,nudg_v,nudg_w,nudg_t,nudg_q
      

      REAL*8 XLVTT
      PARAMETER(XLVTT=2.5008E+6)

      integer nbvar3d_out
      parameter (nbvar3d_out=27)
      character (len=50), dimension(nbvar3d_out) :: varname3d_out
      character*5 type_sol
      character*4 poub

      integer ierr
      
      integer :: altdimout,timedimout
!,timevarout
      integer :: nout

      integer var3didout(27),toto(lm)

!     Read in data
!     read initial profiles

      open(10,file='profil_init.txt')
      read(10,*) psurf2,l
      IF (l.NE.lm) stop('Pas le bon nombre de niveaux lm !')
      do k=1,lm
        read(10,*),z(k),pp(k),theta(k),rv(k),u(k),v(k),Temp(k),qv(k)
      enddo
      read (10,*) z0,z0t
      read (10,*) lati,long
      read (10,*) type_sol
      read (10,*) pclay,psand
      read (10,*) ls2
      IF (ls2.NE.ls) stop('Pas le bon nombre de niveaux ls !')
      do k=1,ls+1
        read (10,*) tg(k),sm(k)
      enddo
      close(10)

!     print *,'Lecture de profil_init.txt: psol lev_init',psurf2,lm2
!     do k=1,lm2
!       print *,'profils',k,z(k),pp(k),theta(k),rv(k),u(k),v(k),Temp(k),qv(k)
!     enddo
!     print *,'z0 z0t ',z0,z0t
!     print *,'lat lon ',lati,long
!     print *,'type_sol ',type_sol
!     print *,'pclay psand ',pclay,psand
!     print *,'lev_sol ',ls2
!     do k=1,ls+1
!       print *,'tg sm',k,tg(k),sm(k)
!     enddo


!read tendencies
!-------------------------------
! TEND A B C D E
! A=1 on utilise tendances sur U
! B=1 on utilise tendances sur V
! C=1 on utilise tendances sur W
! D=1 on utilise tendances sur theta
! E=1 on utilise tendances sur vapeur eau
!-------------------------------
! NUDG A B C D E
! A= U est nudge avec un temps de relaxation de A
! B= V est nudge avec un temps de relaxation de B
! C= W est nudge avec un temps de relaxation de C
! D= theta est nudge avec un temps de relaxation de D
! E= vapeur d'eau est nudgee avec un temps de relaxation de E
! si A,B,C,D,E=0 => pas de nudging
! si A,B=-1 => nudging avec vent geostrophique fourni
!-------------------------------
      open(11,file='profil_tend.txt')
! Lecture des flags de tendances
      read(11,*),poub,tend_u,tend_v,tend_w,tend_t,tend_q
      IF (poub.NE.'TEND') stop('Erreur fichier profil_tend.txt ! manque TEND')
!     print *,'TEND:',poub,tend_u,tend_v,tend_w,tend_t,tend_q
! Lecture des flags de nudging
      read(11,*),poub,nudg_u,nudg_v,nudg_w,nudg_t,nudg_q
      IF (poub.NE.'NUDG') stop('Erreur fichier profil_tend.txt ! manque NUDG')
!     print *,'NUDG:',poub,nudg_u,nudg_v,nudg_w,nudg_t,nudg_q
! Combien de profils de forcing ?
      read(11,*),nb_forc
      IF (nb_forc.NE.nbstep) stop('Erreur fichier profil_tend.txt ! nb_forc NE nbstep')
!     print *,'nb_forc=',nb_forc
! Combien de niveaux ?
      read(11,*),lev_forc
!     print *,'lev_forc=',lev_forc
! Lecture du profil lui-meme
! Detail sur le profil qui suit
      do l=1,nbstep
      read(11,*),poub,an,mois,jour,heure(l)
!     print *,'temps:',an,mois,jour,heure(l)/3600.
        do k=1,lm2
!       read(11,*),zz(k,l),pp1(k,l),vitu(k,l),vitv(k,l),vitw(k,l),dth(k,l),drv(k,l),dt(k,l),dq(k,l)
        read(11,*),tend1(k,l,:)
!       print *,'tendances= ',l,k,tend1(k,l,:)
        enddo
      enddo
      close(11)

!flux definition       
      open(12,file='flux.txt')
      read(12,*) poub
      do k=1,nbstep2
!     read(12,*),hflux(k),flat(k),sens(k) ! fichiers C.Rio
      read(12,*),hflux(k),sens(k),flat(k) ! fichiers Fleur
!     print *,'Flux lus: sens flat',hflux(k)/3600.,sens(k),flat(k)
      enddo
      close(12)

!definition de l'axe des temps
      do l=1,nbstep2
         time2(l)=10.+float(l)*30./60./24.
!        print *,'l time2=',l,time2(l)
      enddo
      do l=1,nbstep
         time(l)=10.25+float(l)*3./24.
!        print *,'l time=',l,time(l)
      enddo


! interpolation en temps des tendances
! sur la meme grille temporelle que les flux
! 12,18,24... sont exprimes en demi heures et correspondent a
! 6,9,12... heures
      do l=1,nbstep
         do k=1,lm
           tend2(k,l,:)=0.
         enddo
      enddo
      do k=1,lm
       do jk=1,nbt
         tend2(k,12,jk)=tend1(k,1,jk)  ! 6h
         tend2(k,18,jk)=tend1(k,2,jk)  ! 9h
         tend2(k,24,jk)=tend1(k,3,jk)  ! 12h
         tend2(k,30,jk)=tend1(k,4,jk)  ! 15h
         tend2(k,36,jk)=tend1(k,5,jk)  ! 18h
         tend2(k,42,jk)=tend1(k,6,jk)  ! 21h
         tend2(k,48,jk)=tend1(k,7,jk)  ! 24h
       enddo
      enddo

      do l=13,17
        do k=1,lm
          do jk=1,nbt
            tend2(k,l,jk)=(tend2(k,18,jk)-tend2(k,12,jk))/6.*float(l)+(18.*tend2(k,12,jk)-12.*tend2(k,18,jk))/6.
          enddo
        enddo
      enddo
      do l=19,23
         do k=1,lm
          do jk=1,nbt
            tend2(k,l,jk)=(tend2(k,24,jk)-tend2(k,18,jk))/6.*float(l)+(24.*tend2(k,18,jk)-18.*tend2(k,24,jk))/6.
          enddo
         enddo
      enddo
      do l=25,29
         do k=1,lm
          do jk=1,nbt
            tend2(k,l,jk)=(tend2(k,30,jk)-tend2(k,24,jk))/6.*float(l)+(30.*tend2(k,24,jk)-24.*tend2(k,30,jk))/6.
          enddo
         enddo
      enddo
      do l=31,35
         do k=1,lm
          do jk=1,nbt
            tend2(k,l,jk)=(tend2(k,36,jk)-tend2(k,30,jk))/6.*float(l)+(36.*tend2(k,30,jk)-30.*tend2(k,36,jk))/6.
          enddo
         enddo
      enddo
      do l=37,41
         do k=1,lm
          do jk=1,nbt
            tend2(k,l,jk)=(tend2(k,42,jk)-tend2(k,36,jk))/6.*float(l)+(42.*tend2(k,36,jk)-36.*tend2(k,42,jk))/6.
          enddo
         enddo
      enddo
      do l=43,47
         do k=1,lm
          do jk=1,nbt
            tend2(k,l,jk)=(tend2(k,48,jk)-tend2(k,42,jk))/6.*float(l)+(48.*tend2(k,42,jk)-42.*tend2(k,48,jk))/6.
          enddo
         enddo
      enddo
      do jk=12,48
      call interp_vertical(lm,z,pp&
  ,lm,tend2(:,jk,1),tend2(:,jk,2),tend2(:,jk,3),tend2(:,jk,4)&
  ,tend2(:,jk,5),tend2(:,jk,6),tend2(:,jk,7),tend2(:,jk,8),tend2(:,jk,9)&
  ,tend_final(:,jk,1),tend_final(:,jk,2),tend_final(:,jk,3),tend_final(:,jk,4)&
  ,tend_final(:,jk,5),tend_final(:,jk,6),tend_final(:,jk,7),tend_final(:,jk,8),tend_final(:,jk,9))
      enddo
      do jk=1,lm
!      print *,'jk z pp z_tend p_tend z_final p_final',jk,z(jk),pp(jk),tend2(jk,12,1),tend2(jk,12,2),&
!      tend_final(jk,12,1),tend_final(jk,12,2)
       print *,'jk pp vitw p_final dw',jk,tend2(jk,28,2),tend2(jk,28,5),tend_final(jk,28,2),tend_final(jk,28,5)
      enddo

      do l=1,nbstep2
!        time_out(l)=(time2(l)-10)*86400.
         time_out(l)=time2(l)
      enddo
!on ecrit la pression en Pa
      do k=1,lm
         pp(k)=pp(k)*100.
      enddo


!ecriture des resultats dans un fichier nc

         print*,'avant ecriture ok'
 
         call initiate ("amma.nc",z,time_out,&
                        altdimout,timedimout,nbstep2,lm)

         print*,'apres initiate'
         call def_var(nout,"zz","z","m",1,(/altdimout/),var3didout(1),ierr)
         call def_var(nout,"pp","pp","hPa",1,(/altdimout/),var3didout(2),ierr)
         call def_var(nout,"theta","theta","K",1,(/altdimout/),var3didout(3),ierr)
         call def_var(nout,"rv","rv","kg/kg",1,(/altdimout/),var3didout(4),ierr)
         call def_var(nout,"u","u","m/s",1,(/altdimout/),var3didout(5),ierr)
         call def_var(nout,"v","v","m/s",1,(/altdimout/),var3didout(6),ierr)
         call def_var(nout,"temp","temp","K",1,(/altdimout/),var3didout(7),ierr)
         call def_var(nout,"qv","qv","kg/kg",1,(/altdimout/),var3didout(8),ierr)
         call def_var(nout,"sens","sens","W/m2",1,(/timedimout/),var3didout(9),ierr)
         call def_var(nout,"flat","flat","W/m2",1,(/timedimout/),var3didout(10),ierr)
         call def_var(nout,"du","du","m/s",2,(/altdimout,timedimout/),var3didout(11),ierr)
         call def_var(nout,"dv","dv","m/s",2,(/altdimout,timedimout/),var3didout(12),ierr)
         call def_var(nout,"dw","dw","m/s",2,(/altdimout,timedimout/),var3didout(13),ierr)
         call def_var(nout,"dth","dth","K/s",2,(/altdimout,timedimout/),var3didout(14),ierr)
         call def_var(nout,"drv","drv","kg/kg/s",2,(/altdimout,timedimout/),var3didout(15),ierr)
         call def_var(nout,"dt","dt","K/s",2,(/altdimout,timedimout/),var3didout(16),ierr)
         call def_var(nout,"dq","dq","kg/kg/s",2,(/altdimout,timedimout/),var3didout(17),ierr)
         call def_var(nout,"tend_u","tend_u","-",1,(/altdimout/),var3didout(18),ierr)
         call def_var(nout,"tend_v","tend_v","-",1,(/altdimout/),var3didout(19),ierr)
         call def_var(nout,"tend_w","tend_w","-",1,(/altdimout/),var3didout(20),ierr)
         call def_var(nout,"tend_t","tend_t","-",1,(/altdimout/),var3didout(21),ierr)
         call def_var(nout,"tend_q","tend_q","-",1,(/altdimout/),var3didout(22),ierr)
         call def_var(nout,"nudg_u","nudg_u","-",1,(/altdimout/),var3didout(23),ierr)
         call def_var(nout,"nudg_v","nudg_v","-",1,(/altdimout/),var3didout(24),ierr)
         call def_var(nout,"nudg_w","nudg_w","-",1,(/altdimout/),var3didout(25),ierr)
         call def_var(nout,"nudg_t","nudg_t","-",1,(/altdimout/),var3didout(26),ierr)
         call def_var(nout,"nudg_q","nudg_q","-",1,(/altdimout/),var3didout(27),ierr)

! A FAIRE: ajouter a la suite les scalaires qui sont fournis par profil_init.txt
    
#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(1),z)
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(1),z)
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_1"
         endif
 
#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(2),pp)
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(2),pp)
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_2"
         endif

#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(3),theta)
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(3),theta)
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_3"
         endif

#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(4),rv)
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(4),rv)
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_4"
         endif

#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(5),u)
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(5),u)
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_5"
         endif
     
#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(6),v)
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(6),v)
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_6"
         endif

#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(7),Temp)
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(7),Temp)
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_7"
         endif

#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(8),qv)
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(8),qv)
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_8"
         endif

#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(9),sens)
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(9),sens)
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_9"
         endif

#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(10),flat)
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(10),flat)
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_10"
         endif

        do jk=3,nbt
#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(10+jk-2),tend_final(:,:,jk))
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(10+jk-2),tend_final(:,:,jk))
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_tendances"
         endif
        enddo

         toto(:)=tend_u
         ierr= NF_PUT_VAR_INT(nout,var3didout(18),toto)
         toto(:)=tend_v
         ierr= NF_PUT_VAR_INT(nout,var3didout(19),toto)
         toto(:)=tend_w
         ierr= NF_PUT_VAR_INT(nout,var3didout(20),toto)
         toto(:)=tend_t
         ierr= NF_PUT_VAR_INT(nout,var3didout(21),toto)
         toto(:)=tend_q
         ierr= NF_PUT_VAR_INT(nout,var3didout(22),toto)
         toto(:)=nudg_u
         ierr= NF_PUT_VAR_INT(nout,var3didout(23),toto)
         toto(:)=nudg_v
         ierr= NF_PUT_VAR_INT(nout,var3didout(24),toto)
         toto(:)=nudg_w
         ierr= NF_PUT_VAR_INT(nout,var3didout(25),toto)
         toto(:)=nudg_t
         ierr= NF_PUT_VAR_INT(nout,var3didout(26),toto)
         toto(:)=nudg_q
         ierr= NF_PUT_VAR_INT(nout,var3didout(27),toto)

ierr=NF_CLOSE(nout)       
       
!
      

       contains
!**********************************************************
Subroutine initiate (filename,z,time,&
                     altdimout,timedimout,&
                     nj,lm)

implicit none

include "/d3/crilmd/LMDZ20091210.trunk/netcdf-4.0.1/include/netcdf.inc"

character (len=*), intent(in):: filename
real*8, dimension(:), intent(in):: z,time
!integer, intent(out):: nout,timevarout
integer, intent(out):: altdimout,timedimout

integer :: altdim,timedim
integer :: newvarid
integer :: nvarid,ierr,i,kloop,nj,lm
!real, dimension(nj), intent(out) :: time

write(*,*) "creating "//trim(adjustl(filename))//'...'
ierr = NF_CREATE(filename,NF_CLOBBER, nout)
if (ierr.NE.NF_NOERR) THEN
   WRITE(*,*)'ERROR: Impossible to create the file.'
   WRITE(*,*) NF_STRERROR(ierr)
   stop ""
endif

!ierr = NF_DEF_DIM (nout, "lat", size(lat), latdimout)
!ierr = NF_DEF_DIM (nout, "lon", size(lon), londimout)
ierr = NF_DEF_DIM (nout, "lev", lm, altdimout)
ierr = NF_DEF_DIM (nout, "time", nj, timedimout)

ierr = NF_ENDDEF(nout)

 call def_var(nout,"lev","lev","m",1,(/altdimout/),nvarid,ierr)
#ifdef NC_DOUBLE
ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,z(1:lm))
#else
ierr = NF_PUT_VAR_REAL (nout,nvarid,z(1:lm))
#endif

 call def_var(nout,"time","time","days since 2006-07-10 00:00:00",1,(/timedimout/),nvarid,ierr)
!ierr = NF_PUT_VAR_INT (nout,timevarout,time(1:nj))
#ifdef NC_DOUBLE
!ierr = NF_PUT_VAR_DOUBLE (nout,timevarout,time(1:nj))
ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,time(1:nj))
#else
!ierr = NF_PUT_VAR_REAL (nout,timevarout,time(1:nj))
ierr = NF_PUT_VAR_REAL (nout,nvarid,time(1:nj))
#endif


end subroutine initiate

!*************************************************************
subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)

implicit none

include "/d3/crilmd/LMDZ20091210.trunk/netcdf-4.0.1/include/netcdf.inc"

character (len=*) :: title,units,name
integer :: nid,nbdim,nvarid,ierr
integer, dimension(nbdim) :: dim

ierr=NF_REDEF(nid)
#ifdef NC_DOUBLE
ierr = NF_DEF_VAR (nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid)
#else
ierr = NF_DEF_VAR (nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
#endif
if(ierr/=NF_NOERR) then
   write(*,*) NF_STRERROR(ierr)
   stop "in def_var"
endif
ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", len_trim(adjustl(title)),adjustl(title))
if(ierr/=NF_NOERR) then
   write(*,*) NF_STRERROR(ierr)
   stop "in def_var"
endif
ierr = NF_PUT_ATT_TEXT (nid, nvarid, "units", len_trim(adjustl(units)),adjustl(units))
if(ierr/=NF_NOERR) then
   write(*,*) NF_STRERROR(ierr)
   stop "in def_var"
endif
ierr = NF_ENDDEF(nid)

end subroutine def_var

!*************************************************************
subroutine catchaxis(nid,iim,jjm,llm,nj,xc,yc,zc,time,ierr)

!include "/home/distrib/local/netcdf-3.6.1/include/netcdf.inc"
include "/d3/crilmd/LMDZ20091210.trunk/netcdf-4.0.1/include/netcdf.inc"
integer, intent(in) :: nid,iim,jjm,llm,nj
real*8, dimension(iim), intent(out) :: xc
real*8, dimension(jjm), intent(out) :: yc
real*8, dimension(llm), intent(out) :: zc
real*8, dimension(nj), intent(out) :: time
integer, intent(out) :: ierr

integer :: i
integer :: latvar,lonvar,altvar,timevar
integer :: latlen,lonlen,altlen,timelen
integer :: londimin,latdimin,altdimin,timedimin

! Control & lecture on dimensions
! ===============================
   ierr=NF_INQ_DIMID(nid,"yc",latdimin)
   ierr=NF_INQ_VARID(nid,"yc",latvar)
   if (ierr.NE.NF_NOERR) then
      write(*,*) 'ERROR: Field <lat> is missing'
      stop ""  
   endif
   ierr=NF_INQ_DIMLEN(nid,latdimin,latlen)

   ierr=NF_INQ_DIMID(nid,"xc",londimin)
   ierr=NF_INQ_VARID(nid,"xc",lonvar)
   if (ierr.NE.NF_NOERR) then
      write(*,*) 'ERROR: Field <lon> is lacking'
      stop "" 
   endif
  ierr=NF_INQ_DIMLEN(nid,londimin,lonlen)

   ierr=NF_INQ_DIMID(nid,"zc",altdimin)
   ierr=NF_INQ_VARID(nid,"zc",altvar)
   if (ierr.NE.NF_NOERR) then
      write(*,*) 'ERROR: Field <presnivs> is lacking'
      stop ""
   endif
   ierr=NF_INQ_DIMLEN(nid,altdimin,altlen)

!test:on rajoute la dimension temps
   ierr=NF_INQ_DIMID(nid,"time",timedimin)
   ierr=NF_INQ_VARID(nid,"time",timevar)
   if (ierr.NE.NF_NOERR) then
      write(*,*) 'ERROR: Field <time> is lacking'
      stop ""
   endif
   ierr=NF_INQ_DIMLEN(nid,timedimin,timelen)

   if((latlen/=jjm).or.(lonlen/=iim).or.(altlen/=llm)) then
      write(*,*) 'ERROR: Not the good lenght for axis'
      write(*,*) 'longitude: ',lonlen,iim+1
      write(*,*) 'latitude: ',latlen,jjm
      write(*,*) 'presniv: ',altlen,llm
!      stop ""  
   endif

#ifdef NC_DOUBLE
   ierr = NF_GET_VAR_DOUBLE(nid,latvar,yc)
   ierr = NF_GET_VAR_DOUBLE(nid,lonvar,xc)
   ierr = NF_GET_VAR_DOUBLE(nid,altvar,zc)
   ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
#else
   ierr = NF_GET_VAR_REAL(nid,latvar,yc)
   ierr = NF_GET_VAR_REAL(nid,lonvar,xc)
   ierr = NF_GET_VAR_REAL(nid,altvar,zc)
   ierr = NF_GET_VAR_REAL(nid,timevar,time)
#endif

end subroutine catchaxis


!*************************************************************

!*************************************************************
subroutine giveatt(nid,mean)

Implicit none

!include "/home/distrib/local/netcdf-3.6.1/include/netcdf.inc"
include "/d3/crilmd/LMDZ20091210.trunk/netcdf-4.0.1/include/netcdf.inc"
integer nid,ierr
character (len=100):: att
logical,intent(in) :: mean

ierr = NF_REDEF (nid)
att="GDT 1.3"
ierr= NF_PUT_ATT_TEXT(nid,NF_GLOBAL,'conventions',len_trim(att),att)
if(ierr/=NF_NOERR) then
   write(*,*) NF_STRERROR(ierr)
   stop "in giveatt"
endif

if (mean) then
   att="Mars Climate Database v4.0 - Mean variables"
else
   att="Mars Climate Database v4.0 - Standard deviations"
endif

ierr= NF_PUT_ATT_TEXT(nid,NF_GLOBAL,'history',len_trim(att),att)
if(ierr/=NF_NOERR) then
   write(*,*) NF_STRERROR(ierr)
   stop "in giveatt"
endif

att="LMD-AOPP-ESA-CNES"
ierr= NF_PUT_ATT_TEXT(nid,NF_GLOBAL,'institution',len_trim(att),att)
if(ierr/=NF_NOERR) then
   write(*,*) NF_STRERROR(ierr)
   stop "in giveatt"
endif
ierr = NF_ENDDEF (nid)

end subroutine giveatt

!*************************************************************
subroutine missing_value(nout,nvarid,missing)
IMPLICIT NONE

! useful to watch results with ferret

!include "/home/distrib/local/netcdf-3.6.1/include/netcdf.inc"    
include "/d3/crilmd/LMDZ20091210.trunk/netcdf-4.0.1/include/netcdf.inc"
                                                                                   
INTEGER :: nout,nvarid,ierr
REAL :: missing
ierr = NF_REDEF (nout)
#ifdef NC_DOUBLE
ierr= NF_PUT_ATT_DOUBLE(nout,nvarid,'missing_value',NF_DOUBLE,1,missing)
#else
ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
#endif
iF (ierr.NE.NF_NOERR) THEN
     PRINT*, 'anl_NC: missing value attribution failed'
     WRITE(*,*) 'NF_NOERR', NF_NOERR
     CALL abort
ENDIF                                
ierr=NF_ENDDEF(nout)
return
                                                    
end subroutine missing_value

!****************************************************************
subroutine interp_vertical(nlev_init,z_init,p_init&
  ,nlev_tend,z_tend,p_tend,du_tend,dv_tend,dw_tend&
  ,dth_tend,drv_tend,dt_tend,dq_tend&
  ,z_final,p_final,du_final,dv_final,dw_final&
  ,dth_final,drv_final,dt_final,dq_final)

IMPLICIT NONE

!-----------------------------------------------------------------------------------
! Vertical interpolation of tendencies forcing profiles onto initial profiles levels
!-----------------------------------------------------------------------------------
! profil_init.txt nb of levels=nlev_init
! contents:
! z_init, p_init, theta_init, rv_init,u_init,v_init,temp_init,qv_init
!-----------------------------------------------------------------------------------
! profil_tend.txt nb of levels=nlev_tend
!                 nb of forcing times= nbt
! contents:
! z_tend, p_tend, du_tend, dv_tend, dw_tend, dth_tend, drv_tend, dt_tend,dq_tend
!-----------------------------------------------------------------------------------

       INTEGER nlevmax
       PARAMETER (nlevmax=41)
       INTEGER nlev_init,nlev_tend,mxcalc

       REAL*8 z_init(nlev_init), p_init(nlev_init)
       REAL*8 z_tend(nlev_tend), p_tend(nlev_tend)
       REAL*8 du_tend(nlev_tend),dv_tend(nlev_tend)
       REAL*8 dw_tend(nlev_tend),dth_tend(nlev_tend)
       REAL*8 drv_tend(nlev_tend),dt_tend(nlev_tend)
       REAL*8 dq_tend(nlev_tend)

       REAL*8 z_final(nlev_tend), p_final(nlev_tend)
       REAL*8 du_final(nlev_tend),dv_final(nlev_tend)
       REAL*8 dw_final(nlev_tend),dth_final(nlev_tend)
       REAL*8 drv_final(nlev_tend),dt_final(nlev_tend)
       REAL*8 dq_final(nlev_tend)

       INTEGER l,k,k1,k2,kp
       REAL*8 aa,frac,frac1,frac2,fact

       do l = 1, nlev_init
!      print *,'p_init,p_tend',l,p_init(l),p_tend(l)

        if (p_init(l).ge.p_tend(nlev_tend)) then

         k1=0
         k2=0

         if (p_init(l).le.p_tend(1)) then

         do k = 1, nlev_tend-1
          if (p_init(l).le.p_tend(k) .and. p_init(l).ge.p_tend(k+1)) then
            k1=k
            k2=k+1
!           print *,'p_init,p_tend,l,k,k1,k2',l,k,k1,k2,p_init(l),p_tend(k),p_tend(k+1)
          endif
         enddo

         if (k1.eq.0 .or. k2.eq.0) then
          write(*,*) 'PB! k1, k2 = ',k1,k2
          write(*,*) 'l,p_init(l) = ',l,p_init(l)
         do k = 1, nlev_tend-1
          write(*,*) 'k,p_tend(k) = ',k,p_tend(k)
         enddo
         endif

         frac = (p_tend(k2)-p_init(l))/(p_tend(k2)-p_tend(k1))
         p_final(l)  = p_init(l)
         z_final(l)  = z_tend(k2)   - frac*(z_tend(k2)-z_tend(k1))
         du_final(l) = du_tend(k2)  - frac*(du_tend(k2)-du_tend(k1))
         dv_final(l) = dv_tend(k2)  - frac*(dv_tend(k2)-dv_tend(k1))
         dw_final(l) = dw_tend(k2)  - frac*(dw_tend(k2)-dw_tend(k1))
         dth_final(l)= dth_tend(k2) - frac*(dth_tend(k2)-dth_tend(k1))
         drv_final(l)= drv_tend(k2) - frac*(drv_tend(k2)-drv_tend(k1))
         dt_final(l) = dt_tend(k2)  - frac*(dt_tend(k2)-dt_tend(k1))
         dq_final(l) = dq_tend(k2)  - frac*(dq_tend(k2)-dq_tend(k1))

         else !p_init>p_tend(1)

         k1=1
         k2=2
         frac1 = (p_init(l)-p_tend(k2))/(p_tend(k1)-p_tend(k2))
         frac2 = (p_init(l)-p_tend(k1))/(p_tend(k1)-p_tend(k2))
         p_final(l)  = p_init(l)
         z_final(l)  = frac1*z_tend(k1)   - frac2*z_tend(k2)
         du_final(l) = frac1*du_tend(k1)  - frac2*du_tend(k2)
         dv_final(l) = frac1*dv_tend(k1)  - frac2*dv_tend(k2)
         dw_final(l) = frac1*dw_tend(k1)  - frac2*dw_tend(k2)
         dth_final(l)= frac1*dth_tend(k1) - frac2*dth_tend(k2)
         drv_final(l)= frac1*drv_tend(k1) - frac2*drv_tend(k2)
         dt_final(l) = frac1*dt_tend(k1)  - frac2*dt_tend(k2)
         dq_final(l) = frac1*dq_tend(k1)  - frac2*dq_tend(k2)

         endif ! p_init<p_tend(1)
       endif

       enddo 

!      do l = 1,nlev_tend
!      print *,'t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l) ',
!    $        l,t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l)
!      enddo

return
end subroutine interp_vertical

!****************************************************************************

      end