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

  SUBROUTINE output_save(nb,tect,u,uz,x,v,vz,sal,temp,        &
                         vis2d,vish_xe,vish_phi,nz,kz,boz,    &


                                                      cv_wat, &

                                                             hexm,heym,xem,kzm)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE output_save  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description :
!&E
!&E ** Called by : output_mng.F90, mpimod.F90
!&E
!&E ** External calls :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  11-2004
!&E       !  01-2006 (V. Garnier introduction of different output types)
!&E       !  2008-03 (C. Struski) introduction of parameter l_out_subs to save the substance in output files
!&E       !  2008-10 (V. Garnier) output managing relative to displacement of init_fromfile from step to main
!&E       !  2008-10 (M. Jouan, V. Garnier) save conservative (in flux meanning) variable uz and vz (for offline purpose)
!&E       !  2008-09 (G. Galibert / ALTRAN) attributes of output files that respect NetCDF OCO v1.1 file format
!&E       !  2009-03 (V. Garnier) attributes according to NetCDF OCO v1.2 file format
!&E                               saving of generalized sigma parameters
!&E       !  2009-05  (V. Garnier)  Save cumulative atmospherical fluxes
!&E       !  2009-09  (P. Garreau)  Save viscosity and diffusivity variables
!&E       !  2009-12  (S. Michel, HOCER)  Save buoyancy variable
!&E       !  2010-02  (G. Simon) Save waves statistics: Hs, direction, foc
!&E       !  2010-07 (R. Verney) Introduction of orbital velocity variable
!&E       !  2011-01 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
!&E       !  2012-02  (P. Garreau) Saving of physical vertical velocity (key_diag_physwz)
!&E       !  2014-02  (B.Thouvenin) Modification of sand output if 1
!&E       !  2014-05  (R. Verney) Introduction of module flocculation (key_flocmod)
!&E       !  2017-09  (F. Ganthy) Turbulence advection (k-epsilon)
!&E
!&E---------------------------------------------------------------------





























































!! Modules used
   USE netcdf
   USE ionc4,       ONLY : ionc4_createfile,ionc4_write_lon,ionc4_gatt,ionc4_createvar,&
                           ionc4_sync,ionc4_write,ionc4_write_sig,ionc4_write_sigw,&
                           ionc4_write_xy,ionc4_write_zxy,ionc4_write_time,&
                           ionc4_write_xyt,ionc4_write_zxyt,ionc4_close,&
                           ionc4_vatt_chain
   USE parameters,  ONLY : rsh,rlg,riosh,kmax,imin,imax,limin,limax,ljmin,ljmax,  &
                           rg_valmanq_io,ljminp1,riolg,valmanq,lchain,grav,pi
   USE output,      ONLY : fileresu,file_output,irlos_cum,latlos_cum,senslos_cum, &
                           rad_cum,evap_cum,rain_cum,tosx_cum,tosy_cum,dt_cum,dtn,&

                           l_filebydate,tsort,tfinsor,stepi,stepj,stepk,dt_ave,   &
                           iextmin,iextmax,jextmin,jextmax,kextmin,kextmax,nbave, &
                           igextmin,igextmax,jgextmin,jgextmax,nbatm,l_out_zaxis_w,&
                           l_out_ssh,l_out_sal,l_out_temp,l_out_ubt,l_out_u3d,    &
                           l_out_visc,l_out_turb,l_out_buoy,l_out_atmsflx,        &
                           l_out_timeave,fileresu_verif,                          &
                           l_out_wind,l_out_off,l_out_wave,l_out_nc4par,l_out_pack
   USE output,      ONLY : ssh_ave,vis2d_ave,u_ave,v_ave,uz_ave,vz_ave, &
                           temp_ave,sal_ave,vish_xe_ave,vish_phi_ave,   &
                           nz_ave,kz_ave,boz_ave


   USE outnetcdf,   ONLY : out_netcdf_conv,out_zaxis_w
   USE outtools,    ONLY : out_write_xytvar,out_write_zxytvar,out_grid
   USE comvars2d,   ONLY : lon2d,lon2d_g,lat2d_g,lat2d,fu,fv
   USE comvars2d,   ONLY : nbouc,date_ref,h0fond,date_startshift,iwarnlog,z0b,          &
                           h0,hx,hy,name_out_hx,name_out_hy,name_out_h0,name_out_vz,name_out_sal,name_out_temp, &
                           name_out_vis2d,name_out_vishxe,name_out_vishphi,name_out_kz,name_out_nz,name_out_bz, &
                           name_out_xe,name_out_u,name_out_v,name_out_uz,l_modele2d,l_initfromfile,             &
                           ig,id,jb,jh,name_out_tosx,name_out_tosy,out_dir,suffix,iscreen,                      &
                           name_out_mhex,name_out_mhey,name_out_mxe,name_out_mkz,                               &
                           name_out_mrad,name_out_mirlos,name_out_mlatlos,name_out_msenslos,                    &
                           name_out_mevap,name_out_mrain,name_out_mtosx,name_out_mtosy,                         &
                           riog_valid_min_h0,riog_valid_max_h0,riog_valid_min_u,riog_valid_max_u,               &
                           riog_valid_min_temp,riog_valid_max_temp,riog_valid_min_sal,riog_valid_max_sal,       &
                           riog_valid_min_xe,riog_valid_max_xe,riog_valid_min_wstress,riog_valid_max_wstress
   USE comvars2d,   ONLY : riog_valid_min_flux,riog_valid_max_flux,riog_valid_min_kz,riog_valid_max_kz,t,       &
                           riog_valid_min_visc,riog_valid_max_visc,riog_valid_min_nz,riog_valid_max_nz,         &
                           riog_valid_min_bz,riog_valid_max_bz,riog_valid_min_diag,riog_valid_max_diag
!MD 08/06/2018
   USE output,      ONLY : l_out_ws
   USE output,      ONLY : l_out_varint
!END MD 08/06/2018




   USE wave,        ONLY : name_out_wave_dir,riog_valid_min_wave_dir,riog_valid_max_wave_dir, &
                           name_out_wave_hs,riog_valid_min_wave_hs,riog_valid_max_wave_hs,    &
                           name_out_wave_foc,riog_valid_min_wave_foc,riog_valid_max_wave_foc, &
                           name_out_wave_ubot,riog_valid_min_wave_ubot,riog_valid_max_wave_ubot, &
                           name_out_wave_ubotx,name_out_wave_uboty, &
                           name_out_wave_t0m1,riog_valid_min_wave_t0m,riog_valid_max_wave_t0m,  &
                           name_out_wave_tauox,name_out_wave_tauoy,name_out_wave_tauwx,         &
                           name_out_wave_tauwy,riog_valid_min_wave_two,riog_valid_max_wave_two, &
                           name_out_wave_tusx,name_out_wave_tusy,name_out_wave_ussx,name_out_wave_ussy, &
                           riog_valid_min_wave_stokes,riog_valid_max_wave_stokes, &
                           name_out_wave_wlm,riog_valid_min_wave_wlm,riog_valid_max_wave_wlm, &
                           name_out_wave_bhd,riog_valid_min_wave_bhd,riog_valid_max_wave_bhd, &
                           name_out_wave_usx,name_out_wave_usy,name_out_wave_usz,             &
                           name_out_wave_tbbx,name_out_wave_tbby,name_out_wave_fbb,           &
                           name_out_wave_sx,name_out_wave_sy

   USE comvars3d,   ONLY : sig,sigw


   USE cometdyn,    ONLY : tosx,tosy

   USE wave,        ONLY : wave_dir,wave_hs,wave_phioc,wave_ubot,l_wave_rdir,l_wave_rhs,l_wave_rfoc,l_wave_rubot, &
                           l_wave_ruubot, l_wave_rvubot,wave_tauox,wave_tauoy,wave_bhd,wave_wlm,wave_t0m1, &
                           l_wave_rfbb,l_wave_rtbb,l_wave_rtauwx,l_wave_rtauwy,wave_tauwx,wave_tauwy, &
                           wave_tusx,wave_tusy,wave_ussx,wave_ussy,l_wave_rtus,l_wave_russ,wave_uubot,wave_vubot, &
                           l_wave_rt0m1,l_wave_rwlm,l_wave_rtauox,l_wave_rtauoy,l_wave_rbhd,l_wave_rfbb,l_wave_rtbb,&
                           l_wave_rradstr,wave_sx,wave_sy


   USE wave,        ONLY : wave_km,wave_km_iter,wave_usz,wave_vsz,wave_wsz,wave_taubbl_x,wave_taubbl_y,wave_phibbl

   USE sflxatm,     ONLY : rad,irlos,latlos,senslos
   USE comfluxdo,   ONLY : evap,rain


   USE parameters,  ONLY : nb_var
   USE comsubstance,ONLY : irk_fil,name_var,unit_var,nv_state,l_out_log,l_out_subs,   &
                           long_name_var,standard_name_var,valid_min_var,valid_max_var, &
                           unit_modif_mudbio_N2dw,l_subs2D,                           &
!MD 08/06/2018
                           name_var_sf,long_name_var_sf,standard_name_var_sf,          &
                           imud1,imud2,ws3,                                             &
                           name_var_sf1,long_name_var_sf1,standard_name_var_sf1
!END MD 08/06/2018



   USE comMUSTANG,   ONLY : rouse2D,sum_tmp,l_outsandrouse
   USE comsubstance,ONLY : isand1,isand2
   USE toolgeom,    ONLY : f_zu,f_zw
   USE comvars2d,   ONLY : hm

!!

   USE obccombine,  ONLY : l_obc_tide
   USE toolmpi, ONLY : ADD_TAG,IPROC,NPROCS,ADD_ALL_MPI_INT, MPI_COMM_MARS



   IMPLICIT NONE


   INCLUDE 'mpif.h'


!! * Arguments
   INTEGER       ,INTENT(in) :: nb
   REAL(KIND=rlg),INTENT(in) :: tect
   REAL(KIND=rsh),DIMENSION(     limin:limax,ljmin:ljmax),INTENT(in)          :: vis2d
   REAL(KIND=rsh),DIMENSION(     limin:limax,ljmin:ljmax),INTENT(in)          :: u,x,v
   REAL(KIND=rsh),DIMENSION(kmax,limin:limax,ljmin:ljmax),INTENT(in)          :: uz,vz,sal,temp
   REAL(KIND=rsh),DIMENSION(kmax,limin:limax,ljmin:ljmax),INTENT(in)          :: vish_xe,vish_phi
   REAL(KIND=rsh),DIMENSION(kmax,limin:limax,ljmin:ljmax),INTENT(in)          :: boz
   REAL(KIND=rsh),DIMENSION(0:kmax,limin:limax,ljmin:ljmax),INTENT(in)        :: nz,kz
   REAL(KIND=rsh),DIMENSION(     limin:limax,ljmin:ljmax),INTENT(in),OPTIONAL :: hexm,heym,xem
   REAL(KIND=rsh),DIMENSION(0:kmax,limin:limax,ljmin:ljmax),INTENT(in),OPTIONAL :: kzm

   REAL(KIND=rsh),DIMENSION(nb_var,kmax,limin:limax,ljmin:ljmax),INTENT(in) :: cv_wat


!! * Local declarations
   INTEGER               :: i,j,k,l
   INTEGER               :: grp_world,grp_zoom,comm_zoom,ncpu,ip,ierr_mpi
   INTEGER, DIMENSION(:), ALLOCATABLE  :: numrank,listrank
   INTEGER,DIMENSION(1) :: listrank_nb
   CHARACTER(LEN=lchain) :: fileres,gridfile,name_out
   CHARACTER(LEN=lchain) :: chain
   CHARACTER (len=19)    :: tool_sectodat,cdateout
   LOGICAL,         DIMENSION(limin:limax,ljmin:ljmax)      :: mask_h0,mask_hx,mask_hy
   REAL(KIND=rsh)        :: xx,yy
   REAL(KIND=riosh)      :: vecrit
   REAL(KIND=riosh),DIMENSION(limin:limax,ljmin:ljmax)      :: tecrit
   REAL(KIND=riosh),DIMENSION(kmax,limin:limax,ljmin:ljmax) :: tzecrit
   REAL(KIND=rsh),DIMENSION(limin:limax,ljmin:ljmax)        :: tmp2d
   REAL(KIND=rsh),DIMENSION(kmax,limin:limax,ljmin:ljmax)   :: tmp3d
