!
! $Header$
!
!-------------------------------------------------------------------
! Make relationship between any variable and w500 at resolution 2.5x2.5
! -> mean seasonal cycle of the relationship over the period 1985-1989
! -> interannual anomalies (by removing the mean seasonal cycle)
! -> decomposition of interannual anomalies into thermo, dynam and mix components
!
! Specify:
!	- undefined value of the data
!	- min/max/step of the w500 bins
!	- nb of months in timeseries
!	- montly index of January 1985
!
! Sandrine Bony, March 2004
! Repris pour faire du netcdf, L. Fairhead 2005/02
!-------------------------------------------------------------------

  implicit none

#include "netcdf.inc"

  integer imax,jmax,an,mmax,nbinmax,noc
  parameter (imax=360,jmax=180,mmax=1700,an=12,nbinmax=200,noc=5)

  integer nrun
  parameter (nrun=1)

  real zmsk(imax,jmax)
  real surfmax,zsurf(imax,jmax)
  logical ok_msk(imax,jmax)

  character*20 namegcm, run

  real surftest
  real lat0, undef_x
!  parameter (lat0 = 30.0) ! limits tropical belt
  parameter (undef_x = 1.e+20 ) ! undefined value

  real, allocatable, dimension(:,:) :: surf, msk, w, x
  real, allocatable, dimension(:,:) :: x_binw
  real, allocatable, dimension(:,:) :: w_binw
  real, allocatable, dimension(:,:) :: xpdf,nx_binw,sx_binw
  real, allocatable, dimension(:)   :: xpdfmean

  integer nb_bin, nbr_inf, nbr_sup
  real min_bin, max_bin, step_bin, w_mult, x_mult
  real xpdftot(mmax)
  real xpdfmeantot



  INTEGER i,j,m,ir,im,jm,lm,itime,n, il, l
  real undef, undef_w, pi, x1, w1, msk1
  integer :: iostat
  character (len=512)        :: line_read
  character (len=20)         :: first_part
  character (len=128)       :: second_part  
  character cnat*4, bintype*4

! netcdf stuff
  real         :: time
  real, dimension(:), allocatable :: var_dim
  INTEGER      :: ierr, nvdims, lunout, in_file_id, ierrs, ndims, natts
  integer      :: lonid, latid, in_time_id, varid, out_file_id, ngatts, levid
  integer      :: nbin_id, var_bin_id, time_id, var_time_id, in_var_time_id
  INTEGER      :: var_pdf_id, var_nx_id, var_w_id, var_x_id, outlev_id
  integer      :: var_lev_id
  integer, dimension(4)  :: vdims, start, count
  character (len=80) :: long_name, varname, attname
  character (len=132) :: in_file
  real, dimension(:), allocatable ::lat

  logical :: var_3d = .false.

  integer :: iargc
  external iargc  
!---------------------------------------------------------------------
!-- General specifications :
!---------------------------------------------------------------------

!******************************* FIN INTERFACE *********************
! Lecture du fichier de configuration
  open (20, IOSTAT=iostat, file='config.def',form='formatted')
  if (iostat /= 0) THEN
    write(*,*)'Erreur ouverture du fichier config.def'
    stop
  endif

  config: do
    read(20,'(A)',iostat=iostat)line_read
    if (iostat /= 0) exit
    line_read = trim(line_read)
    IF (INDEX(line_read, '#') /= 1) THEN
      first_part = trim(line_read(1:INDEX(line_read, '=')-1))
      second_part = trim(line_read(INDEX(line_read, '=')+1:))
      selectcase(first_part)
        case('bintype')
          bintype = trim(second_part)
        case('cnat')
          cnat = trim(second_part)
        case('lunout')
          read(second_part,'(i)') lunout
        case('lat0')
          read(second_part,'(f)') lat0
      end select
    endif
  enddo config
  if (iostat > 0) then
    write(lunout,*) &
 &  'Probleme de lecture du fichier config.def, iostat = ',iostat
    stop
  endif
  close(20)
  
! undefined value for the output:
  undef  = 999.999 ! for output only

!------------------------------------------------------------
! Definition of the w500 bins:
!------------------------------------------------------------

  if (bintype.eq.'w500') then
    min_bin = -200.0  
    max_bin =  200.0
    step_bin =  10.0
    step_bin =  5.0
    w_mult =  864.0 ! Pa/s -> hPa/day
!   w_mult =  1.0 ! Pa/s -> hPa/day
  endif

