 MODULE outnetcdf
   !!======================================================================
   !!                   ***  MODULE  outnetcdf  ***
   !! Outputs : tools to complete ionc4
   !!
   !!   History :
   !!        !  2017-07 (V. Garnier)
   !!
   !!======================================================================

#include "toolcpp.h"

   !! * Modules used

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC out_netcdf_conv,       &    ! routine called by output_save.F90
          out_zaxis_w

   !! * Shared or public module variables

   !! * Private variables

 CONTAINS

   !!======================================================================

  SUBROUTINE out_netcdf_conv(filename,nb,comm_zoom,ioutmin,ioutmax,joutmin,joutmax)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE out_netcdf_conv  ***
   !&E
   !&E ** Purpose : add attributes to nedCDF file in order to follow
   !&E              NetCDF OCO convention.
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : output_save.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays :
   !&E
   !&E ** Modified variables :
   !&E
   !&E ** Reference : Stephane Tarot (Csv $ NetCDF OCO user s manual,
   !&E                09/12/2008 - COR-DO/DTI-MUT/07.159
   !&E
   !&E ** History :
   !&E       !  2008-12  (V. Garnier)  Original code
   !&E       !  2009-03 (V. Garnier) attributes according to NetCDF OCO v1.2 file format
   !&E       !  2011-01 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters, ONLY : lchain,rsh,rlg,riosh,riolg,rg_valmanq_io,           &
                          dg_valmanq_io,valmanq,kmax,imax,jmax,               &
                          limin,limax,ljmin,ljmax,imin,jmin,jshift_obc_south, &
                          ishift_obc_west,ishift_obc_east,jshift_obc_north
   USE comvars2d, ONLY  : cg_conv_institution,cg_conv_institution_references, &
                          cg_conv_data_centre,cg_conv_data_centre_references, &
                          cg_conv_title,cg_conv_comment,cg_conv_source,       &
                          cg_conv_forecast_range,cg_conv_forecast_type,       &
                          cg_conv_field_type,                                 &
                          cg_conv_operational_status,cg_conv_area,dfiac,dgac, &
                          h0,hx,hy,h0fond,lon2d_g,lat2d_g,name_out_xe,        &
                          name_out_h0,ig,id,jb,jh,l_obc_cycl_x,l_obc_cycl_y
   USE output, ONLY     : tdebsor,tfinsor,pasor,iextmin
   USE ionc4, ONLY      : ionc4_gatt_conv,ionc4_gatt,ionc4_vatt_conventions,  &
                          ionc4_vatt_valid_min,ionc4_vatt_valid_max,ionc4_vatt_fill,&
                          ionc4_vatt_positive,ionc4_vatt_chain
   USE obccombine, ONLY : l_obc_east,l_obc_west,l_obc_north,l_obc_south
   USE_MPI toolmpi,ONLY : min_all_mpi_rsh,max_all_mpi_rsh,MASTER,MPI_COMM_MARS

   IMPLICIT NONE

   !! * Arguments
   CHARACTER(LEN=lchain), INTENT( in ) :: filename
   INTEGER              , INTENT( in ) :: nb,ioutmin,ioutmax,joutmin,joutmax
   INTEGER              , INTENT( in ) :: comm_zoom

   !! * Local declarations
#ifdef gfortran
   INTEGER                            :: stime,communicator
#else
   INTEGER                            :: stime,time,communicator
#endif
   CHARACTER(LEN=19)                  :: comdate,tool_sectodat
   CHARACTER(LEN=lchain)              :: chain,lim_t,lim_u,lim_v,lim_f, &
                                         limi_t,limi_u,limi_v,limi_f,lima_t,lima_u,lima_v,lima_f
   REAL(KIND=rlg)                     :: torigin,tool_datosec
   REAL(KIND=riosh)                   :: xy
   REAL(KIND=rsh)                     :: ww
   REAL(KIND=rsh),DIMENSION(limin:limax,ljmin:ljmax) :: tmp2d_rsh

   !!----------------------------------------------------------------------
   !! * Executable part

     PRINT_DBG*, 'enter OUT_NETCDF_CONV routine',TRIM(filename)

     CALL ionc4_gatt_conv(filename)
     CALL ionc4_gatt(filename,'institution',TRIM(cg_conv_institution))
     CALL ionc4_gatt(filename,'institution_references',TRIM(cg_conv_institution_references))
     CALL ionc4_gatt(filename,'data_centre',TRIM(cg_conv_data_centre))
     CALL ionc4_gatt(filename,'data_centre_references',TRIM(cg_conv_data_centre_references))
     CALL ionc4_gatt(filename,'title',TRIM(cg_conv_title))

#if ! defined key_castest

     CALL getenv("PWD",chain)
     chain=TRIM(ADJUSTL(ADJUSTR(chain)))//'/'//TRIM(filename)
     CALL ionc4_gatt(filename,'product_name',TRIM(chain))

     stime = time()
!     print*,'Seconds since 01/01/01 00:00:00 = ', stime
     torigin=tool_datosec('01/01/1970 00:00:00')-tool_datosec('01/01/1900 00:00:00')
     comdate=tool_sectodat(REAL(stime+torigin,rlg))
     chain=comdate(7:10)//'-'//comdate(4:5)//'-'//comdate(1:2)
     chain=TRIM(chain)//'T'//comdate(12:19)//'Z'
     CALL ionc4_gatt(filename,'creation_date',TRIM(chain))
     CALL ionc4_gatt(filename,'run_time',TRIM(chain))
     chain=TRIM(chain)//': creation'
     CALL ionc4_gatt(filename,'history',TRIM(chain))

#endif

#ifdef key_MPI_2D
     communicator=comm_zoom
#endif

     CALL ionc4_gatt(filename,'model_name','MARS')
     CALL ionc4_gatt(filename,'arakawa_grid_type','C1')
     CALL ionc4_gatt(filename,'comment',TRIM(cg_conv_comment))
     CALL ionc4_gatt(filename,'source',TRIM(cg_conv_source))
     CALL ionc4_gatt(filename,'area',TRIM(cg_conv_area))
     !xy=lat2d_g(imin,jmin)
     xy=MINVAL(lat2d_g,lat2d_g>-valmanq)
     write(chain,123) xy
     CALL ionc4_gatt(filename,'southernmost_latitude',TRIM(ADJUSTL(chain)))
     !xy=lat2d_g(imax,jmax)
     xy=MAXVAL(lat2d_g,lat2d_g>-valmanq)
     write(chain,123) xy
     CALL ionc4_gatt(filename,'northernmost_latitude',TRIM(ADJUSTL(chain)))
     write(chain,*) dfiac
     CALL ionc4_gatt(filename,'latitude_resolution',TRIM(ADJUSTL(chain)))
     !xy=lon2d_g(imin,jmin)
     xy=MINVAL(lon2d_g,lon2d_g>-valmanq)
     write(chain,123) xy
     CALL ionc4_gatt(filename,'westernmost_longitude',TRIM(ADJUSTL(chain)))
     !xy=lon2d_g(imax,jmax)
     xy=MAXVAL(lon2d_g,lon2d_g>-valmanq)
     write(chain,123) xy
     CALL ionc4_gatt(filename,'easternmost_longitude',TRIM(ADJUSTL(chain)))
     write(chain,*) dgac
     CALL ionc4_gatt(filename,'longitude_resolution',TRIM(ADJUSTL(chain)))
     ! minimum and maximum depth
     tmp2d_rsh(limin:limax,ljmin:ljmax)=h0(limin:limax,ljmin:ljmax)
     WHERE (tmp2d_rsh(:,:)==-valmanq) tmp2d_rsh(:,:)=valmanq
     ww=MINVAL(tmp2d_rsh(:,:)) 
     CALL_MPI MIN_ALL_MPI_RSH(ww,comm_active=communicator)
     tmp2d_rsh(limin:limax,ljmin:ljmax)=hx(limin:limax,ljmin:ljmax)+h0fond
     WHERE (tmp2d_rsh(:,:)==-valmanq) tmp2d_rsh(:,:)=valmanq
     ww=MIN(MINVAL(tmp2d_rsh(:,:)),ww)
     CALL_MPI MIN_ALL_MPI_RSH(ww,comm_active=communicator)
     tmp2d_rsh(limin:limax,ljmin:ljmax)=hy(limin:limax,ljmin:ljmax)+h0fond
     WHERE (tmp2d_rsh(:,:)==-valmanq) tmp2d_rsh(:,:)=valmanq
     ww=MIN(MINVAL(tmp2d_rsh(:,:)),ww)
     CALL_MPI MIN_ALL_MPI_RSH(ww,comm_active=communicator)
     xy=REAL(ww,riosh)
     write(chain,*) xy
     CALL ionc4_gatt(filename,'minimum_depth',TRIM(ADJUSTL(chain)))
     ww=MAX(MAXVAL(h0(:,:)),MAXVAL(hx(:,:)),MAXVAL(hy(:,:)))
     CALL_MPI MAX_ALL_MPI_RSH(ww,comm_active=communicator)
     xy=REAL(ww,riosh)
     write(chain,*) xy
     CALL ionc4_gatt(filename,'maximum_depth',TRIM(ADJUSTL(chain)))
     CALL ionc4_gatt(filename,'depth_resolution',"n/a")

123 FORMAT(f8.4)

!     CALL ionc4_gatt(filename,'',)

     comdate=tool_sectodat(tdebsor(nb))
     chain=comdate(7:10)//'-'//comdate(4:5)//'-'//comdate(1:2)
     chain=TRIM(chain)//'T'//comdate(12:19)//'Z'
     CALL ionc4_gatt(filename,'start_date',TRIM(chain))
     comdate=tool_sectodat(tfinsor(nb))
     chain=comdate(7:10)//'-'//comdate(4:5)//'-'//comdate(1:2)
     chain=TRIM(chain)//'T'//comdate(12:19)//'Z'
     CALL ionc4_gatt(filename,'stop_date',TRIM(chain))
     IF ( TRIM(ADJUSTL(ADJUSTR(cg_conv_field_type)))=="n/a") THEN
       write(chain,*) pasor(nb)
       chain=TRIM(chain)//' s'
     ELSE
       chain=TRIM(ADJUSTL(ADJUSTR(cg_conv_field_type)))
     END IF
     CALL ionc4_gatt(filename,'field_type',TRIM(chain))

     CALL ionc4_gatt(filename,'forecast_range',TRIM(cg_conv_forecast_range))
     CALL ionc4_gatt(filename,'forecast_type',TRIM(cg_conv_forecast_type))
     CALL ionc4_gatt(filename,'operational_status',TRIM(cg_conv_operational_status))

!     CALL ionc4_vatt_fill(filename,'time',dg_valmanq_io)
     CALL ionc4_vatt_conventions(filename,'time','relative number of seconds with no decimal part')

!     CALL ionc4_vatt_fill(filename,'latitude',dg_valmanq_io)
     CALL ionc4_vatt_valid_min(filename,'latitude',-90.0_riolg)
     CALL ionc4_vatt_valid_max(filename,'latitude',90.0_riolg)
!     CALL ionc4_vatt_fill(filename,'longitude',dg_valmanq_io)
     CALL ionc4_vatt_valid_min(filename,'longitude',-180.0_riolg)
     CALL ionc4_vatt_valid_max(filename,'longitude',180.0_riolg)

     ! Horizontal axes : indexes of the variables estimated by the model (dynamical domain inside the obc)
     ! fortran convention : the fisrt point in the netcdf file is 1
     ! COMODO 1.0 convention
     ! ---------------------------------------------------------------------------------------------------
     IF (l_obc_cycl_x) THEN
       WRITE(limi_t,*) MAX((imin+ishift_obc_west+3)-ioutmin,1)
       WRITE(limi_u,*) MAX((imin+ishift_obc_west+3)-ioutmin,1)
       WRITE(limi_v,*) MAX((imin+ishift_obc_west+3)-ioutmin,1)
       WRITE(limi_f,*) MAX((imin+ishift_obc_west+3)-ioutmin,1)
     ELSE IF (l_obc_west) THEN
       WRITE(limi_t,*) MAX((imin+ishift_obc_west+3)-ioutmin,1)
       WRITE(limi_u,*) MAX((imin+ishift_obc_west+2)-ioutmin,1)
       WRITE(limi_v,*) MAX((imin+ishift_obc_west+3)-ioutmin,1)
       WRITE(limi_f,*) MAX((imin+ishift_obc_west+2)-ioutmin,1)
     ELSE
       WRITE(limi_t,*) MAX((MINVAL(ig)+2)-ioutmin,1)
       WRITE(limi_u,*) MAX((MINVAL(ig)+2)-ioutmin,1)
       WRITE(limi_v,*) MAX((MINVAL(ig)+2)-ioutmin,1)
       WRITE(limi_f,*) MAX((MINVAL(ig)+2)-ioutmin,1)
     END IF
     IF (l_obc_east .OR. l_obc_cycl_x) THEN
       WRITE(lima_t,*) MIN(ioutmax,imax-ishift_obc_east-1)-ioutmin+1
       WRITE(lima_u,*) MIN(ioutmax,imax-ishift_obc_east-1)-ioutmin+1
       WRITE(lima_v,*) MIN(ioutmax,imax-ishift_obc_east-1)-ioutmin+1
       WRITE(lima_f,*) MIN(ioutmax,imax-ishift_obc_east-1)-ioutmin+1
     ELSE
       ! a wall must have 2 cells at least
       WRITE(lima_t,*) MIN(ioutmax,MAXVAL(id)-2)-ioutmin+1
       WRITE(lima_u,*) MIN(ioutmax,MAXVAL(id)-3)-ioutmin+1
       WRITE(lima_v,*) MIN(ioutmax,MAXVAL(id)-2)-ioutmin+1
       WRITE(lima_f,*) MIN(ioutmax,MAXVAL(id)-3)-ioutmin+1
     END IF
     lim_t=TRIM(ADJUSTL(ADJUSTR(limi_t)))//':'//TRIM(ADJUSTL(ADJUSTR(lima_t)))
     lim_u=TRIM(ADJUSTL(ADJUSTR(limi_u)))//':'//TRIM(ADJUSTL(ADJUSTR(lima_u)))
     lim_v=TRIM(ADJUSTL(ADJUSTR(limi_v)))//':'//TRIM(ADJUSTL(ADJUSTR(lima_v)))
     lim_f=TRIM(ADJUSTL(ADJUSTR(limi_f)))//':'//TRIM(ADJUSTL(ADJUSTR(lima_f)))
     CALL ionc4_vatt_chain(filename,'ni','c_grid_dynamic_range',TRIM(ADJUSTL(ADJUSTR(lim_t)))) 
     CALL ionc4_vatt_chain(filename,'ni_u','c_grid_dynamic_range',TRIM(ADJUSTL(ADJUSTR(lim_u)))) 
     CALL ionc4_vatt_chain(filename,'ni_v','c_grid_dynamic_range',TRIM(ADJUSTL(ADJUSTR(lim_v)))) 
     CALL ionc4_vatt_chain(filename,'ni_f','c_grid_dynamic_range',TRIM(ADJUSTL(ADJUSTR(lim_f)))) 

     IF (l_obc_cycl_y) THEN
       WRITE(limi_t,*) MAX((jmin+jshift_obc_south+3)-joutmin,1)
       WRITE(limi_u,*) MAX((jmin+jshift_obc_south+3)-joutmin,1)
       WRITE(limi_v,*) MAX((jmin+jshift_obc_south+3)-joutmin,1)
       WRITE(limi_f,*) MAX((jmin+jshift_obc_south+3)-joutmin,1)
     ELSE IF (l_obc_south) THEN
       WRITE(limi_t,*) MAX((jmin+jshift_obc_south+3)-joutmin,1)
       WRITE(limi_u,*) MAX((jmin+jshift_obc_south+3)-joutmin,1)
       WRITE(limi_v,*) MAX((jmin+jshift_obc_south+2)-joutmin,1)
       WRITE(limi_f,*) MAX((jmin+jshift_obc_south+2)-joutmin,1)
     ELSE
       WRITE(limi_t,*) MAX((MINVAL(jb)+2)-joutmin,1)
       WRITE(limi_u,*) MAX((MINVAL(jb)+2)-joutmin,1)
       WRITE(limi_v,*) MAX((MINVAL(jb)+2)-joutmin,1)
       WRITE(limi_f,*) MAX((MINVAL(jb)+2)-joutmin,1)
     END IF
     IF (l_obc_north .OR. l_obc_cycl_y) THEN
       WRITE(lima_t,*) MIN(joutmax,jmax-jshift_obc_north-1)-joutmin+1
       WRITE(lima_u,*) MIN(joutmax,jmax-jshift_obc_north-1)-joutmin+1
       WRITE(lima_v,*) MIN(joutmax,jmax-jshift_obc_north-1)-joutmin+1
       WRITE(lima_f,*) MIN(joutmax,jmax-jshift_obc_north-1)-joutmin+1
     ELSE
       ! a wall must have 2 cells at least
       WRITE(lima_t,*) MIN(joutmax,MAXVAL(jh)-2)-joutmin+1
       WRITE(lima_u,*) MIN(joutmax,MAXVAL(jh)-2)-joutmin+1
       WRITE(lima_v,*) MIN(joutmax,MAXVAL(jh)-3)-joutmin+1
       WRITE(lima_f,*) MIN(joutmax,MAXVAL(jh)-3)-joutmin+1
     END IF
     lim_t=TRIM(ADJUSTL(ADJUSTR(limi_t)))//':'//TRIM(ADJUSTL(ADJUSTR(lima_t)))
     lim_u=TRIM(ADJUSTL(ADJUSTR(limi_u)))//':'//TRIM(ADJUSTL(ADJUSTR(lima_u)))
     lim_v=TRIM(ADJUSTL(ADJUSTR(limi_v)))//':'//TRIM(ADJUSTL(ADJUSTR(lima_v)))
     lim_f=TRIM(ADJUSTL(ADJUSTR(limi_f)))//':'//TRIM(ADJUSTL(ADJUSTR(lima_f)))
     CALL ionc4_vatt_chain(filename,'nj','c_grid_dynamic_range',TRIM(ADJUSTL(ADJUSTR(lim_t)))) 
     CALL ionc4_vatt_chain(filename,'nj_u','c_grid_dynamic_range',TRIM(ADJUSTL(ADJUSTR(lim_u)))) 
     CALL ionc4_vatt_chain(filename,'nj_v','c_grid_dynamic_range',TRIM(ADJUSTL(ADJUSTR(lim_v)))) 
     CALL ionc4_vatt_chain(filename,'nj_f','c_grid_dynamic_range',TRIM(ADJUSTL(ADJUSTR(lim_f)))) 

     ! Vertical axes
     ! -------------
!     CALL ionc4_vatt_fill(filename,'nk',rg_valmanq_io)
     CALL ionc4_vatt_positive(filename,'level','up')
     CALL ionc4_vatt_valid_min(filename,'level',-1.0_riosh)
     CALL ionc4_vatt_valid_max(filename,'level',0.0_riosh)

#ifdef key_siggen_g2
     chain='ocean_s_coordinate_g2'
     CALL ionc4_vatt_chain(filename,'level','standard_name',TRIM(chain))
     chain='s: level C: Csu_sig eta: '//TRIM(name_out_xe)//' depth: '//TRIM(name_out_h0)//' depth_c: hc'
     CALL ionc4_vatt_chain(filename,'level','formula_terms',TRIM(chain))
     CALL ionc4_vatt_chain(filename,'level','formula_definition', &
         'z(n,k,j,i) = eta(n,j,i) + (eta(n,j,i)+depth(j,i))*S(k,j,i) ; ' // &
         'S(k,j,i) = (depth_c*s(k) + depth(j,i)*C(k)) / (depth_c + depth(j,i))')
#elif key_siggen
     chain='ocean_s_coordinate'
     CALL ionc4_vatt_chain(filename,'level','standard_name',TRIM(chain))
     chain='s: level eta: '//TRIM(name_out_xe)//' depth: '//TRIM(name_out_h0)//' a: theta b: b depth_c: hc'
     CALL ionc4_vatt_chain(filename,'level','formula_terms',TRIM(chain))
     CALL ionc4_vatt_chain(filename,'level','formula_definition', &
          'z(n,k,j,i) = eta(n,j,i)*(1+s(k)) + depth_c*s(k) + (depth(j,i)-depth_c)*C(k)')
#else
     chain='ocean_sigma_coordinate'
     CALL ionc4_vatt_chain(filename,'level','standard_name',TRIM(chain))
     chain='sigma: level eta: '//TRIM(name_out_xe)//' depth: '//TRIM(name_out_h0)
     CALL ionc4_vatt_chain(filename,'level','formula_terms',TRIM(chain))
     CALL ionc4_vatt_chain(filename,'level','formula_definition', &
          'z(n,k,j,i) = eta(n,j,i) + sigma(k)*(depth(j,i)+eta(n,j,i))')
#endif

     PRINT_DBG*, 'exit OUT_NETCDF_CONV routine',TRIM(filename)

  END SUBROUTINE out_netcdf_conv

   !!===========================================================================

  SUBROUTINE out_zaxis_w(filename,lkmi,lkma,pask,comm_zoom)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE out_zaxis_w  ***
   !&E
   !&E ** Purpose : add attributes to nedCDF file in order to follow
   !&E              NetCDF OCO convention.
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : output_save.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays :
   !&E
   !&E ** Modified variables :
   !&E
   !&E ** Reference : Stephane Tarot (Csv $ NetCDF OCO user s manual,
   !&E                09/12/2008 - COR-DO/DTI-MUT/07.159
   !&E
   !&E ** History :
   !&E       !  2008-12  (V. Garnier)  Original code
   !&E       !  2009-03 (V. Garnier) attributes according to NetCDF OCO v1.2 file format
   !&E       !  2011-01 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE netcdf
   USE parameters, ONLY : lchain,rsh,rlg,riosh,riolg,rg_valmanq_io,           &
                          dg_valmanq_io,valmanq,kmax
   USE comionc4,  ONLY  : ionc_pnomz_w
   USE ionc4, ONLY      : ionc4_vatt_positive,ionc4_vatt_chain,ionc4_corres,  &
                          ionc4_vatt_valid_min,ionc4_vatt_valid_max
   USE comvars2d, ONLY  : name_out_xe,name_out_h0
   USE comvars3d, ONLY  : sigw
 
   USE_MPI toolmpi, ONLY : MPI_COMM_MARS

   IMPLICIT NONE

   !! * Arguments
   CHARACTER(LEN=lchain), INTENT( in ) :: filename
   INTEGER,               INTENT( in ) :: lkma,lkmi,pask
   INTEGER, OPTIONAL    , INTENT( in ) :: comm_zoom

   !! * Local declarations
   INTEGER               :: communicator,dim_k,dim_w_id,nc_id,nc_err,sig_w_id
   INTEGER               :: ik,indice
   CHARACTER(LEN=lchain) :: chain
   REAL(KIND=4),DIMENSION(:),  ALLOCATABLE :: rsig

   !!----------------------------------------------------------------------
   !! * Executable part

     PRINT_DBG*, 'enter OUT_ZAXIS_W routine',TRIM(filename)

#ifdef key_MPI_2D
     IF ( PRESENT(comm_zoom) ) THEN
       communicator=comm_zoom
     ELSE
       communicator=MPI_COMM_MARS
     END IF
#endif

     dim_k   = (lkma-lkmi)/pask +1
     IF (dim_k ==  kmax) THEN
       dim_k = dim_k + 1
       ALLOCATE(rsig(0:dim_k))
       rsig(0) = sigw(0)
     ELSE
       ALLOCATE(rsig(dim_k))
     END IF
     indice = 1
     DO ik = lkmi,lkma,pask
       rsig(indice) = sigw(ik)
       indice = indice + 1
     END DO

     CALL ionc4_corres(TRIM(filename),nc_id)
     IF (nc_id == 0) THEN
       WRITE(*,*) 'Pb to open ',TRIM(filename),' and write zaxis at w location'
     ELSE

       nc_err = nf90_redef(nc_id) 
       IF (nc_err /= 0) WRITE(*,*) 'Pb to open for redefinition'

       ! create the dimension
       IF (dim_k ==  kmax) THEN
         nc_err = nf90_def_dim(nc_id,TRIM(ionc_pnomz_w),dim_k+1,dim_w_id)
       ELSE
         nc_err = nf90_def_dim(nc_id,TRIM(ionc_pnomz_w),dim_k,dim_w_id)
       END IF
       IF (nc_err /= 0) WRITE(*,*) 'Pb to create the dimension',TRIM(ionc_pnomz_w)
       ! create the variable-dimension
       nc_err = nf90_def_var(nc_id, TRIM(ionc_pnomz_w), nf90_real, dim_w_id, sig_w_id)
       IF (nc_err /= 0) WRITE(*,*) 'Pb to create the variable',TRIM(ionc_pnomz_w)

       nc_err = nf90_put_att(nc_id,sig_w_id, 'long_name', 'sigma level at the interface')
       nc_err = nf90_put_att(nc_id,sig_w_id,'axis','Z')
       nc_err = nf90_put_att(nc_id,sig_w_id,'c_grid_axis_shift',-0.5)

       nc_err = nf90_enddef(nc_id) 

       CALL ionc4_vatt_positive(filename,'level_w','up')
       CALL ionc4_vatt_valid_min(filename,'level_w',-1.1_riosh)
       CALL ionc4_vatt_valid_max(filename,'level_w',0.1_riosh)

#ifdef key_siggen_g2
     chain='ocean_s_coordinate_g2_at_w_location'
     CALL ionc4_vatt_chain(filename,'level_w','standard_name',TRIM(chain))
     chain='s: level_w C: Csw_sig eta: '//TRIM(name_out_xe)//' depth: '//TRIM(name_out_h0)//' depth_c: hc'
     CALL ionc4_vatt_chain(filename,'level_w','formula_terms',TRIM(chain))
     CALL ionc4_vatt_chain(filename,'level_w','formula_definition', &
         'z(n,k,j,i) = eta(n,j,i) + (eta(n,j,i)+depth(j,i))*S(k,j,i) ; ' // &
         'S(k,j,i) = (depth_c*s(k) + depth(j,i)*C(k)) / (depth_c + depth(j,i))')
#elif key_siggen
       chain='ocean_s_coordinate_at_w_location'
       CALL ionc4_vatt_chain(filename,'level_w','standard_name',TRIM(chain))
       chain='s: level_w eta: '//TRIM(name_out_xe)//' depth: '//TRIM(name_out_h0)//' a: theta b: b depth_c: hc'
       CALL ionc4_vatt_chain(filename,'level_w','formula_terms',TRIM(chain))
       CALL ionc4_vatt_chain(filename,'level_w','formula_definition', &
          'z(n,k,j,i) = eta(n,j,i)*(1+s_w(k)) + depth_c*s_w(k) + (depth(j,i)-depth_c)*C_w(k)')
#else
       chain='ocean_sigma_coordinate_at_w_location'
       CALL ionc4_vatt_chain(filename,'level_w','standard_name',TRIM(chain))
       chain='sigma: level_w eta: '//TRIM(name_out_xe)//' depth: '//TRIM(name_out_h0)
       CALL ionc4_vatt_chain(filename,'level_w','formula_terms',TRIM(chain))
       CALL ionc4_vatt_chain(filename,'level_w','formula_definition', &
          'z(n,k,j,i) = eta(n,j,i) + sigma_w(k)*(depth(j,i)+eta(n,j,i))')
#endif

       IF (dim_k == kmax+1) THEN
       ! strange, probably because netcdf starts arrays by index 1
         nc_err = nf90_put_var(nc_id,sig_w_id,rsig(0:kmax))
       ELSE
         nc_err = nf90_put_var(nc_id,sig_w_id,rsig)
       END IF
       IF (nc_err /= 0) WRITE(*,*) 'Pb to write the variable sig_w'


     END IF

     PRINT_DBG*, 'exit OUT_ZAXIS_W routine',TRIM(filename)

  END SUBROUTINE out_zaxis_w

   !!======================================================================

 END MODULE outnetcdf