! toseelater virer variables tmp3du,tmp3dv,tmp3dw quand variables dz integrees globalement
   REAL(KIND=rsh),DIMENSION(kmax,limin:limax,ljmin:ljmax)   :: tmp3du,tmp3dv
   REAL(KIND=rsh),DIMENSION(0:kmax,limin:limax,ljmin:ljmax) :: tmp3dw
   REAL(KIND=riosh),DIMENSION(0:kmax+1)                     :: sigecrit
   REAL(KIND=riolg),DIMENSION(imin:imax)                    :: lon
!MD 08/06/2018
!   REAL(KIND=rsh),DIMENSION(kmax,limin:limax,ljmin:ljmax)       :: ws_mud
!END MD 08/06/2018

   REAL(KIND=rsh)        :: htott,zk
   LOGICAL               :: l_outsandrouse_l




   CHARACTER(LEN=lchain), PARAMETER :: name_out_km='KM'
   CHARACTER(LEN=lchain), PARAMETER :: name_out_km_iter='iter'

   CHARACTER(len=lchain) :: name_out_z0b="Z0B"

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

   !PRINT*, 'enter OUTPUT_SAVE routine'

! Masks
! -----

   mask_h0(:,:)=.TRUE.
   DO j=ljminp1,ljmax
     DO i=MAX(limin,ig(j)),MIN(limax,id(j))
! For heights smaller than h0fond, the cell is considered as wet
       IF ( (j >= jb(i) .AND. j <= jh(i)) .AND. (h0(i,j)+x(i,j) > h0fond) ) mask_h0(i,j)=.FALSE.
       IF (i==31 ) mask_h0(i,j)=.false. 
     END DO
   END DO

   mask_hx(:,:)=.TRUE.
   DO j=ljminp1,ljmax
     DO i=MAX(limin,ig(j)-1),MIN(limax,id(j))
       IF ( (j >= jb(i) .AND. j <= jh(i)) .AND. (hx(i,j) /= -valmanq -h0fond) ) mask_hx(i,j)=.FALSE.
     END DO
   END DO

   mask_hy(:,:)=.TRUE.
   DO j=ljmin,ljmax
     DO i=MAX(limin,ig(j)),MIN(limax,id(j))
       IF ( (j >= jb(i)-1 .AND. j <= jh(i)) .AND. (hy(i,j) /= -valmanq -h0fond) ) mask_hy(i,j)=.FALSE.
     END DO
   END DO

   comm_zoom=0

   comm_zoom=MPI_COMM_MARS
!IF ( ( (nbouc == 0 .AND. .NOT.l_filebydate(nb)) .OR. (l_filebydate(nb) .AND. tect >= tsort(nb)) ) .AND. l_out_nc4par ) THEN
   IF ( (nbouc == 0 .AND. .NOT.l_filebydate(nb)) .OR. (l_filebydate(nb) .AND. tect >= tsort(nb)) ) THEN
! what a code ! Not a clever code probably but not clear in my head at the moment
! Count the number of processors active over the area to be saved
     ncpu=0
     IF ( jextmin(nb) <= jextmax(nb) .AND. iextmin(nb) <= iextmax(nb) ) ncpu=ncpu+1
     CALL ADD_ALL_MPI_INT(ncpu)
! Save the number of active processors
     ALLOCATE (numrank(ncpu))
     ALLOCATE (listrank(nprocs))
     numrank(:)=0   ! must be initialized if the proc 0 is active on the saved zoom
     listrank(:)=0
! which cpu is active ?
     IF ( jextmin(nb) <= jextmax(nb) .AND. iextmin(nb) <= iextmax(nb) ) listrank(1)=iproc
! all cpus get the number of active cpus plus sleeping cpus (set to 0)

     listrank_nb(1) = listrank(1)
     CALL MPI_ALLGATHER(listrank_nb(1),1,MPI_INTEGER,listrank,1,MPI_INTEGER,MPI_COMM_MARS,ierr_mpi)

! all cpus keep the active cpus only
     IF (ncpu/=nprocs) THEN
       ip=0
       DO i=1,nprocs 
         IF ( listrank(i) /= 0 ) THEN
         ip=ip+1
         numrank(ip)=listrank(i)
         ENDIF
       END DO
! Group relative to the MPI_COMM_MARS
       CALL MPI_COMM_GROUP(MPI_COMM_MARS,grp_world,ierr_mpi)
! Create the group of active cpus
       CALL MPI_GROUP_INCL(grp_world,ncpu,numrank,grp_zoom,ierr_mpi)
! Create the group of active cpus
       CALL MPI_COMM_CREATE(MPI_COMM_MARS,grp_zoom,comm_zoom,ierr_mpi)
     ENDIF
   ENDIF



   IF ( jextmin(nb) <= jextmax(nb) .AND. iextmin(nb) <= iextmax(nb) ) THEN

!-----------------------
! create the ouput files
!-----------------------

   IF ( (nbouc == 0 .AND. .NOT.l_filebydate(nb)) .OR. (l_filebydate(nb) .AND. tect >= tsort(nb)) ) THEN
!   IF ( nbouc == 0 .OR. l_filebydate(nb) ) THEN

     IF (l_filebydate(nb)) THEN
       cdateout = tool_sectodat(tect)
       fileres=TRIM(out_dir)//TRIM(fileresu(nb))//'_'//TRIM(suffix)
       fileres=TRIM(fileres)//'_'//cdateout(7:10)//cdateout(4:5)//cdateout(1:2)//'T'
       fileres=TRIM(fileres)//cdateout(12:13)//cdateout(15:16)//cdateout(18:19)//'Z'
       fileres=TRIM(fileres)//'.nc'
       IF (.NOT. l_out_nc4par) CALL ADD_TAG(fileres,LEN(fileres))
     ELSE
       fileres=TRIM(fileresu(nb))
     END IF
     fileresu_verif(nb)=fileres

!     DO i=imin,imax   ! to tina, yes it is correct, at least up to now.
!       lon(i)=gwac+(i-1)*dgac  !+360.0
!     END DO
     sigecrit(:)=sig(:)
     IF (l_out_nc4par) THEN
! To allow the different cpus to write at the same time into a global
! file, set l_out_nc4par= .TRUE.
! The dimension of the file must be global (igextmin)
       CALL ionc4_createfile(fileres,lon2d_g(igextmin(nb):igextmax(nb):stepi(nb),jgextmin(nb):jgextmax(nb):stepj(nb)), &
                                     lat2d_g(igextmin(nb):igextmax(nb):stepi(nb),jgextmin(nb):jgextmax(nb):stepj(nb)), &
                                     igextmin(nb),igextmax(nb),stepi(nb), &
                                     jgextmin(nb),jgextmax(nb),stepj(nb), &
                                     kextmin(nb),kextmax(nb),stepk(nb),   &
                                     kmax,0,l_out_nc4par= .TRUE.,comm_active=comm_zoom, &
                                     zaxis_var=sigecrit(1:kmax))
     ELSE
! If none l_out_nc4par argument, l_out_nc4par is considered as false
! Each cpu create its own output file inside the local domain (iextmin)
       CALL ionc4_createfile(fileres,lon2d(iextmin(nb):iextmax(nb):stepi(nb),jextmin(nb):jextmax(nb):stepj(nb)), &
                       lat2d(iextmin(nb):iextmax(nb):stepi(nb),jextmin(nb):jextmax(nb):stepj(nb)), &
                       iextmin(nb),iextmax(nb),stepi(nb),jextmin(nb),jextmax(nb),stepj(nb), &
                       kextmin(nb),kextmax(nb),stepk(nb),kmax,0,zaxis_var=sigecrit(1:kmax))
     ENDIF
!     CALL ionc4_write_lon(fileres,lon(iextmin(nb):iextmax(nb)),iextmin(nb),iextmax(nb))

!     CALL ionc4_torigine(fileres,date_ref)

! save attribute of the spatial extension of the domain
! -----------------------------------------------------
     CALL ionc4_gatt(fileres,'global_imin',igextmin(nb))
     CALL ionc4_gatt(fileres,'global_imax',igextmax(nb))
     CALL ionc4_gatt(fileres,'global_jmin',jgextmin(nb))
     CALL ionc4_gatt(fileres,'global_jmax',jgextmax(nb))

! save longitude, latitude fields and variables relative to the mesh grid
! -----------------------------------------------------------------------
     CALL out_grid(fileres)

! Add the vertical axis at w location
     IF (l_out_zaxis_w(nb)) CALL out_zaxis_w(fileres,kextmin(nb),kextmax(nb),stepk(nb),comm_zoom=comm_zoom)

     CALL ionc4_createvar(fileres,'SIG',' ','sigma variable',&

                          standard_name= 'ocean_sigma_variable',&

                          l_auxcoordinate=.false.,&
                          dims='z',valid_min= -1.0_riosh,valid_max= 0.0_riosh,fill_value= rg_valmanq_io)
     CALL ionc4_write_sig(fileres,'SIG',sigecrit,kmax)

     IF (l_out_zaxis_w(nb)) THEN
       CALL ionc4_createvar(fileres,'SIG_W',' ','sigma variable at the interface',&

                          standard_name= 'ocean_sigma_variable_at_w_location',&

                          l_auxcoordinate=.false.,location_in_cell='w',&
                          dims='z',valid_min= -1.0_riosh,valid_max= 0.0_riosh,fill_value= rg_valmanq_io)
       sigecrit(:)=sigw(:)
       CALL ionc4_write_sigw(fileres,'SIG_W',sigecrit,kmax)
     END IF



! Add attributes to follow OCO netCDF convention
     CALL out_netcdf_conv(fileres,nb,comm_zoom,igextmin(nb),igextmax(nb),jgextmin(nb),jgextmax(nb))

     IF (l_out_nc4par .AND. comm_zoom/=MPI_COMM_MARS) CALL MPI_COMM_FREE(comm_zoom,ierr_mpi)


! Save the bathymetry
! -------------------

! h0
! --
     IF (l_obc_tide) THEN
       chain='model_sea_floor_depth_below_sea_level'
     ELSE
       chain='model_sea_floor_depth_below_geoid'
     END IF
     CALL ionc4_createvar(fileres,name_out_h0,units='m',&
                           long_name='bathymetry relative to the mean level',&
                           standard_name = TRIM(chain),&
                           fill_value = -rg_valmanq_io,&
                           valid_min = riog_valid_min_h0,&
                           valid_max = riog_valid_max_h0,&
                           dims = 'xy',&
                           l_pack = l_out_pack,&
                           l_out_nc4par = l_out_nc4par) 
     tecrit(:,:)=h0(limin:limax,ljmin:ljmax)
     CALL ionc4_write_xy(fileres,TRIM(name_out_h0),tecrit,limin,limax,ljmin,ljmax,&
                          fill_value=-rg_valmanq_io,l_out_nc4par = l_out_nc4par )


! hx
! --
     IF (l_obc_tide) THEN
       chain='model_sea_floor_depth_below_sea_level_at_u_location'
     ELSE
       chain='model_sea_floor_depth_below_geoid_at_u_location'
     END IF
     CALL ionc4_createvar(fileres,name_out_hx,units='m',location_in_cell='u',&
                           long_name = 'bathymetry at u-location relative to the mean level',&
                           standard_name = TRIM(chain),&
                           fill_value = -REAL(rg_valmanq_io+h0fond,riosh),&
                           valid_min = riog_valid_min_h0,&
                           valid_max = riog_valid_max_h0,&
                           dims = 'xy',&
                           l_pack = l_out_pack,&
                           l_out_nc4par = l_out_nc4par)
     tecrit(:,:)=hx(limin:limax,ljmin:ljmax)
     CALL ionc4_write_xy(fileres,TRIM(name_out_hx),tecrit,limin,limax,ljmin,ljmax ,&
                         fill_value = -REAL(rg_valmanq_io+h0fond,riosh),l_out_nc4par = l_out_nc4par)