! scale factor for surface temperature, precip and variable x:
  x_mult =  1.0 

  undef_w = undef

!---------------------------------------------------------------------
!-- preliminaries:
!---------------------------------------------------------------------

  if (cnat.ne.'ocea'.and.cnat.ne.'land'.and.cnat.ne.'glob'.and.cnat.ne.'mixt')&
& then
    write(*,*) 'erreur cnat ',cnat
    stop
  endif

  nb_bin = (max_bin-min_bin)/step_bin + 1
  if (nb_bin.gt.nbinmax) then
    write(*,*) 'augmenter le nb de bins'
    write(*,*) 'nbinmax, nb_bin = ',nbinmax,nb_bin
    stop
  endif

! Il faut deux arguments a l appel: fichier et variable
  CALL getarg(1, in_file)
  if (iargc() == 0 .OR. iargc() /= 2 ) THEN
    write(lunout,*)' '
    write(lunout,*)' Utilisation de ce programme: '
    write(lunout,*)' ./geo2reg nom_de_fichier [variable]'
    write(lunout,*)                                    &
    &  '        ou nom_de_fichier est le nom du fichier a traiter'
    write(lunout,*)                                             &
    &  '        et variable la variable a traiter [optionel]'
       write(lunout,*)' '
    write(lunout,*)' ./geo2reg -h sort ce message'
    write(lunout,*)' '
    stop
  endif
  CALL getarg(2, varname)
 
! Ouverture du fichier a traiter (histmth)
  ierr = NF_OPEN(in_file, NF_NOwrite, in_file_id)
  if (ierr /= NF_NOERR) then
    write(lunout,*)NF_STRERROR(ierr)
    stop
  endif
!
! lire im, jm, lm et temps
  ierrs = 0
  ierr = NF_INQ_DIMID(in_file_id, 'lon', lonid)  
  ierrs = ierrs + ierr
  ierr = NF_INQ_DIMLEN(in_file_id, lonid, im)
  ierrs = ierrs + ierr 
  if (ierrs /= 2 * NF_NOERR) THEN
    write(lunout,*)'Pb. avec la lecture de la longitude'
    stop
  endif
  ierrs = 0
  ierr = NF_INQ_DIMID(in_file_id, 'lat', latid)  
  ierrs = ierrs + ierr
  ierr = NF_INQ_DIMLEN(in_file_id, latid, jm)
  ierrs = ierrs + ierr
  allocate (lat(jm)) 
  ierr = NF_INQ_VARID(in_file_id,'lat', latid)
  ierrs = ierrs + ierr
  ierr = NF_GET_VAR_REAL(in_file_id, latid, lat)
  ierrs = ierrs + ierr
  if (ierrs /= 4 * NF_NOERR) THEN
    write(lunout,*)'Pb. avec la lecture de la latitude'
    stop
  endif
  ierrs = 0
  ierr = NF_INQ_DIMID(in_file_id, 'presnivs', levid)
  ierrs = ierrs + ierr
  ierr = NF_INQ_DIMLEN(in_file_id, levid, lm)
  ierrs = ierrs + ierr 
  if (ierrs /= 2 * NF_NOERR) THEN
    write(lunout,*)'Pb. avec la lecture des niveaux verticaux'
    stop
  endif
  

  ierrs = 0
  ierr = NF_INQ_DIMID(in_file_id, 'time_counter', in_time_id)
  ierrs = ierrs + ierr
  ierr = NF_INQ_DIMLEN(in_file_id, in_time_id, itime)
  ierrs = ierrs + ierr 
  if (ierrs /= 2 * NF_NOERR) THEN
    write(lunout,*)'Pb. avec la lecture du temps'
    stop
  endif

! lecture de l aire des mailles
  allocate (surf(im,jm))
  ierrs = 0
  ierr = NF_INQ_VARID(in_file_id, 'aire', varid)
  ierrs = ierrs + ierr
  ierr = NF_GET_VAR_REAL(in_file_id, varid, surf)
  ierrs = ierrs + ierr
  if (ierrs /= 2 * NF_NOERR) THEN
    write(lunout,*)'Pb. avec la lecture de l''aire'
    stop
  endif
 
!-- compute gridbox areas:


  pi = acos(-1.)
  surfmax=0.
  surftest=0.0
  do i = 1, im
    do j = 1, jm
      surftest=surftest+surf(i,j)
      if(surf(i,j).gt.surfmax) surfmax=surf(i,j)
    enddo
  enddo
  write(*,*) 'sfce totale = 4pi ',surftest/pi

!---------------------------------------------------------------
! Lecture du masque
!
  allocate (msk(im,jm))
  if (cnat.ne.'glob') then
    ierrs = 0
    ierr = NF_INQ_VARID(in_file_id, 'pourc_ter', varid)
    ierrs = ierrs + ierr
    ierr = NF_GET_VAR_REAL(in_file_id, varid, msk)
    ierrs = ierrs + ierr
    if (ierrs /= 2 * NF_NOERR) THEN
      write(lunout,*)'Pb. avec la lecture du masque'
      stop
    endif
  else
    do i = 1, im
      do j = 1, jm
        msk(i,j) = undef
      enddo
    enddo
  endif


  do i = 1, im
    do j = 1, jm
      ok_msk(i,j) = .FALSE.
      if (cnat.eq.'ocea' .and. msk(i,j).le.0.01) ok_msk(i,j) = .TRUE.
      if (cnat.eq.'mixt' .and. msk(i,j).gt.0.01.and. msk(i,j).le.0.99) &
  &                                             ok_msk(i,j) = .TRUE.
      if (cnat.eq.'land' .and. msk(i,j).gt.0.99) ok_msk(i,j) = .TRUE.
      if (cnat.eq.'glob') ok_msk(i,j) = .TRUE.
    enddo
  enddo

!
! Pour savoir si la variable a traiter est 3D
!
  ierrs = 0
  ierr = NF_INQ_VARID(in_file_id, varname, varid)
  ierrs = ierrs + ierr
  ierrs = NF_INQ_VARNDIMS(in_file_id, varid, ndims)
  if (ierrs /= 2 * NF_NOERR) THEN
    write(lunout,*)'Pb. avec la lecture de ', varname
    stop
  endif
  if (ndims == 4) var_3d = .true.    

!--------------------------------------------------------------------
! creation du fichier de sortie et entete
!--------------------------------------------------------------------


  ierr = NF_CREATE('out.nc', NF_CLOBBER, out_file_id)
  if (ierr /= NF_NOERR) then
    write(lunout,*)NF_STRERROR(ierr)
    stop
  endif
  ierr = NF_INQ_NATTS(in_file_id, ngatts)
  do i = 1, ngatts
    ierr = NF_INQ_ATTNAME(in_file_id, NF_GLOBAL, i, attname)
    ierr = NF_COPY_ATT(in_file_id, NF_GLOBAL, attname, out_file_id, NF_GLOBAL)
  enddo
!
! definition des dimensions (nbin, level, temps)
! nbin:
  ierrs = 0
  ierr = NF_DEF_DIM(out_file_id, 'nbin', nb_bin, nbin_id)
  ierrs = ierrs + ierr
  nvdims = 1
  vdims(1) = nbin_id
  ierr = NF_DEF_VAR(out_file_id, 'nbin', NF_FLOAT, nvdims, vdims, var_bin_id)
  ierrs = ierrs + ierr
  ierr = NF_PUT_ATT_TEXT(out_file_id, var_bin_id, 'units', 5, 'hPa/j')
  ierrs = ierrs + ierr
  ierr = NF_PUT_ATT_TEXT(out_file_id, var_bin_id, 'long_name', 15, &
 &       'number of bins')
  ierrs = ierrs + ierr
  if (ierrs /= 4 * NF_NOERR) THEN
    write(lunout,*)'Pb. dans la definition de nbin'
    stop
  endif
!
! y bidon pour grads
  ierr = NF_DEF_DIM(out_file_id, 'lat', 1, latid)
  nvdims = 1
  vdims(1) = latid
  ierr = NF_DEF_VAR(out_file_id, 'lat', NF_FLOAT, nvdims, vdims,latid)
  ierr = NF_PUT_ATT_TEXT(out_file_id, latid, 'units', 13, 'degrees_north')