! hy
! --
     IF (l_obc_tide) THEN
       chain='model_sea_floor_depth_below_sea_level_at_v_location'
     ELSE
       chain='model_sea_floor_depth_below_geoid_at_v_location'
     END IF
     CALL ionc4_createvar(fileres,name_out_hy,units='m',location_in_cell='v',&
                           long_name = 'bathymetry at v-location relative to the mean level',&
                           standard_name = TRIM(chain),&
                           fill_value = -REAL(rg_valmanq_io+h0fond,riosh),&
                           valid_min = riog_valid_min_h0,&
                           valid_max = riog_valid_max_h0,&
                           dims = 'xy',&
                           l_pack = l_out_pack,&
                           l_out_nc4par = l_out_nc4par)
     tecrit(:,:)=hy(limin:limax,ljmin:ljmax)
     CALL ionc4_write_xy(fileres,TRIM(name_out_hy),tecrit,limin,limax,ljmin,ljmax,&
                         fill_value = -REAL(rg_valmanq_io+h0fond,riosh),l_out_nc4par = l_out_nc4par)


!Z0b
     name_out_z0b="Z0B"
     CALL ionc4_createvar(fileres,name_out_z0b,units='m',location_in_cell='v',&
                           long_name = 'bottom roughness',&
                           standard_name = "z0b",&
                           fill_value = -REAL(rg_valmanq_io+h0fond,riosh),&
                           valid_min =  0.0_riosh,&
                           valid_max = 1.0_riosh,&
                           dims = 'xy',&
                           l_pack = l_out_pack,&
                           l_out_nc4par = l_out_nc4par)
     tecrit(:,:)=z0b(limin:limax,ljmin:ljmax)
     CALL ionc4_write_xy(fileres,TRIM(name_out_z0b),tecrit,limin,limax,ljmin,ljmax,&
                         fill_value = -REAL(rg_valmanq_io+h0fond,riosh),l_out_nc4par = l_out_nc4par)





! Create time-depending variables
!--------------------------------

    IF (.NOT. l_out_timeave(nb) ) THEN
     IF (l_out_ssh(nb)) THEN
       IF (l_obc_tide) THEN
         chain='sea_surface_height_above_sea_level'
       ELSE
         chain='sea_surface_height_above_geoid'
       END IF
       CALL ionc4_createvar(fileres,name_out_xe,units= 'm',&
                            long_name = 'sea surface height',&
                            standard_name = TRIM(chain),&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_xe,&
                            valid_max = riog_valid_max_xe,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)    