!
! presnivs
  ierrs = 0
  ierr = NF_DEF_DIM(out_file_id, 'presnivs', lm, outlev_id)
  ierrs = ierrs + ierr
  nvdims = 1
  vdims(1) = outlev_id
  ierr = NF_DEF_VAR(out_file_id,'presnivs',NF_FLOAT, nvdims, vdims, var_lev_id)
  ierrs = ierrs + ierr
  ierr = NF_INQ_VARID(in_file_id, 'presnivs', levid)
  ierrs = ierrs + ierr
  ierr = NF_INQ_VARNATTS(in_file_id, levid, natts)
  ierrs = ierrs + ierr
  do i = 1, natts
    attname = ''
    ierr = NF_INQ_ATTNAME(in_file_id, levid, i, attname)
    ierrs = ierrs + ierr
    ierr = NF_COPY_ATT(in_file_id, levid, attname, out_file_id, var_lev_id) 
    ierrs = ierrs + ierr
  enddo
  IF (ierrs /= (4 + natts) * NF_NOERR) THEN
    write(lunout,*)'Pb. dans la definition de presnivs'
    stop
  endif
  
!
! temps:
  ierrs = 0
  ierr = NF_DEF_DIM(out_file_id, 'time_counter', NF_UNLIMITED, time_id)
  ierrs = ierrs + ierr
  nvdims = 1
  vdims(1) = time_id
  ierr = NF_DEF_VAR(out_file_id, 'time_counter', NF_FLOAT, nvdims, vdims, &
 &                  var_time_id)
  ierrs = ierrs + ierr
  ierr = NF_INQ_VARID(in_file_id, 'time_counter', in_var_time_id)
  if (ierr /= NF_NOERR) &
 &   ierr = NF_INQ_VARID(in_file_id, 't_ave_02592000', in_var_time_id)
  if (ierr == NF_NOERR) then
    ierr = NF_COPY_ATT(in_file_id, in_var_time_id, 'units', out_file_id, &
 &                     var_time_id)
  else
    ierr = NF_PUT_ATT_TEXT(out_file_id, var_time_id, 'units', 34, &
 &                         'seconds since 1860-01-01 00:00:00')
    in_var_time_id = 0
  endif
  ierrs = ierrs + ierr
  ierr = NF_PUT_ATT_TEXT(out_file_id, var_time_id, 'long_name', 9, &
 &       'Time axis')
  ierrs = ierrs + ierr
  if (ierrs /= 4 * NF_NOERR) THEN
    write(lunout,*)'Pb. dans la definition de time_counter'
    stop
  endif
!
! Definition des variables a ecrire:
! pdf:
  ierrs = 0
  nvdims = 2
  vdims(1) = nbin_id
  vdims(2) = time_id
  ierr = NF_DEF_VAR(out_file_id, 'pdf', NF_FLOAT, nvdims, vdims, & 
  &                 var_pdf_id)
  ierrs = ierrs + ierr
  ierr = NF_PUT_ATT_TEXT(out_file_id, var_pdf_id, 'long_name', 16, &
 &       'monthly PDF w500')
  ierrs = ierrs + ierr
  ierr = NF_PUT_ATT_REAL(out_file_id, var_pdf_id, 'missing_value', NF_FLOAT, &
 &       1, undef)
  ierrs = ierrs + ierr
  if (ierrs /= 3 * NF_NOERR) THEN
    write(lunout,*)'Pb. dans la definition de pdf'
    stop
  endif
! nx:
  ierrs = 0
  nvdims = 2
  vdims(1) = nbin_id
  vdims(2) = time_id
  ierr = NF_DEF_VAR(out_file_id, 'nx', NF_FLOAT, nvdims, vdims, & 
  &                 var_nx_id)
  ierrs = ierrs + ierr
  ierr = NF_PUT_ATT_TEXT(out_file_id, var_nx_id, 'long_name', 16, &
 &       'nb of points in w500 bin')
  ierrs = ierrs + ierr
  if (ierrs /= 2 * NF_NOERR) THEN
    write(lunout,*)'Pb. dans la definition de nx'
    stop
  endif
!
! w:
  ierrs = 0
  nvdims = 2
  vdims(1) = nbin_id
  vdims(2) = time_id
  ierr = NF_DEF_VAR(out_file_id, 'w', NF_FLOAT, nvdims, vdims, & 
  &                 var_w_id)
  ierrs = ierrs + ierr
  long_name = 'mean relationship w(binw)'
  ierr = NF_PUT_ATT_TEXT(out_file_id, var_w_id, 'long_name', 25, &
 &       long_name)
  ierrs = ierrs + ierr
  ierr = NF_PUT_ATT_REAL(out_file_id, var_w_id, 'missing_value', NF_FLOAT, &
 &       1, undef)
  ierrs = ierrs + ierr
  if (ierrs /= 3 * NF_NOERR) THEN
    write(lunout,*)'Pb. dans la definition de w'
    stop
  endif
!
! x:
  ierrs = 0
  nvdims = 2
  vdims(1) = nbin_id
  vdims(2) = time_id
  if (var_3d) then
    nvdims = 3
    vdims(2) = outlev_id
    vdims(3) = time_id
  endif
  ierr = NF_DEF_VAR(out_file_id, varname, NF_FLOAT, nvdims, vdims, & 
  &                 var_x_id)
  ierrs = ierrs + ierr
  long_name = 'mean relationship '//trim(varname)
  ierr = NF_PUT_ATT_TEXT(out_file_id, var_x_id, 'long_name', len(long_name), &
 &       long_name)
  ierrs = ierrs + ierr
  ierr = NF_PUT_ATT_REAL(out_file_id, var_x_id, 'missing_value', NF_FLOAT, &
 &       1, undef)
  ierrs = ierrs + ierr
  if (ierrs /= 3 * NF_NOERR) THEN
    write(lunout,*)'Pb. dans la definition de ',varname
    stop
  endif
  ierr = NF_ENDDEF(out_file_id)

!
! On passe aux calculs
!
! initialisations:

  allocate(xpdfmean(nb_bin))
  allocate(xpdf(nb_bin, itime))
  allocate(nx_binw(nb_bin, itime))
  allocate(sx_binw(nb_bin, itime))
  allocate(w_binw(nb_bin, itime))
  allocate(x_binw(nb_bin, itime))

!-- temporal loop:
  il = 1
  if (var_3d) then
    il = lm
  endif

  do l = 1, il    
  xpdfmeantot = 0.
  xpdfmean = 0.0
  x_binw = 0.0
  xpdf = 0.0
  nx_binw = 0.0
  sx_binw = 0.0
  w_binw = 0.0
  xpdftot = 0.
  do m = 1, itime ! loop over timeperiod of x_data

! read data:
    start(1) = 1
    start(2) = 1
    start(3) = m
    count(1) = im
    count(2) = jm
    count(3) = 1
! vit. verticale à 500hP
    if (.not. allocated(w)) allocate (w(im,jm))
    ierrs = 0
    ierr = NF_INQ_VARID(in_file_id, 'w500', varid)
    ierrs = ierrs + ierr
    ierr = NF_GET_VARA_REAL(in_file_id, varid, start, count, w)
    ierrs = ierrs + ierr
    if (ierrs /= 2* NF_NOERR) THEN
      write(lunout,*)'Pb. avec la lecture de la vit. vert. à 500hP'
      stop
    endif
! 
! variable a traiter
    if (.not. allocated(x)) allocate (x(im,jm))
    if (var_3d) then
      start(3) = l
      count(3) = 1
      start(4) = m
      count(4) = 1
    endif
    ierrs = 0
    ierr = NF_INQ_VARID(in_file_id, varname, varid)
    ierrs = ierrs + ierr
    ierr = NF_GET_VARA_REAL(in_file_id, varid, start, count, x)
    ierrs = ierrs + ierr
    if (ierrs /= 2 * NF_NOERR) THEN
      write(lunout,*)'Pb. avec la lecture de ', varname
      stop
    endif
!
    do j=1,jm
      do i=1,im
        zmsk(i,j)=0.
      enddo
    enddo
    nbr_inf = 0; nbr_sup = 0 
! tropical belt:
    do j = 1, jm
      if (ABS(lat(j)).le.lat0) then
! loop over grid points:
        do i = 1, im
          if (x(i,j).ne.undef_x .and. w(i,j).ne.undef_w     &
  &           .and. ok_msk(i,j) ) then
            zmsk(i,j)=1.
            w1 = w(i,j)*w_mult
            x1 = x(i,j)*x_mult
! bin w500:
            if (w1 < min_bin) then
              ir = 1
              nbr_inf = nbr_inf + 1
            else if (w1 > max_bin) then
              ir = nb_bin
              nbr_sup = nbr_sup + 1
            else
              ir = INT((w1-min_bin)/step_bin) + 1
            endif
! monthly PDF :
            xpdfmeantot = xpdfmeantot + surf(i,j)
            xpdfmean(ir) = xpdfmean(ir) + surf(i,j)
            xpdftot(m) = xpdftot(m) + surf(i,j)
            xpdf(ir,m) = xpdf(ir,m) + surf(i,j)