!       CALL ionc4_vatt_axis(fileres,TRIM(name_out_xe),'Z')
!       CALL ionc4_vatt_positive(fileres,TRIM(name_out_xe),'up')
     END IF

     IF (l_out_ubt(nb)) THEN
       CALL ionc4_createvar(fileres,name_out_u,units= 'm s-1',&
                            long_name = 'barotropic zonal velocity',&
                            standard_name = 'barotropic_sea_water_x_velocity_at_u_location',&
                            location_in_cell='u',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_u,&
                            valid_max = riog_valid_max_u,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
       CALL ionc4_createvar(fileres,name_out_v,units= 'm s-1',&
                            long_name = 'barotropic meridional velocity',&
                            standard_name = 'barotropic_sea_water_y_velocity_at_v_location',&
                            location_in_cell='v',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_u,&
                            valid_max = riog_valid_max_u,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)                              
     END IF

     IF (l_out_u3d(nb).AND..NOT.l_modele2d) THEN
       CALL ionc4_createvar(fileres,name_out_uz,units= 'm s-1',&
                            long_name = '3d zonal velocity',&
                            standard_name = 'sea_water_x_velocity_at_u_location',&
                            location_in_cell='u',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_u,&
                            valid_max = riog_valid_max_u,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
       CALL ionc4_createvar(fileres,name_out_vz,units= 'm s-1',&
                            long_name = '3d meridional velocity',&
                            standard_name = 'sea_water_y_velocity_at_v_location',&
                            location_in_cell='v',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_u,&
                            valid_max = riog_valid_max_u,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)

     END IF

     IF (l_out_sal(nb).AND.l_modele2d) THEN
       CALL ionc4_createvar(fileres,name_out_sal,units= '1e-3',&
                            long_name = 'salinity',&
                            standard_name = 'sea_water_salinity',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_sal,&
                            valid_max = riog_valid_max_sal,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par) 
     ELSE IF (l_out_sal(nb).AND..NOT.l_modele2d) THEN
       CALL ionc4_createvar(fileres,name_out_sal,units= '1e-3',&
                            long_name = 'salinity',&
                            standard_name = 'sea_water_salinity',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_sal,&
                            valid_max = riog_valid_max_sal,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par) 
     END IF

     IF (l_out_temp(nb).AND.l_modele2d) THEN
       CALL ionc4_createvar(fileres,name_out_temp,units= 'degrees_Celsius',&
                            long_name = 'temperature',&
                            standard_name = 'sea_water_potential_temperature',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_temp,&
                            valid_max = riog_valid_max_temp,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par) 
     ELSE IF (l_out_temp(nb).AND..NOT.l_modele2d) THEN
       CALL ionc4_createvar(fileres,name_out_temp,units= 'degrees_Celsius',&
                            long_name = 'temperature',&
                            standard_name = 'sea_water_potential_temperature',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_temp,&
                            valid_max = riog_valid_max_temp,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par) 
     END IF

     IF (l_out_visc(nb).AND.l_modele2d) THEN
       CALL ionc4_createvar(fileres,name_out_vis2d,units= 'm2 s-1',&
                            long_name = 'viscosity',&
                            standard_name = 'ocean_tracer_xy_laplacian_diffusivity',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_visc,&
                            valid_max = riog_valid_max_visc,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par) 
     ELSE IF (l_out_visc(nb).AND..NOT.l_modele2d) THEN
       CALL ionc4_createvar(fileres,name_out_vishxe,units= 'm2 s-1',&
                            long_name = 'horizontal viscosity',&
                            standard_name = 'ocean_momentum_xy_laplacian_diffusivity',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_visc,&
                            valid_max = riog_valid_max_visc,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)      
       CALL ionc4_createvar(fileres,name_out_vishphi,units= 'm2 s-1',&
                            long_name = 'horizontal viscosity',&
                            standard_name = 'ocean_momentum_xy_laplacian_diffusivity_at_f_location',&
                            location_in_cell='f',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_visc,&
                            valid_max = riog_valid_max_visc,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par) 
     END IF

     IF (l_out_turb(nb).AND..NOT.l_modele2d) THEN
       CALL ionc4_createvar(fileres,name_out_nz,units= 'm2 s-1',&
                            long_name = 'vertical viscosity',&
                            standard_name = 'ocean_upward_momentum_diffusivity_at_w_location',&
                            location_in_cell='w',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_nz,&
                            valid_max = riog_valid_max_nz,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par) 
       CALL ionc4_createvar(fileres,name_out_kz,units= 'm2 s-1',&
                            long_name = 'vertical diffusivity',&
                            standard_name = 'ocean_upward_tracer_diffusivity_at_w_location',&
                            location_in_cell='w',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_kz,&
                            valid_max = riog_valid_max_kz,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)


     END IF

     IF (l_out_buoy(nb).AND..NOT.l_modele2d) THEN
       CALL ionc4_createvar(fileres,name_out_bz,units= 'kg m-3',&
                            long_name = 'density of water at atmospheric pressure minus 1000 kg m-3',&
                            standard_name = 'sea_water_sigma_t',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_bz,&
                            valid_max = riog_valid_max_bz,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
     END IF

     IF (l_out_wind(nb)) THEN
       CALL ionc4_createvar(fileres,name_out_tosx,units= 'N m-2',&
                            long_name = 'U-component of Surface Wind Stress',&
                            standard_name = 'surface_downward_x_stress_at_u_location',&
                            location_in_cell='u',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wstress,&
                            valid_max = riog_valid_max_wstress,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
       CALL ionc4_createvar(fileres,name_out_tosy,units= 'N m-2',&
                            long_name = 'V-component of Surface Wind Stress',&
                            standard_name = 'surface_downward_y_stress_at_v_location',&
                            location_in_cell='v',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wstress,&
                            valid_max = riog_valid_max_wstress,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
     END IF






     IF (l_out_wave(nb)) THEN
       IF (l_wave_rdir) &
       CALL ionc4_createvar(fileres,name_out_wave_dir,units= 'degrees',&
                            long_name = 'Direction of wind waves and swell (oceanographic convention : from)',&
!standard_name = 'sea_surface_wave_to_direction_due_to_wind_and_swell_waves',&
                            standard_name = 'sea_surface_wave_to_direction',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_dir,&
                            valid_max = riog_valid_max_wave_dir,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
       IF (l_wave_rhs) &
       CALL ionc4_createvar(fileres,name_out_wave_hs,units= 'm',&
                            long_name = 'Significant height of combined wind waves and swell',&
!standard_name = 'sea_surface_wave_significant_height_due_to_wind_and_swell_waves',&
                            standard_name = 'sea_surface_wave_significant_height',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_hs,&
                            valid_max = riog_valid_max_wave_hs,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
       IF (l_wave_rfoc) &
       CALL ionc4_createvar(fileres,name_out_wave_foc,units= 'W m-2',&
                            long_name = 'Wave mixing kinetic turbulent energy due to surface breaking wave',&
                            standard_name = 'turbulent_kinetic_energy_flux_into_sea_water_due_to_surface_dissipation_wave',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_foc,&
                            valid_max = riog_valid_max_wave_foc,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)       
       IF (l_wave_rtauox .AND. l_wave_rtauoy) THEN
         CALL ionc4_createvar(fileres,name_out_wave_tauox,units= 'm2 s-2',&
                            long_name = 'zonal component of the surface wave-ocean momentum flux',&
                            standard_name = 'x_component_of_surface_wave_ocean_momentum_flux',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_two,&
                            valid_max = riog_valid_max_wave_two,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)       
         CALL ionc4_createvar(fileres,name_out_wave_tauoy,units= 'm2 s-2',&
                            long_name = 'meridional component of the surface wave-ocean momentum flux',&
                            standard_name = 'y_component_of_surface_wave_ocean_momentum_flux',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_two,&
                            valid_max = riog_valid_max_wave_two,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)       
       END IF
       IF (l_wave_rubot) &
       CALL ionc4_createvar(fileres,name_out_wave_ubot,units= 'm s-1',&
                            long_name = 'Near bottom rms velocities',&
                            standard_name = 'rms_amplitude_of_orbital_bottom_velocity_due_to_wave',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_ubot,&
                            valid_max = riog_valid_max_wave_ubot,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
       IF (l_wave_ruubot .AND. l_wave_rvubot) THEN
          CALL ionc4_createvar(fileres,name_out_wave_ubotx,units= 'm s-1',&
                            long_name = 'Near bottom rms zonal velocities',&
                            standard_name = 'rms_amplitude_of_zonal_orbital_bottom_velocity_due_to_wave',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_ubot,&
                            valid_max = riog_valid_max_wave_ubot,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
          CALL ionc4_createvar(fileres,name_out_wave_uboty,units= 'm s-1',&
                            long_name = 'Near bottom rms meridional velocities',&
                            standard_name = 'rms_amplitude_of_meridional_orbital_bottom_velocity_due_to_wave',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_ubot,&
                            valid_max = riog_valid_max_wave_ubot,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
       END IF
       IF (l_wave_rt0m1) &
       CALL ionc4_createvar(fileres,name_out_wave_t0m1,units= 's',&
                            long_name = 'Mean wave period',&
                            standard_name = 'mean_wave_period',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_t0m,&
                            valid_max = riog_valid_max_wave_t0m,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)       
       IF (l_wave_rtauwx .AND. l_wave_rtauwy) THEN
         CALL ionc4_createvar(fileres,name_out_wave_tauwx,units= 'm2 s-2',&
                            long_name = 'zonal component of the wind stress supported by wave',&
                            standard_name = 'x_component_of_wave_supported_wind_stress',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_two,&
                            valid_max = riog_valid_max_wave_two,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)       
         CALL ionc4_createvar(fileres,name_out_wave_tauwy,units= 'm2 s-2',&
                            long_name = 'meridional component of the wind stress supported by wave',&
                            standard_name = 'y_component_of_wave_supported_wind_stress',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_two,&
                            valid_max = riog_valid_max_wave_two,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)       
       END IF
       IF (l_wave_rfbb) &
       CALL ionc4_createvar(fileres,name_out_wave_fbb,units= 'm3 s-3',&
                            long_name = 'Wave mixing kinetic turbulent energy due to wave bottom dissipation',&
                            standard_name = 'turbulent_kinetic_energy_flux_into_sea_water_due_to_wave_bottom_dissipation',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_foc,&
                            valid_max = riog_valid_max_wave_foc,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)       
       IF (l_wave_rtbb) THEN
         CALL ionc4_createvar(fileres,name_out_wave_tbbx,units= 'm2 s-2',&
                            long_name = 'zonal component of the bottom wave-ocean momentum flux',&
                            standard_name = 'x_component_of_bottom_wave_ocean_momentum_flux',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_two,&
                            valid_max = riog_valid_max_wave_two,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)       
         CALL ionc4_createvar(fileres,name_out_wave_tbby,units= 'm2 s-2',&
                            long_name = 'meridional component of the bottom wave-ocean momentum flux',&
                            standard_name = 'y_component_of_bottom_wave_ocean_momentum_flux',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_two,&
                            valid_max = riog_valid_max_wave_two,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)       
       END IF
       IF (l_wave_rtus) THEN
         CALL ionc4_createvar(fileres,name_out_wave_tusx,units= 'm2 s-1',&
                            long_name = 'eastward stokes transport',&
                            standard_name = 'eastward_stokes_transport',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_stokes,&
                            valid_max = riog_valid_max_wave_stokes,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)       
         CALL ionc4_createvar(fileres,name_out_wave_tusy,units= 'm2 s-1',&
                            long_name = 'northward stokes transport',&
                            standard_name = 'northward_stokes_transport',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_stokes,&
                            valid_max = riog_valid_max_wave_stokes,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)       
       END IF
       IF (l_wave_russ) THEN
         CALL ionc4_createvar(fileres,name_out_wave_ussx,units= 'm s-1',&
                            long_name = 'eastward surface stokes drift',&
                            standard_name = 'eastward_surface_stokes_drift',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_stokes,&
                            valid_max = riog_valid_max_wave_stokes,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)       
         CALL ionc4_createvar(fileres,name_out_wave_ussy,units= 'm s-1',&
                            long_name = 'northward surface stokes drift',&
                            standard_name = 'northward_surface_stokes_drift',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_stokes,&
                            valid_max = riog_valid_max_wave_stokes,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)       

       END IF
       IF (l_wave_rwlm) &
       CALL ionc4_createvar(fileres,name_out_wave_wlm,units= 'm',&
                            long_name = 'Mean wave length',&
                            standard_name = 'mean_wave_length',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_wlm,&
                            valid_max = riog_valid_max_wave_wlm,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)       

       IF (l_wave_rbhd) &
       CALL ionc4_createvar(fileres,name_out_wave_bhd,units= 'N m-1',&
                            long_name = 'Wave-induced Bernoulli head pressure',&
                            standard_name = 'bernoulli head pressure',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wave_bhd,&
                            valid_max = riog_valid_max_wave_bhd,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)       
       IF (l_wave_rradstr) THEN
         CALL ionc4_createvar(fileres,name_out_wave_sx,units= 'N m-1',&
                            long_name = 'zonal component of the radiation stress',&
                            standard_name = 'x_component_of_radiation_stress',&
                            fill_value = rg_valmanq_io,&
                            valid_min = -3000.0_riosh,&
                            valid_max = 3000.0_riosh,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
         CALL ionc4_createvar(fileres,name_out_wave_sy,units= 'N m-1',&
                            long_name = 'meridional component of the radiation stres',&
                            standard_name = 'y_component_of_radiation_stress',&
                            fill_value = rg_valmanq_io,&
                            valid_min = -3000.0_riosh,&
                            valid_max = 3000.0_riosh,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
       END IF

     END IF


    END IF   ! end test on .NOT.l_out_timeave(nb)

     IF (nbave /= 0 .AND. l_out_timeave(nb) ) THEN
       IF (l_out_ssh(nb)) THEN
         chain='sea_surface_height_above_geoid'
         CALL ionc4_createvar(fileres,name_out_xe,units= 'm',&
                            long_name = 'time averaged sea surface height',&
                            standard_name = TRIM(chain),&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_xe,&
                            valid_max = riog_valid_max_xe,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)

!         CALL ionc4_vatt_axis(fileres,TRIM(name_out_xe),'Z')
!         CALL ionc4_vatt_positive(fileres,TRIM(name_out_xe),'up')
       END IF
       IF (l_out_ubt(nb)) THEN
         CALL ionc4_createvar(fileres,name_out_u,units= 'm s-1',&
                            long_name = 'time averaged barotropic zonal velocity',&
                            standard_name = 'barotropic_sea_water_x_velocity_at_u_location',&
                            location_in_cell='u',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_u,&
                            valid_max = riog_valid_max_u,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
         CALL ionc4_createvar(fileres,name_out_v,units= 'm s-1',&
                            long_name = 'time averaged barotropic meridional velocity',&
                            standard_name = 'barotropic_sea_water_y_velocity_at_v_location',&
                            location_in_cell='v',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_u,&
                            valid_max = riog_valid_max_u,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)                         
       END IF
       IF (l_out_u3d(nb).AND..NOT.l_modele2d) THEN
         CALL ionc4_createvar(fileres,name_out_uz,units= 'm s-1',&
                            long_name = 'time averaged 3d zonal velocity',&
                            standard_name = 'sea_water_x_velocity_at_u_location',&
                            location_in_cell='u',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_u,&
                            valid_max = riog_valid_max_u,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
         CALL ionc4_createvar(fileres,name_out_vz,units= 'm s-1',&
                            long_name = 'time averaged 3d meridional velocity',&
                            standard_name = 'sea_water_y_velocity_at_v_location',&
                            location_in_cell='v',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_u,&
                            valid_max = riog_valid_max_u,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)

       END IF
       IF (l_out_sal(nb).AND.l_modele2d) THEN
         CALL ionc4_createvar(fileres,name_out_sal,units= '1e-3',&
                            long_name = 'time averaged salinity',&
                            standard_name = 'sea_water_salinity',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_sal,&
                            valid_max = riog_valid_max_sal,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
       ELSE IF (l_out_sal(nb).AND..NOT.l_modele2d) THEN
         CALL ionc4_createvar(fileres,name_out_sal,units= '1e-3',&
                            long_name = 'time averaged salinity',&
                            standard_name = 'sea_water_salinity',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_sal,&
                            valid_max = riog_valid_max_sal,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
       END IF
       IF (l_out_temp(nb).AND.l_modele2d) THEN
         CALL ionc4_createvar(fileres,name_out_temp,units= 'degrees_Celsius',&
                            long_name = 'time averaged temperature',&
                            standard_name = 'sea_water_potential_temperature',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_temp,&
                            valid_max = riog_valid_max_temp,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
       ELSE IF (l_out_temp(nb).AND..NOT.l_modele2d) THEN
         CALL ionc4_createvar(fileres,name_out_temp,units= 'degrees_Celsius',&
                            long_name = 'time averaged temperature',&
                            standard_name = 'sea_water_potential_temperature',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_temp,&
                            valid_max = riog_valid_max_temp,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
       END IF
       IF (l_out_visc(nb).AND.l_modele2d) THEN
         CALL ionc4_createvar(fileres,name_out_vis2d,units= 'm2 s-1',&
                            long_name = 'time averaged viscosity',&
                            standard_name = 'ocean_tracer_xy_laplacian_diffusivity',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_visc,&
                            valid_max = riog_valid_max_visc,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
       ELSE IF (l_out_visc(nb).AND..NOT.l_modele2d) THEN
         CALL ionc4_createvar(fileres,name_out_vishxe,units= 'm2 s-1',&
                            long_name = 'time averaged horizontal viscosity',&
                            standard_name = 'ocean_momentum_xy_laplacian_diffusivity',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_visc,&
                            valid_max = riog_valid_max_visc,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
         CALL ionc4_createvar(fileres,name_out_vishphi,units= 'm2 s-1',&
                            long_name = 'time average horizontal viscosity',&
                            standard_name = 'ocean_momentum_xy_laplacian_diffusivity_at_f_location',&
                            location_in_cell='f',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_visc,&
                            valid_max = riog_valid_max_visc,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
       END IF
       IF (l_out_turb(nb).AND..NOT.l_modele2d) THEN
         CALL ionc4_createvar(fileres,name_out_nz,units= 'm2 s-1',&
                            long_name = 'time averaged vertical viscosity',&
                            standard_name = 'ocean_upward_momentum_diffusivity_at_w_location',&
                            location_in_cell='w',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_nz,&
                            valid_max = riog_valid_max_nz,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
         CALL ionc4_createvar(fileres,name_out_kz,units= 'm2 s-1',&
                            long_name = 'time averaged vertical diffusivity',&
                            standard_name = 'ocean_upward_tracer_diffusivity_at_w_location',&
                            location_in_cell='w',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_kz,&
                            valid_max = riog_valid_max_kz,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)


       END IF
       IF (l_out_buoy(nb).AND..NOT.l_modele2d) THEN
         CALL ionc4_createvar(fileres,name_out_bz,units= 'kg m-3',&
                            long_name = 'time averaged density of water at atmospheric pressure minus 1000 kg m-3',&
                            standard_name = 'sea_water_sigma_t',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_bz,&
                            valid_max = riog_valid_max_bz,&
                            dims= 'zxyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)
       END IF


     END IF


     IF (nbatm /= 0 .AND. l_out_atmsflx(nb) ) THEN
       CALL ionc4_createvar(fileres,name_out_mrad,units= 'W m-2',&
                            long_name = 'radiative flux',&
                            standard_name = 'average_of_surface_net_downward_shortwave_flux',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_flux,&
                            valid_max = riog_valid_max_flux,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)      
       CALL ionc4_createvar(fileres,name_out_mirlos,units= 'W m-2',&
                            long_name = 'thermal flux',&
                            standard_name = 'average_of_surface_net_downward_longwave_flux',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_flux,&
                            valid_max = riog_valid_max_flux,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)                               
       CALL ionc4_createvar(fileres,name_out_mlatlos,units= 'W m-2',&
                            long_name = 'latent flux',&
                            standard_name = 'average_of_surface_downward_latent_heat_flux',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_flux,&
                            valid_max = riog_valid_max_flux,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)                               
       CALL ionc4_createvar(fileres,name_out_msenslos,units= 'W m-2',&
                            long_name = 'sensible flux',&
                            standard_name = 'average_of_surface_downward_sensible_heat_flux',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_flux,&
                            valid_max = riog_valid_max_flux,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)                               
       CALL ionc4_createvar(fileres,name_out_mevap,units= 'm',&
                            long_name = 'Evaporation',&
                            standard_name = 'lwe_thickness_of_upward_water_evaporation_amount',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_flux,&
                            valid_max = riog_valid_max_flux,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)                               
       CALL ionc4_createvar(fileres,name_out_mrain,units= 'm',&
                            long_name = 'rain',&
                            standard_name = 'lwe_thickness_of_downward_precipitation_amount',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_flux,&
                            valid_max = riog_valid_max_flux,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)                               
       CALL ionc4_createvar(fileres,name_out_mtosx,units= 'N m-2',&
                            long_name = 'Averaged U-component of Surface Wind Stress',&
                            standard_name = 'average_of_surface_downward_x_stress_at_u_location',&
                            location_in_cell='u',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wstress,&
                            valid_max = riog_valid_max_wstress,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)                               
       CALL ionc4_createvar(fileres,name_out_mtosy,units= 'N m-2',&
                            long_name = 'Averaged V-component of Surface Wind Stress',&
                            standard_name = 'average_of_surface_downward_y_stress_at_v_location',&
                            location_in_cell='v',&
                            fill_value = rg_valmanq_io,&
                            valid_min = riog_valid_min_wstress,&
                            valid_max = riog_valid_max_wstress,&
                            dims= 'xyt',&
                            l_pack = l_out_pack,&
                            l_out_nc4par = l_out_nc4par)                               


     END IF





     DO l = 1,nv_state
      IF (l_out_subs(l)) THEN

       IF (l_out_log) THEN
         unit_var(irk_fil(l))= 'log('//TRIM(unit_var(irk_fil(l)))//')' 
       ELSE
         unit_var(irk_fil(l))=TRIM(unit_var(irk_fil(l)))
       ENDIF
       IF (l_modele2d) THEN
         CALL ionc4_createvar(fileres,name_var(irk_fil(l)),units= TRIM(unit_var(irk_fil(l))),&
                              long_name = TRIM(long_name_var(irk_fil(l))),&
                              standard_name = TRIM(standard_name_var(irk_fil(l))),&
                              fill_value = rg_valmanq_io,&
                              valid_min = REAL(valid_min_var(irk_fil(l)),riosh),&
                              valid_max = REAL(valid_max_var(irk_fil(l)),riosh),&
                              dims= 'xyt',&
                              l_pack = .false.,&
                              l_out_nc4par = l_out_nc4par) 
       ELSE
         CALL ionc4_createvar(fileres,name_var(irk_fil(l)),units= TRIM(unit_var(irk_fil(l))),&
                              long_name = TRIM(long_name_var(irk_fil(l))),&
                              standard_name = TRIM(standard_name_var(irk_fil(l))),&
                              fill_value = rg_valmanq_io,&
                              valid_min = REAL(valid_min_var(irk_fil(l)),riosh),&
                              valid_max = REAL(valid_max_var(irk_fil(l)),riosh),&
                              dims= 'zxyt',&
                              l_pack = .false.,&
                              l_out_nc4par = l_out_nc4par)
!MD 08/06/2018
!         write(100,*) 'rentre dans IF'
         IF (l_out_ws(nb)) THEN
             IF (l >= imud1 .AND. l <= imud2) THEN
!             write(101,*) name_var(irk_fil(l))
             name_var_sf = TRIM('WS_'//name_var(irk_fil(l)))
             long_name_var_sf = TRIM('WS_'//long_name_var(irk_fil(l)))
             standard_name_var_sf = TRIM('WS_'//standard_name_var(irk_fil(l)))
             CALL ionc4_createvar(fileresu(nb),name_var_sf,units= 'm s-1',&
                                  long_name = long_name_var_sf,&
                                  standard_name = standard_name_var_sf,&
                                  fill_value = rg_valmanq_io,&
                                  valid_min = REAL(-1.0e9,riosh),&
                                  valid_max = REAL(1.0e9,riosh),&
                                  dims= 'zxyt',&
                                  l_pack = .false.,&
                                  l_out_nc4par = l_out_nc4par)
             END IF
         END IF
         IF (l_out_varint(nb)) THEN
             name_var_sf1 = TRIM('C_int_'//name_var(irk_fil(l)))
             long_name_var_sf1 = TRIM('C_int_'//long_name_var(irk_fil(l)))
             standard_name_var_sf1 = TRIM('C_int_'//standard_name_var(irk_fil(l)))
             CALL ionc4_createvar(fileresu(nb),name_var_sf1,units= 'kg m-3',&
                                  long_name = long_name_var_sf1,&
                                  standard_name = standard_name_var_sf1,&
                                  fill_value = rg_valmanq_io,&
                                  valid_min = REAL(-1.0e9,riosh),&
                                  valid_max = REAL(1.0e9,riosh),&
                                  dims= 'xyt',&
                                  l_pack = .false.,&
                                  l_out_nc4par = l_out_nc4par)
         ENDIF
!END MD 08/06/2018
       ENDIF
      ENDIF
     ENDDO
     

     
     


!------------------------------------
! sauvegarde de la condition initiale
!------------------------------------

!     IF ( (.not.l_initfromfile .and. date_startshift == 0.0_rlg) .OR. (l_filebydate(nb) .AND. nbouc /= 0) .OR. (l_initfromfile .AND. nbouc== 0) ) THEN
     IF (tect >= tsort(nb).AND.tect <= tfinsor(nb)) THEN
       CALL ionc4_write_time(fileres,0,tect)
       
      IF (.NOT. l_out_timeave(nb) ) THEN   ! instantaneous outputs

       IF (l_out_ssh(nb)) CALL out_write_xytvar(fileres,name_out_xe,rg_valmanq_io,mask_h0, &
                                                x(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)

       IF (l_out_ubt(nb)) THEN
         CALL out_write_xytvar(fileres,name_out_u,rg_valmanq_io,mask_hx, &
                               u(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         CALL out_write_xytvar(fileres,name_out_v,rg_valmanq_io,mask_hy, &
                               v(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       END IF


       IF (l_out_u3d(nb) .AND. .NOT.l_modele2d) THEN
         CALL out_write_zxytvar(fileres,name_out_uz,rg_valmanq_io,mask_hx, &
                                uz(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
         CALL out_write_zxytvar(fileres,name_out_vz,rg_valmanq_io,mask_hy, &
                                vz(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)

       END IF

       IF (l_out_sal(nb) .AND. l_modele2d) THEN
         CALL out_write_xytvar(fileres,name_out_sal,rg_valmanq_io,mask_h0, &
                               sal(1,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       ELSE IF (l_out_sal(nb).AND..NOT.l_modele2d) THEN
         CALL out_write_zxytvar(fileres,name_out_sal,rg_valmanq_io,mask_h0, &
                                sal(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
       END IF

       IF (l_out_temp(nb) .AND. l_modele2d) THEN
         CALL out_write_xytvar(fileres,name_out_temp,rg_valmanq_io,mask_h0, &
                               temp(1,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       ELSE IF (l_out_temp(nb) .AND. .NOT.l_modele2d) THEN
         CALL out_write_zxytvar(fileres,name_out_temp,rg_valmanq_io,mask_h0, &
                                temp(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
       END IF

       IF (l_out_visc(nb) .AND. l_modele2d) THEN
         CALL out_write_xytvar(fileres,name_out_vis2d,rg_valmanq_io,mask_h0, &
                               vis2d(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       ELSE IF (l_out_visc(nb) .AND. .NOT.l_modele2d) THEN
         CALL out_write_zxytvar(fileres,name_out_vishxe,rg_valmanq_io,mask_h0, &
                                vish_xe(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
         CALL out_write_zxytvar(fileres,name_out_vishphi,rg_valmanq_io,mask_h0, &
                                vish_phi(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
       END IF

       IF (l_out_turb(nb) .AND. .NOT.l_modele2d) THEN
         CALL out_write_zxytvar(fileres,name_out_nz,rg_valmanq_io,mask_h0, &
                                nz(:,limin:limax,ljmin:ljmax),0,0,l_out_nc4par = l_out_nc4par)
         CALL out_write_zxytvar(fileres,name_out_kz,rg_valmanq_io,mask_h0, &
                                kz(:,limin:limax,ljmin:ljmax),0,0,l_out_nc4par = l_out_nc4par)


       END IF

       IF (l_out_buoy(nb) .AND. .NOT.l_modele2d) THEN
         tmp3d(:,:,:) = 1000.0_rsh - 1000.0_rsh*boz(:,limin:limax,ljmin:ljmax)/grav
         CALL out_write_zxytvar(fileres,name_out_bz,rg_valmanq_io,mask_h0, &
                                tmp3d,0,1,l_out_nc4par = l_out_nc4par)
       END IF

       IF (l_out_wind(nb)) THEN
         CALL out_write_xytvar(fileres,name_out_tosx,rg_valmanq_io,mask_h0, &
                               tosx(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         CALL out_write_xytvar(fileres,name_out_tosy,rg_valmanq_io,mask_h0, &
                               tosy(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       END IF







       IF (l_out_wave(nb)) THEN
         IF (l_wave_rdir) THEN
! trigonometric to oceanographic conventions : from north = 0, from west = 270
           tmp2d(limin:limax,ljmin:ljmax) = MOD(-wave_dir(limin:limax,ljmin:ljmax)*180.0_rsh/pi+270.0_rsh,360.0_rsh)
           CALL out_write_xytvar(fileres,name_out_wave_dir,rg_valmanq_io,mask_h0, &
                                 wave_dir(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_wave_rhs) &
         CALL out_write_xytvar(fileres,name_out_wave_hs,rg_valmanq_io,mask_h0, &
                               wave_hs(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         IF (l_wave_rfoc) &
         CALL out_write_xytvar(fileres,name_out_wave_foc,rg_valmanq_io,mask_h0, &
                               wave_phioc(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         IF (l_wave_rtauox .AND. l_wave_rtauoy) THEN
           CALL out_write_xytvar(fileres,name_out_wave_tauox,rg_valmanq_io,mask_h0, &
                               wave_tauox(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
           CALL out_write_xytvar(fileres,name_out_wave_tauoy,rg_valmanq_io,mask_h0, &
                               wave_tauoy(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_wave_rubot) &
         CALL out_write_xytvar(fileres,name_out_wave_ubot,rg_valmanq_io,mask_h0, &
                               wave_ubot(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         IF (l_wave_ruubot .AND. l_wave_rvubot) THEN
           CALL out_write_xytvar(fileres,name_out_wave_ubotx,rg_valmanq_io,mask_h0, &
                               wave_uubot(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
           CALL out_write_xytvar(fileres,name_out_wave_uboty,rg_valmanq_io,mask_h0, &
                               wave_vubot(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_wave_rt0m1) &
         CALL out_write_xytvar(fileres,name_out_wave_t0m1,rg_valmanq_io,mask_h0, &
                               wave_t0m1(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         IF (l_wave_rtauwx .AND. l_wave_rtauwy) THEN
           CALL out_write_xytvar(fileres,name_out_wave_tauwx,rg_valmanq_io,mask_h0, &
                               wave_tauwx(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
           CALL out_write_xytvar(fileres,name_out_wave_tauwy,rg_valmanq_io,mask_h0, &
                               wave_tauwy(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_wave_rfbb) &
         CALL out_write_xytvar(fileres,name_out_wave_fbb,rg_valmanq_io,mask_h0, &
                               wave_phibbl(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         IF (l_wave_rtbb) THEN
           CALL out_write_xytvar(fileres,name_out_wave_tbbx,rg_valmanq_io,mask_h0, &
                               wave_taubbl_x(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
           CALL out_write_xytvar(fileres,name_out_wave_tbby,rg_valmanq_io,mask_h0, &
                               wave_taubbl_y(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_wave_rtus) THEN
           CALL out_write_xytvar(fileres,name_out_wave_tusx,rg_valmanq_io,mask_h0, &
                               wave_tusx(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
           CALL out_write_xytvar(fileres,name_out_wave_tusy,rg_valmanq_io,mask_h0, &
                               wave_tusy(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_wave_russ) THEN
           CALL out_write_xytvar(fileres,name_out_wave_ussx,rg_valmanq_io,mask_h0, &
                               wave_ussx(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
           CALL out_write_xytvar(fileres,name_out_wave_ussy,rg_valmanq_io,mask_h0, &
                               wave_ussy(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)

         END IF
         IF (l_wave_rbhd) &
         CALL out_write_xytvar(fileres,name_out_wave_bhd,rg_valmanq_io,mask_h0, &
                               wave_bhd(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         IF (l_wave_rwlm) &
         CALL out_write_xytvar(fileres,name_out_wave_wlm,rg_valmanq_io,mask_h0, &
                               wave_wlm(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)

       END IF


    END IF   ! end test on .NOT.l_out_timeave(nb)

       IF (nbave /= 0 .AND. l_out_timeave(nb) ) THEN
         IF (l_out_ssh(nb)) CALL out_write_xytvar(fileres,name_out_xe,rg_valmanq_io,mask_h0, &
                                                ssh_ave(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)

         IF (l_out_ubt(nb)) THEN
           CALL out_write_xytvar(fileres,name_out_u,rg_valmanq_io,mask_hx, &
                               u_ave(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
           CALL out_write_xytvar(fileres,name_out_v,rg_valmanq_io,mask_hy, &
                               v_ave(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_out_u3d(nb) .AND. .NOT.l_modele2d) THEN
           CALL out_write_zxytvar(fileres,name_out_uz,rg_valmanq_io,mask_hx, &
                                uz_ave(nb,:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
           CALL out_write_zxytvar(fileres,name_out_vz,rg_valmanq_io,mask_hy, &
                                vz_ave(nb,:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)

         END IF

         IF (l_out_sal(nb) .AND. l_modele2d) THEN
           CALL out_write_xytvar(fileres,name_out_sal,rg_valmanq_io,mask_h0, &
                               sal_ave(nb,1,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         ELSE IF (l_out_sal(nb).AND..NOT.l_modele2d) THEN
           CALL out_write_zxytvar(fileres,name_out_sal,rg_valmanq_io,mask_h0, &
                                sal_ave(nb,:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_out_temp(nb) .AND. l_modele2d) THEN
           CALL out_write_xytvar(fileres,name_out_temp,rg_valmanq_io,mask_h0, &
                               temp_ave(nb,1,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         ELSE IF (l_out_temp(nb) .AND. .NOT.l_modele2d) THEN
           CALL out_write_zxytvar(fileres,name_out_temp,rg_valmanq_io,mask_h0, &
                                temp_ave(nb,:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
         END IF

         IF (l_out_visc(nb) .AND. l_modele2d) THEN
           CALL out_write_xytvar(fileres,name_out_vis2d,rg_valmanq_io,mask_h0, &
                               vis2d_ave(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         ELSE IF (l_out_visc(nb) .AND. .NOT.l_modele2d) THEN
           CALL out_write_zxytvar(fileres,name_out_vishxe,rg_valmanq_io,mask_h0, &
                                vish_xe_ave(nb,:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
           CALL out_write_zxytvar(fileres,name_out_vishphi,rg_valmanq_io,mask_h0, &
                                vish_phi_ave(nb,:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_out_turb(nb) .AND. .NOT.l_modele2d) THEN
           CALL out_write_zxytvar(fileres,name_out_nz,rg_valmanq_io,mask_h0, &
                                nz_ave(nb,:,limin:limax,ljmin:ljmax),0,0,l_out_nc4par = l_out_nc4par)
           CALL out_write_zxytvar(fileres,name_out_kz,rg_valmanq_io,mask_h0, &
                                kz_ave(nb,:,limin:limax,ljmin:ljmax),0,0,l_out_nc4par = l_out_nc4par)


         END IF

         IF (l_out_buoy(nb) .AND. .NOT.l_modele2d) THEN
           tmp3d(:,:,:) = 1000.0_rsh - 1000.0_rsh*boz_ave(nb,:,limin:limax,ljmin:ljmax)/grav
           CALL out_write_zxytvar(fileres,name_out_bz,rg_valmanq_io,mask_h0, &
                                tmp3d,0,1,l_out_nc4par = l_out_nc4par)
         END IF

! initialize averaged dynamical and thermodynamical variables
! -----------------------------------------------------------

         dtn = t - tsort(nb)
         IF (l_out_ssh(nb)) &
         ssh_ave(nb,:,:) = x(limin:limax,ljmin:ljmax)*dtn
         IF (l_out_ubt(nb)) THEN
           u_ave(nb,:,:) = u(limin:limax,ljmin:ljmax)*dtn
           v_ave(nb,:,:) = v(limin:limax,ljmin:ljmax)*dtn
         END IF
         IF (l_out_u3d(nb).AND..NOT.l_modele2d) THEN
           uz_ave(nb,:,:,:) = dtn * uz(:,limin:limax,ljmin:ljmax)
           vz_ave(nb,:,:,:) = dtn * vz(:,limin:limax,ljmin:ljmax)
         END IF
         IF (l_out_temp(nb)) &
           temp_ave(nb,:,:,:) = dtn * temp(:,limin:limax,ljmin:ljmax)
         IF (l_out_sal(nb)) &
           sal_ave(nb,:,:,:) = dtn * sal(:,limin:limax,ljmin:ljmax)
         IF (l_out_visc(nb).AND.l_modele2d) THEN
           vis2d_ave(nb,:,:) = dtn * vis2d(limin:limax,ljmin:ljmax)
         ELSE IF (l_out_visc(nb).AND..NOT.l_modele2d) THEN
           vish_xe_ave(nb,:,:,:) = dtn * vish_xe(:,limin:limax,ljmin:ljmax)
           vish_phi_ave(nb,:,:,:) = dtn * vish_phi(:,limin:limax,ljmin:ljmax)
         END IF
         IF (l_out_turb(nb).AND..NOT.l_modele2d) THEN
           nz_ave(nb,:,:,:) = dtn * nz(:,limin:limax,ljmin:ljmax)
           kz_ave(nb,:,:,:) = dtn * kz(:,limin:limax,ljmin:ljmax)
         END IF
         IF (l_out_buoy(nb).AND..NOT.l_modele2d) &
           boz_ave(nb,:,:,:) = dtn * boz(:,limin:limax,ljmin:ljmax)
         dt_ave(nb) =  dtn
       END IF

       IF (nbatm /= 0 .AND. l_out_atmsflx(nb) ) THEN
         CALL out_write_xytvar(fileres,name_out_mrad,rg_valmanq_io,mask_h0,     &
                               rad_cum(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         CALL out_write_xytvar(fileres,name_out_mirlos,rg_valmanq_io,mask_h0,   &
                               irlos_cum(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         CALL out_write_xytvar(fileres,name_out_mlatlos,rg_valmanq_io,mask_h0,  &
                               latlos_cum(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         CALL out_write_xytvar(fileres,name_out_msenslos,rg_valmanq_io,mask_h0, &
                               senslos_cum(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         CALL out_write_xytvar(fileres,name_out_mevap,rg_valmanq_io,mask_h0, &
                               evap_cum(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         CALL out_write_xytvar(fileres,name_out_mrain,rg_valmanq_io,mask_h0, &
                               rain_cum(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         CALL out_write_xytvar(fileres,name_out_mtosx,rg_valmanq_io,mask_h0, &
                               tosx_cum(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         CALL out_write_xytvar(fileres,name_out_mtosy,rg_valmanq_io,mask_h0, &
                               tosy_cum(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)


! initialize cumulative atmospherical surface fluxes
! --------------------------------------------------

         dtn = t - tsort(nb)
         rad_cum(nb,:,:)     = rad(limin:limax,ljmin:ljmax)    *dtn
         irlos_cum(nb,:,:)   = irlos(limin:limax,ljmin:ljmax)  *dtn
         latlos_cum(nb,:,:)  = latlos(limin:limax,ljmin:ljmax) *dtn
         senslos_cum(nb,:,:) = senslos(limin:limax,ljmin:ljmax)*dtn
         evap_cum(nb,:,:)    = evap(limin:limax,ljmin:ljmax)   *dtn
         rain_cum(nb,:,:)    = rain(limin:limax,ljmin:ljmax)   *dtn
         tosx_cum(nb,:,:)    = tosx(limin:limax,ljmin:ljmax)   *dtn
         tosy_cum(nb,:,:)    = tosy(limin:limax,ljmin:ljmax)   *dtn

         dt_cum(nb) =  dtn

       END IF




       DO l = 1,nv_state
        IF (l_out_subs(l)) THEN
         IF (l_modele2d) THEN
           tecrit(:,:)=rg_valmanq_io
           DO j=ljmin,ljmax
             DO i=MAX(limin,ig(j)),MIN(limax,id(j))
               IF ((j.GE.jb(i)).AND.(j.LE.jh(i))) THEN
                 IF (l_out_log) THEN
                   tecrit(i,j) = LOG(MAX(cv_wat(l,1,i,j)/unit_modif_mudbio_N2dw(irk_fil(l)),1.0e-30))/LOG(10.0_rsh)
                 ELSE
!tecrit(i,j) =MAX(cv_wat(l,1,i,j)/unit_modif_mudbio_N2dw(irk_fil(l)),1.0e-30)
                   tecrit(i,j) =cv_wat(l,1,i,j)/unit_modif_mudbio_N2dw(irk_fil(l))
                 END IF
                 IF (h0(i,j)+x(i,j) <= h0fond) tecrit(i,j)=rg_valmanq_io
               END IF
             END DO
           END DO
           CALL ionc4_write_xyt(fileres,TRIM(name_var(irk_fil(l))),tecrit,limin,limax,ljmin,ljmax,0,l_out_nc4par = l_out_nc4par)
         ELSE
           tzecrit(:,:,:)=rg_valmanq_io

! modif septembre 2014 : separation du test pour eviter depassement de tableaux sur l
           IF (l>=isand1 .AND. l<=isand2 ) THEN
             l_outsandrouse_l=l_outsandrouse(l)
           ELSE
             l_outsandrouse_l=.FALSE.
           ENDIF
           IF (l>=isand1 .AND. l<=isand2 .AND. l_outsandrouse_l .AND. l_subs2D(l)) THEN
             DO j=ljmin,ljmax
               DO i=MAX(limin,ig(j)),MIN(limax,id(j))
                 IF ((j.GE.jb(i)).AND.(j.LE.jh(i))) THEN
                   IF (l_out_log) THEN
                     tzecrit(:,i,j) = LOG(MAX(cv_wat(l,:,i,j),1.0e-30))/LOG(10.0_rsh)
                   ELSE
                     htott=-f_zw(h0(i,j),x(i,j),0,i,j)+x(i,j)
                     IF (htott > h0fond .and. sum_tmp(l,i,j) > 0.0_rsh) THEN
                       DO k=1,kmax
                         zk=-f_zu(h0(i,j),x(i,j),k,i,j)+x(i,j)
                         IF (zk > 0.0_rsh) THEN
                           zk=htott-zk
                           tzecrit(k,i,j)=cv_wat(l,1,i,j)*htott*((htott-zk)/zk)**rouse2D(l,i,j)
                           tzecrit(k,i,j)=tzecrit(k,i,j)/sum_tmp(l,i,j)
                         ELSE
                           tzecrit(k,i,j)=0.0_rsh
                         ENDIF
                       ENDDO 
                     ELSE
                       tzecrit(:,i,j)=0.0_rsh
                     ENDIF
                   END IF
!MD 08/06/2018
!                   DO k=1,kmax
!                      ws_mud(k,i,j)=ws3(k,5,i,j)
!                   END DO
!END MD 08/06/2018
                   IF (h0(i,j)+x(i,j) <= h0fond) tzecrit(:,i,j)=rg_valmanq_io
                 END IF
               END DO
             END DO
           ELSE

             DO j=ljmin,ljmax
             DO i=MAX(limin,ig(j)),MIN(limax,id(j))
               IF ((j.GE.jb(i)).AND.(j.LE.jh(i))) THEN
                 IF (l_out_log) THEN
                   tzecrit(:,i,j) = LOG(MAX(cv_wat(l,:,i,j)/unit_modif_mudbio_N2dw(irk_fil(l)),1.0e-30))/LOG(10.0_rsh)
                 ELSE
                   tzecrit(:,i,j) = cv_wat(l,:,i,j)/unit_modif_mudbio_N2dw(irk_fil(l))
                 END IF
                 IF (h0(i,j)+x(i,j) <= h0fond) tzecrit(:,i,j)=rg_valmanq_io
               END IF
             END DO
             END DO

           ENDIF

           CALL ionc4_write_zxyt(fileres,TRIM(name_var(irk_fil(l))),tzecrit,limin,limax,ljmin,ljmax,1,kmax,0, &
                                         l_out_nc4par = l_out_nc4par)
!MD 08/06/2018
!           write(102,*) 'RENTRE DANS IF 2'
           IF (l_out_ws(nb)) THEN
             IF (l >= imud1 .AND. l <= imud2) THEN
               name_var_sf = TRIM('WS_'//name_var(irk_fil(l)))
               DO j=ljmin,ljmax
               DO i=MAX(limin,ig(j)),MIN(limax,id(j))
                 IF ((j.GE.jb(i)).AND.(j.LE.jh(i))) THEN
                   IF (l_out_log) THEN
                     tzecrit(:,i,j) = LOG(ws3(:,l,i,j))/LOG(10.0_rsh)
                   ELSE
                     tzecrit(:,i,j) = ws3(:,l,i,j)
                   END IF
                   IF (h0(i,j)+x(i,j) <= h0fond) tzecrit(:,i,j)=rg_valmanq_io
                 END IF
               END DO
               END DO
               CALL ionc4_write_zxyt(fileres,name_var_sf,tzecrit,limin,limax,ljmin,ljmax,1,kmax,0, &
                                           l_out_nc4par = l_out_nc4par)
             END IF
           END IF
           IF (l_out_varint(nb)) THEN
               tecrit(:,:)=rg_valmanq_io 
               DO j=ljmin,ljmax
                   DO i=MAX(limin,ig(j)),MIN(limax,id(j))
                       IF ((j.GE.jb(i)).AND.(j.LE.jh(i))) THEN
!                   DO i=limin,limax
                           tecrit(i,j)=SUM(cv_wat(l,:,i,j))
                           IF (h0(i,j)+x(i,j) <= h0fond) tecrit(i,j)=rg_valmanq_io
                       END IF
                   END DO
               END DO
!               CALL ionc4_write_xyt(fileresu(nb),TRIM(name_var(irk_fil(l))),tecrit,limin,limax,ljmin,ljmax,0,l_out_nc4par = l_out_nc4par)
               name_var_sf1 = TRIM('C_int_'//name_var(irk_fil(l)))
               CALL ionc4_write_xyt(fileres,name_var_sf1,tecrit,limin,limax,ljmin,ljmax,0,l_out_nc4par = l_out_nc4par)
           END IF
!END MD 08/06/2018
         END IF
        END IF
       END DO
       

       
       


       IF (l_filebydate(nb)) THEN
         CALL ionc4_close(fileres)
       ELSE
! To write the data on the disk and not loose data in case of run crash
         CALL ionc4_sync(fileres)
       ENDIF

     END IF


   ELSE IF ( .NOT. l_filebydate(nb) ) THEN


     CALL ionc4_write_time(fileresu(nb),0,tect)

    IF (.NOT. l_out_timeave(nb) ) THEN  ! instantaneous outputs

     IF (l_out_ssh(nb)) CALL out_write_xytvar(fileresu(nb),name_out_xe,rg_valmanq_io,mask_h0, &
                                              x(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)

     IF (l_out_ubt(nb)) THEN
       CALL out_write_xytvar(fileresu(nb),name_out_u,rg_valmanq_io,mask_hx, &
                             u(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       CALL out_write_xytvar(fileresu(nb),name_out_v,rg_valmanq_io,mask_hy, &
                             v(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
     END IF

     IF (l_out_u3d(nb) .AND. .NOT.l_modele2d) THEN
       CALL out_write_zxytvar(fileresu(nb),name_out_uz,rg_valmanq_io,mask_hx, &
                              uz(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
       CALL out_write_zxytvar(fileresu(nb),name_out_vz,rg_valmanq_io,mask_hy, &
                              vz(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)

     END IF

     IF (l_out_sal(nb) .AND. l_modele2d) THEN
       CALL out_write_xytvar(fileresu(nb),name_out_sal,rg_valmanq_io,mask_h0, &
                             sal(1,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
     ELSE IF (l_out_sal(nb) .AND. .NOT.l_modele2d) THEN
       CALL out_write_zxytvar(fileresu(nb),name_out_sal,rg_valmanq_io,mask_h0, &
                              sal(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     END IF

     IF (l_out_temp(nb) .AND. l_modele2d) THEN
       CALL out_write_xytvar(fileresu(nb),name_out_temp,rg_valmanq_io,mask_h0, &
                             temp(1,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
     ELSE IF (l_out_temp(nb) .AND. .NOT.l_modele2d) THEN
       CALL out_write_zxytvar(fileresu(nb),name_out_temp,rg_valmanq_io,mask_h0, &
                              temp(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     END IF

     IF (l_out_visc(nb) .AND. l_modele2d) THEN
       CALL out_write_xytvar(fileresu(nb),name_out_vis2d,rg_valmanq_io,mask_h0, &
                             vis2d(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
     ELSE IF (l_out_visc(nb) .AND. .NOT.l_modele2d) THEN
       CALL out_write_zxytvar(fileresu(nb),name_out_vishxe,rg_valmanq_io,mask_h0, &
                              vish_xe(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
       CALL out_write_zxytvar(fileresu(nb),name_out_vishphi,rg_valmanq_io,mask_h0, &
                              vish_phi(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     END IF

     IF (l_out_turb(nb) .AND. .NOT.l_modele2d) THEN
       CALL out_write_zxytvar(fileresu(nb),name_out_nz,rg_valmanq_io,mask_h0, &
                              nz(:,limin:limax,ljmin:ljmax),0,0,l_out_nc4par = l_out_nc4par)
       CALL out_write_zxytvar(fileresu(nb),name_out_kz,rg_valmanq_io,mask_h0, &
                              kz(:,limin:limax,ljmin:ljmax),0,0,l_out_nc4par = l_out_nc4par)


     END IF



     IF (l_out_buoy(nb) .AND. .NOT.l_modele2d) THEN
       tmp3d(:,:,:) = 1000.0_rsh - 1000.0_rsh*boz(:,limin:limax,ljmin:ljmax)/grav
       CALL out_write_zxytvar(fileresu(nb),name_out_bz,rg_valmanq_io,mask_h0, &
                              tmp3d,0,1,l_out_nc4par = l_out_nc4par)
     END IF

     IF (l_out_wind(nb)) THEN
       CALL out_write_xytvar(fileresu(nb),name_out_tosx,rg_valmanq_io,mask_h0, &
                             tosx(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       CALL out_write_xytvar(fileresu(nb),name_out_tosy,rg_valmanq_io,mask_h0, &
                             tosy(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
     END IF





       IF (l_out_wave(nb)) THEN
         IF (l_wave_rdir) THEN
! trigonometric to oceanographic conventions : from north = 0, from west = 270
           tmp2d(limin:limax,ljmin:ljmax) = MOD(-wave_dir(limin:limax,ljmin:ljmax)*180.0_rsh/pi+270.0_rsh,360.0_rsh)
           CALL out_write_xytvar(fileresu(nb),name_out_wave_dir,rg_valmanq_io,mask_h0, &
                                 tmp2d,0,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_wave_rhs) &
         CALL out_write_xytvar(fileresu(nb),name_out_wave_hs,rg_valmanq_io,mask_h0, &
                               wave_hs(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         IF (l_wave_rfoc) &
         CALL out_write_xytvar(fileresu(nb),name_out_wave_foc,rg_valmanq_io,mask_h0, &
                               wave_phioc(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         IF (l_wave_rtauox .AND. l_wave_rtauoy) THEN
           CALL out_write_xytvar(fileresu(nb),name_out_wave_tauox,rg_valmanq_io,mask_h0, &
                               wave_tauox(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
           CALL out_write_xytvar(fileresu(nb),name_out_wave_tauoy,rg_valmanq_io,mask_h0, &
                               wave_tauoy(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_wave_rubot) &
         CALL out_write_xytvar(fileresu(nb),name_out_wave_ubot,rg_valmanq_io,mask_h0, &
                               wave_ubot(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         IF (l_wave_ruubot .AND. l_wave_rvubot) THEN
           CALL out_write_xytvar(fileresu(nb),name_out_wave_ubotx,rg_valmanq_io,mask_h0, &
                               wave_uubot(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
           CALL out_write_xytvar(fileresu(nb),name_out_wave_uboty,rg_valmanq_io,mask_h0, &
                               wave_vubot(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_wave_rt0m1) &
         CALL out_write_xytvar(fileresu(nb),name_out_wave_t0m1,rg_valmanq_io,mask_h0, &
                               wave_t0m1(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         IF (l_wave_rtauwx .AND. l_wave_rtauwy) THEN
           CALL out_write_xytvar(fileresu(nb),name_out_wave_tauwx,rg_valmanq_io,mask_h0, &
                               wave_tauwx(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
           CALL out_write_xytvar(fileresu(nb),name_out_wave_tauwy,rg_valmanq_io,mask_h0, &
                               wave_tauwy(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_wave_rfbb) &
         CALL out_write_xytvar(fileresu(nb),name_out_wave_fbb,rg_valmanq_io,mask_h0, &
                               wave_phibbl(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         IF (l_wave_rtbb) THEN
           CALL out_write_xytvar(fileresu(nb),name_out_wave_tbbx,rg_valmanq_io,mask_h0, &
                               wave_taubbl_x(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
           CALL out_write_xytvar(fileresu(nb),name_out_wave_tbby,rg_valmanq_io,mask_h0, &
                               wave_taubbl_y(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_wave_rtus) THEN
           CALL out_write_xytvar(fileresu(nb),name_out_wave_tusx,rg_valmanq_io,mask_h0, &
                               wave_tusx(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
           CALL out_write_xytvar(fileresu(nb),name_out_wave_tusy,rg_valmanq_io,mask_h0, &
                               wave_tusy(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_wave_russ) THEN
           CALL out_write_xytvar(fileresu(nb),name_out_wave_ussx,rg_valmanq_io,mask_h0, &
                               wave_ussx(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
           CALL out_write_xytvar(fileresu(nb),name_out_wave_ussy,rg_valmanq_io,mask_h0, &
                               wave_ussy(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)

         END IF
         IF (l_wave_rbhd) &
         CALL out_write_xytvar(fileresu(nb),name_out_wave_bhd,rg_valmanq_io,mask_h0, &
                               wave_bhd(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         IF (l_wave_rwlm) &
         CALL out_write_xytvar(fileresu(nb),name_out_wave_wlm,rg_valmanq_io,mask_h0, &
                               wave_wlm(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)

         IF (l_wave_rradstr) THEN
           CALL out_write_xytvar(fileresu(nb),name_out_wave_sx,rg_valmanq_io,mask_h0, &
                               wave_sx(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
           CALL out_write_xytvar(fileresu(nb),name_out_wave_sy,rg_valmanq_io,mask_h0, &
                               wave_sy(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         END IF
       END IF


    END IF   ! end test on .NOT.l_out_timeave(nb)

     IF (nbave /= 0 .AND. l_out_timeave(nb) ) THEN
         IF (l_out_ssh(nb)) CALL out_write_xytvar(fileresu(nb),name_out_xe,rg_valmanq_io,mask_h0, &
                                                ssh_ave(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)

         IF (l_out_ubt(nb)) THEN
           CALL out_write_xytvar(fileresu(nb),name_out_u,rg_valmanq_io,mask_hx, &
                               u_ave(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
           CALL out_write_xytvar(fileresu(nb),name_out_v,rg_valmanq_io,mask_hy, &
                               v_ave(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_out_u3d(nb) .AND. .NOT.l_modele2d) THEN
           CALL out_write_zxytvar(fileresu(nb),name_out_uz,rg_valmanq_io,mask_hx, &
                                uz_ave(nb,:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
           CALL out_write_zxytvar(fileresu(nb),name_out_vz,rg_valmanq_io,mask_hy, &
                                vz_ave(nb,:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)

         END IF

         IF (l_out_sal(nb) .AND. l_modele2d) THEN
           CALL out_write_xytvar(fileresu(nb),name_out_sal,rg_valmanq_io,mask_h0, &
                               sal_ave(nb,1,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         ELSE IF (l_out_sal(nb).AND..NOT.l_modele2d) THEN
           CALL out_write_zxytvar(fileresu(nb),name_out_sal,rg_valmanq_io,mask_h0, &
                                sal_ave(nb,:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_out_temp(nb) .AND. l_modele2d) THEN
           CALL out_write_xytvar(fileresu(nb),name_out_temp,rg_valmanq_io,mask_h0, &
                               temp_ave(nb,1,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         ELSE IF (l_out_temp(nb) .AND. .NOT.l_modele2d) THEN
           CALL out_write_zxytvar(fileresu(nb),name_out_temp,rg_valmanq_io,mask_h0, &
                                temp_ave(nb,:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
         END IF

         IF (l_out_visc(nb) .AND. l_modele2d) THEN
           CALL out_write_xytvar(fileresu(nb),name_out_vis2d,rg_valmanq_io,mask_h0, &
                               vis2d_ave(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
         ELSE IF (l_out_visc(nb) .AND. .NOT.l_modele2d) THEN
           CALL out_write_zxytvar(fileresu(nb),name_out_vishxe,rg_valmanq_io,mask_h0, &
                                vish_xe_ave(nb,:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
           CALL out_write_zxytvar(fileresu(nb),name_out_vishphi,rg_valmanq_io,mask_h0, &
                                vish_phi_ave(nb,:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
         END IF
         IF (l_out_turb(nb) .AND. .NOT.l_modele2d) THEN
           CALL out_write_zxytvar(fileresu(nb),name_out_nz,rg_valmanq_io,mask_h0, &
                                nz_ave(nb,:,limin:limax,ljmin:ljmax),0,0,l_out_nc4par = l_out_nc4par)
           CALL out_write_zxytvar(fileresu(nb),name_out_kz,rg_valmanq_io,mask_h0, &
                                kz_ave(nb,:,limin:limax,ljmin:ljmax),0,0,l_out_nc4par = l_out_nc4par)

         END IF

         IF (l_out_buoy(nb) .AND. .NOT.l_modele2d) THEN
           tmp3d(:,:,:) = 1000.0_rsh - 1000.0_rsh*boz_ave(nb,:,limin:limax,ljmin:ljmax)/grav
           CALL out_write_zxytvar(fileresu(nb),name_out_bz,rg_valmanq_io,mask_h0, &
                                tmp3d,0,1,l_out_nc4par = l_out_nc4par)
         END IF

! initialize averaged dynamical and thermodynamical variables
! -----------------------------------------------------------

         dtn = t - tsort(nb)
         IF (l_out_ssh(nb)) &
         ssh_ave(nb,:,:) = x(limin:limax,ljmin:ljmax)*dtn
         IF (l_out_ubt(nb)) THEN
           u_ave(nb,:,:) = u(limin:limax,ljmin:ljmax)*dtn
           v_ave(nb,:,:) = v(limin:limax,ljmin:ljmax)*dtn
         END IF
         IF (l_out_u3d(nb).AND..NOT.l_modele2d) THEN
           uz_ave(nb,:,:,:) = dtn * uz(:,limin:limax,ljmin:ljmax)
           vz_ave(nb,:,:,:) = dtn * vz(:,limin:limax,ljmin:ljmax)
         END IF
         IF (l_out_temp(nb)) &
           temp_ave(nb,:,:,:) = dtn * temp(:,limin:limax,ljmin:ljmax)
         IF (l_out_sal(nb)) &
           sal_ave(nb,:,:,:) = dtn * sal(:,limin:limax,ljmin:ljmax)
         IF (l_out_visc(nb).AND.l_modele2d) THEN
           vis2d_ave(nb,:,:) = dtn * vis2d(limin:limax,ljmin:ljmax)
         ELSE IF (l_out_visc(nb).AND..NOT.l_modele2d) THEN
           vish_xe_ave(nb,:,:,:) = dtn * vish_xe(:,limin:limax,ljmin:ljmax)
           vish_phi_ave(nb,:,:,:) = dtn * vish_phi(:,limin:limax,ljmin:ljmax)
         END IF
         IF (l_out_turb(nb).AND..NOT.l_modele2d) THEN
           nz_ave(nb,:,:,:) = dtn * nz(:,limin:limax,ljmin:ljmax)
           kz_ave(nb,:,:,:) = dtn * kz(:,limin:limax,ljmin:ljmax)
         END IF
         IF (l_out_buoy(nb).AND..NOT.l_modele2d) &
           boz_ave(nb,:,:,:) = dtn * boz(:,limin:limax,ljmin:ljmax)
         dt_ave(nb) =  dtn
     END IF

     IF (nbatm /= 0 .AND. l_out_atmsflx(nb) ) THEN
       CALL out_write_xytvar(fileresu(nb),name_out_mrad,rg_valmanq_io,mask_h0,     &
                             rad_cum(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       CALL out_write_xytvar(fileresu(nb),name_out_mirlos,rg_valmanq_io,mask_h0,   &
                             irlos_cum(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       CALL out_write_xytvar(fileresu(nb),name_out_mlatlos,rg_valmanq_io,mask_h0,  &
                             latlos_cum(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       CALL out_write_xytvar(fileresu(nb),name_out_msenslos,rg_valmanq_io,mask_h0, &
                             senslos_cum(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       CALL out_write_xytvar(fileresu(nb),name_out_mevap,rg_valmanq_io,mask_h0, &
                             evap_cum(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       CALL out_write_xytvar(fileresu(nb),name_out_mrain,rg_valmanq_io,mask_h0, &
                             rain_cum(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       CALL out_write_xytvar(fileresu(nb),name_out_mtosx,rg_valmanq_io,mask_h0, &
                             tosx_cum(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       CALL out_write_xytvar(fileresu(nb),name_out_mtosy,rg_valmanq_io,mask_h0, &
                             tosy_cum(nb,limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)


! initialize cumulative atmospherical surface fluxes
! --------------------------------------------------
 
       dtn = t - tsort(nb)
       rad_cum(nb,:,:)     = rad(limin:limax,ljmin:ljmax)    *dtn
       irlos_cum(nb,:,:)   = irlos(limin:limax,ljmin:ljmax)  *dtn
       latlos_cum(nb,:,:)  = latlos(limin:limax,ljmin:ljmax) *dtn
       senslos_cum(nb,:,:) = senslos(limin:limax,ljmin:ljmax)*dtn
       evap_cum(nb,:,:)    = evap(limin:limax,ljmin:ljmax)   *dtn
       rain_cum(nb,:,:)    = rain(limin:limax,ljmin:ljmax)   *dtn
       tosx_cum(nb,:,:)    = tosx(limin:limax,ljmin:ljmax)   *dtn
       tosy_cum(nb,:,:)    = tosy(limin:limax,ljmin:ljmax)   *dtn

       dt_cum(nb) =  dtn

     END IF




     DO l = 1,nv_state
      IF (l_out_subs(l)) THEN
       IF (l_modele2d) THEN
         tecrit(:,:)=rg_valmanq_io
         DO j=ljmin,ljmax
           DO i=MAX(limin,ig(j)),MIN(limax,id(j))
             IF ((j.GE.jb(i)).AND.(j.LE.jh(i))) THEN
               IF (l_out_log) THEN
                 tecrit(i,j) = LOG(MAX(cv_wat(l,1,i,j)/unit_modif_mudbio_N2dw(irk_fil(l)),1.0e-30))/LOG(10.0_rsh)
               ELSE
!tecrit(i,j) =MAX(cv_wat(l,1,i,j)/unit_modif_mudbio_N2dw(irk_fil(l)),1.0e-30)
                 tecrit(i,j) =cv_wat(l,1,i,j)/unit_modif_mudbio_N2dw(irk_fil(l))
               END IF
               IF (h0(i,j)+x(i,j) <= h0fond) tecrit(i,j)=rg_valmanq_io
             END IF
           END DO
         END DO
         CALL ionc4_write_xyt(fileresu(nb),TRIM(name_var(irk_fil(l))),tecrit,limin,limax,ljmin,ljmax,0,l_out_nc4par = l_out_nc4par)
       ELSE
         tzecrit(:,:,:)=rg_valmanq_io

! modif septembre 2014 : separation du test pour eviter depassement de tableaux sur l
           IF (l>=isand1 .AND. l<=isand2 ) THEN
             l_outsandrouse_l=l_outsandrouse(l)
           ELSE
             l_outsandrouse_l=.FALSE.
           ENDIF
           IF (l>=isand1 .AND. l<=isand2 .AND. l_outsandrouse_l .AND. l_subs2D(l)) THEN
           DO j=ljmin,ljmax
             DO i=MAX(limin,ig(j)),MIN(limax,id(j))
               IF ((j.GE.jb(i)).AND.(j.LE.jh(i))) THEN
                 IF (l_out_log) THEN
                   tzecrit(:,i,j) = LOG(MAX(cv_wat(l,:,i,j),1.0e-30))/LOG(10.0_rsh)
                 ELSE
                   htott=-f_zw(h0(i,j),x(i,j),0,i,j)+x(i,j)
                   IF (htott > h0fond .and. sum_tmp(l,i,j) > 0.0_rsh) THEN
                       DO k=1,kmax
                         zk=-f_zu(h0(i,j),x(i,j),k,i,j)+x(i,j)
                         IF (zk > 0.0_rsh) THEN
                           zk=htott-zk
                           tzecrit(k,i,j)=cv_wat(l,1,i,j)*htott*((htott-zk)/zk)**rouse2D(l,i,j)
                           tzecrit(k,i,j)=tzecrit(k,i,j)/sum_tmp(l,i,j)
                         ELSE
                           tzecrit(k,i,j)=0.0_rsh
                         ENDIF
                       ENDDO 
                   ELSE
                       tzecrit(:,i,j)=0.0_rsh
                   ENDIF
                 END IF
!MD 08/06/2018
!                 write(104,*) 'Rentre dans IF 3'
!                 IF (l_out_ws(nb)) THEN
!                   IF (l == imud1) THEN
!                     write(105,*) 'WS'
!                     DO k=1,kmax
!                       write(106,*) ws3(k,5,i,j)
!                       ws_mud(k,i,j)=ws3(k,5,i,j)
!                     END DO
!                   END IF
!                 END IF
!END MD 08/06/2018
                 IF (h0(i,j)+x(i,j) <= h0fond) tzecrit(:,i,j)=rg_valmanq_io
               END IF
             END DO
           END DO
         ELSE

          DO j=ljmin,ljmax
           DO i=MAX(limin,ig(j)),MIN(limax,id(j))
             IF ((j.GE.jb(i)).AND.(j.LE.jh(i))) THEN
               IF (l_out_log) THEN
                 tzecrit(:,i,j) = LOG(MAX(cv_wat(l,:,i,j)/unit_modif_mudbio_N2dw(irk_fil(l)),1.0e-30))/LOG(10.0_rsh)
               ELSE
                 tzecrit(:,i,j) = cv_wat(l,:,i,j)/unit_modif_mudbio_N2dw(irk_fil(l))
               END IF
               IF (h0(i,j)+x(i,j) <= h0fond) tzecrit(:,i,j)=rg_valmanq_io
             END IF
           END DO
          END DO

         ENDIF

         CALL ionc4_write_zxyt(fileresu(nb),TRIM(name_var(irk_fil(l))),tzecrit,limin,limax,ljmin,ljmax,1,kmax,0, &
                                            l_out_nc4par = l_out_nc4par)
!MD 08/06/2018
!         write(107,*) 'RENTRE DANS IF 4'
         IF (l_out_ws(nb)) THEN
             IF (l >= imud1 .AND. l <= imud2) THEN
               name_var_sf = TRIM('WS_'//name_var(irk_fil(l)))
               DO j=ljmin,ljmax
               DO i=MAX(limin,ig(j)),MIN(limax,id(j))
                 IF ((j.GE.jb(i)).AND.(j.LE.jh(i))) THEN
                   IF (l_out_log) THEN
                     tzecrit(:,i,j) = LOG(ws3(:,l,i,j))/LOG(10.0_rsh)
                   ELSE
                     tzecrit(:,i,j) = ws3(:,l,i,j)
                   END IF
                   IF (h0(i,j)+x(i,j) <= h0fond) tzecrit(:,i,j)=rg_valmanq_io
                 END IF
               END DO
               END DO
               CALL ionc4_write_zxyt(fileresu(nb),name_var_sf,tzecrit,limin,limax,ljmin,ljmax,1,kmax,0, &
                                             l_out_nc4par = l_out_nc4par)
             END IF
         END IF
         IF (l_out_varint(nb)) THEN
             tecrit(:,:)=rg_valmanq_io 
             DO j=ljmin,ljmax
                 DO i=MAX(limin,ig(j)),MIN(limax,id(j))
                     IF ((j.GE.jb(i)).AND.(j.LE.jh(i))) THEN
!                   DO i=limin,limax
                         tecrit(i,j)=SUM(cv_wat(l,:,i,j))
                         IF (h0(i,j)+x(i,j) <= h0fond) tecrit(i,j)=rg_valmanq_io
                     END IF
                 END DO
             END DO
             name_var_sf1 = TRIM('C_int_'//name_var(irk_fil(l)))
             CALL ionc4_write_xyt(fileresu(nb),name_var_sf1,tecrit,limin,limax,ljmin,ljmax,0,l_out_nc4par = l_out_nc4par)
         END IF
!END MD 08/06/2018
       END IF
      END IF
     END DO
     
     


     


! To write the data on the disk and not loose data in case of run crash
     CALL ionc4_sync(fileresu(nb))

   END IF

   END IF

   !PRINT*, 'exit OUTPUT_SAVE routine',TRIM(fileresu(nb))

  END SUBROUTINE output_save

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