! monthly x-w relationship:
            nx_binw(ir,m) = nx_binw(ir,m) + 1.0
            sx_binw(ir,m) = sx_binw(ir,m) + surf(i,j)
            x_binw(ir,m) = x_binw(ir,m) + x1*surf(i,j)
            w_binw(ir,m) = w_binw(ir,m) + w1*surf(i,j)

          endif
        enddo ! i
      endif ! lat
    enddo ! j
    if (nbr_inf /=0) &
    & write(lunout,*)'nbre de points ou w500 < ',min_bin,' = ', nbr_inf
    if (nbr_sup /=0) &
    & write(lunout,*)'nbre de points ou w500 > ',max_bin,' = ', nbr_sup 
  enddo ! m

! normalize PDF:

  do ir = 1, nb_bin
    xpdfmean(ir) = xpdfmean(ir)/xpdfmeantot
  enddo

  do m = 1, itime
    if (xpdftot(m).gt.0.) then
      do ir = 1, nb_bin
        xpdf(ir,m) = xpdf(ir,m)/xpdftot(m)
      enddo
    else
      do ir = 1, nb_bin
        xpdf(ir,m) = undef
      enddo
    endif
  enddo


  do ir = 1, nb_bin
    do m = 1, itime
      if (nx_binw(ir,m).gt.1.) then
!      if (nx_binw(ir,m).gt.0.) then
        x_binw(ir,m) = x_binw(ir,m)/sx_binw(ir,m)
        w_binw(ir,m) = w_binw(ir,m)/sx_binw(ir,m)
      else
        x_binw(ir,m) = undef
        w_binw(ir,m) = undef
      endif
    enddo
  enddo ! ir

!
! ecriture du fichier sortie
!

  allocate(var_dim(nb_bin))
  do ir = 1, nb_bin
    var_dim(ir) = min_bin + step_bin * (ir -1)
  enddo
  ierr = NF_PUT_VAR_REAL(out_file_id, var_bin_id, var_dim)
  if (ierr /= NF_NOERR) then
    write(lunout,*)NF_STRERROR(ierr)
    stop
  endif
  deallocate(var_dim)
!
! ecriture y bidon
!
  ierr = NF_PUT_VAR1_REAL(out_file_id, latid, 1, 0.)
  
!
! ecriture niveaux verticaux
!
  allocate(var_dim(lm))
  ierrs = 0
  ierr = NF_GET_VAR_REAL(in_file_id, levid, var_dim)
  ierrs = ierrs + ierr
  ierr = NF_PUT_VAR_REAL(out_file_id, var_lev_id, var_dim)
  ierrs = ierrs + ierr
  if (ierrs /= 2 * NF_NOERR) THEN
    write(lunout,*)'Pb. d''ecriture niveaux verticaux'
    stop
  endif
  deallocate(var_dim)
  
  
    
  do m = 1, itime
    start(1) = 1
    start(2) = m
    count(1) = nb_bin
    count(2) = 1
    ierrs = 0
!
! ecriture temps:
    if (in_var_time_id == 0) then
      time = (m - 1) * 86400. * 30.
    else
      ierr = NF_GET_VAR1_REAL(in_file_id, in_var_time_id, m, time)
    endif
    ierr = NF_PUT_VAR1_REAL(out_file_id, var_time_id, m, time)
    ierrs = ierrs + ierr
!
! pdf, nb_bin, x, w:
    ierr = NF_PUT_VARA_REAL(out_file_id, var_pdf_id, start, count, &
   &                        xpdf(1,m))
    ierrs = ierrs + ierr
    ierr = NF_PUT_VARA_REAL(out_file_id, var_nx_id, start, count, &
   &                        nx_binw(1,m))
    ierrs = ierrs + ierr
    ierr = NF_PUT_VARA_REAL(out_file_id, var_w_id, start, count, &
   &                        w_binw(1,m))
    ierrs = ierrs + ierr
    if (var_3d) then
      start(2) = l
      count(2) = 1
      start(3) = m
      count(3) = 1
    endif
    ierr = NF_PUT_VARA_REAL(out_file_id, var_x_id, start, count, &
   &                        x_binw(1,m))
    ierrs = ierrs + ierr
  if (ierrs /= 5 * NF_NOERR) THEN
    write(lunout,*)'Pb. d''ecriture'
    stop
  endif
  enddo ! m
  enddo ! l

  ierr = NF_CLOSE(out_file_id)

  end