 MODULE output

!!======================================================================
!!                   ***  MODULE  output  ***
!! Outputs :  Initialyze type of outputs
!!            (global, zoom, section, precise date, filtered outputs ...)
!!            Save outputs
!!
!!   History :
!!        !  2011-01 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
!!
!!======================================================================





























































!! * Modules used
   USE parameters, ONLY : riosh,rsh,rlg,imin,imax,jmin,jmax,kmax,     &
                          rg_valmanq_io,dg_valmanq_io,lchain,         &
                          limin,limax,ljmin,ljmax

   USE parameters, ONLY : nb_var


   IMPLICIT NONE
   PRIVATE

!! * Accessibility
   PUBLIC output_alloc                ! routine called by alloc.F90
   PUBLIC output_dealloc              ! routine called by dealloc.F90
   PUBLIC output_init                 ! routine called by main.F90
   PUBLIC output_mng                  ! routine called by main.F90, step.F90

!! * Shared module variables
   LOGICAL, ALLOCATABLE, DIMENSION(:), PUBLIC :: l_out_ssh,l_out_sal,l_out_temp,   &
                                                 l_out_ubt,l_out_u3d,l_out_zaxis_w,&
                                                 l_out_visc,l_out_turb,l_out_buoy, &
                                                 l_out_atmsflx,l_out_wind,         &
                                                 l_out_wave,l_out_timeave
!MD 08/06/2018
   LOGICAL, ALLOCATABLE, DIMENSION(:), PUBLIC :: l_out_ws,l_out_varint
!END MD 08/06/2018
   LOGICAL, ALLOCATABLE, DIMENSION(:), PUBLIC :: l_out_demer,l_out_precdate, &
                                                 l_filebydate,l_out_off
   INTEGER, PUBLIC                            :: nb_out,nbatm,nbave
   INTEGER, ALLOCATABLE, DIMENSION(:), PUBLIC :: iextmin,iextmax,   &
                                                 jextmin,jextmax,   &
                                                 kextmin,kextmax,   &
                                                 stepi,stepj,stepk, &
                                                 igextmin,igextmax, &
                                                 jgextmin,jgextmax
   CHARACTER(LEN=lchain), PUBLIC  :: file_output  ! file defining outputs
   LOGICAL, PUBLIC                :: l_out_nc4par ! logical activating parallel I/O
   LOGICAL, PUBLIC                :: l_out_pack ! logical activating compression of real in short
   CHARACTER(LEN=lchain), ALLOCATABLE, DIMENSION(:), PUBLIC :: fileresu,fileresu_verif
   REAL(KIND=rlg), PUBLIC                     :: dtn               ! next period step
   REAL(KIND=rlg),        ALLOCATABLE, DIMENSION(:), PUBLIC :: tdebsor,tfinsor, &
                                                               tinter,tsort, &
                                                               pasor,dt_sum

! intermediate variables for time averaged outputs of thermo-dynamical variables
   REAL(KIND=rsh),ALLOCATABLE, DIMENSION(:,:,:), PUBLIC   :: ssh_ave,vis2d_ave,         &
                                                             u_ave,v_ave
   REAL(KIND=rsh),ALLOCATABLE, DIMENSION(:,:,:,:), PUBLIC :: uz_ave,vz_ave,         &
                                                             temp_ave,sal_ave,      &
                                                             vish_xe_ave,           &
                                                             vish_phi_ave,          &
                                                             nz_ave,kz_ave,boz_ave


! intermediate variables for outputs of cumulative atmospherical fluxes
   REAL(KIND=rsh),ALLOCATABLE, DIMENSION(:,:,:), PUBLIC   :: rad_cum,irlos_cum,         &
                                                             latlos_cum,senslos_cum,    &
                                                             evap_cum,rain_cum,         &
                                                             tosx_cum,tosy_cum

   REAL(KIND=rlg),ALLOCATABLE, DIMENSION(:), PUBLIC       :: dt_cum            ! period of integration
   REAL(KIND=rlg),ALLOCATABLE, DIMENSION(:), PUBLIC       :: dt_ave            ! period of the time average

!! * Private variables

! intermediate variables for outputs at precise time
! or for filtering on time
   REAL(KIND=rsh),ALLOCATABLE, DIMENSION(:,:)     :: ua,va,xea
   REAL(KIND=rsh),ALLOCATABLE, DIMENSION(:,:)     :: vis2da
   REAL(KIND=rsh),ALLOCATABLE, DIMENSION(:,:,:)   :: uza,vza,sala,tempa
   REAL(KIND=rsh),ALLOCATABLE, DIMENSION(:,:,:)   :: vish_xea,vish_phia,nza,kza,boza

!cval pourquoi garde ?   REAL(KIND=rsh),DIMENSION(nb_var,kmax,imin:imax,jmin:jmax) :: cv_wata
   REAL(KIND=rsh),ALLOCATABLE, DIMENSION(:,:,:,:) :: cv_wata

   REAL(KIND=rlg)                                 :: ta

! intermediate variables for outputs of cumulative atmospherical fluxes
   REAL(KIND=rlg)                                 :: dtp           ! previous period step
   REAL(KIND=rlg)                                 :: dtnext        ! next period step



 CONTAINS

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

  SUBROUTINE output_alloc

   USE comvars2d, ONLY : ierrorlog

   IF (ALLOCATED(ua)) WRITE(ierrorlog,*) 'ERROR: output_alloc already allocated'
   ALLOCATE(   ua(     limin:limax,ljmin:ljmax) )
   ALLOCATE(   va(     limin:limax,ljmin:ljmax) )
   ALLOCATE(  xea(     limin:limax,ljmin:ljmax) )
   ALLOCATE(vis2da(    limin:limax,ljmin:ljmax) )
   ALLOCATE( sala(kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(tempa(kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(vish_xea(kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(vish_phia(kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(nza(0:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(kza(0:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(boza(kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(  uza(kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(  vza(kmax,limin:limax,ljmin:ljmax) )
!cval pourquoi commente (lien ligne 61)

   ALLOCATE(cv_wata(nb_var,kmax,limin:limax,ljmin:ljmax) )



  END SUBROUTINE output_alloc

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

  SUBROUTINE output_dealloc
   DEALLOCATE(ua,va,xea,sala,tempa,uza,vza)
   DEALLOCATE(vish_xea,vish_phia,vis2da,nza,kza,boza)
!cval pourquoi commente (lien lignes 61 et 93)

   DEALLOCATE(cv_wata)


  END SUBROUTINE output_dealloc

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

  SUBROUTINE output_init

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE output_init  ***
!&E
!&E ** Purpose : lecture du fichier descriptif des sorties
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2007-01  (V. Garnier)  Original code
!&E       !  2009-05  (V. Garnier)  Save cumulative atmospherical fluxes
!&E       !  2009-11  (S. Michel, HOCER)  Check the content of the "file_output" file
!&E                                       and write error and warning messages
!&E       !  2010-02  (G. Simon) add l_out_wave
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE comvars2d, ONLY : iscreenlog,ierrorlog,iwarnlog,out_dir,suffix, &
                         tdeb,tfin,l_modele2d,l_initfromfile,lon2d_g,lat2d_g

   USE toolmpi,   ONLY : ADD_TAG,MASTER
   USE toolerror

!! * Arguments

!! * Local declarations
   LOGICAL                     :: ex,l_posit
   LOGICAL,DIMENSION(4)        :: lpara
   INTEGER                     :: eof,nb,nbdp,nboff,nberror,ispace,ntmpmin,ntmpmax
   INTEGER                     :: ncpos,ncdat,ncoff,ncpre,ncdem,ncssh,ncubt,ncwav
   INTEGER                     :: ncu3d,ncsal,nctem,ncwin,ncvis,nctur,ncbuo,ncflx,ncave
   INTEGER                     :: steptmp,stepmax
   INTEGER,DIMENSION(4)        :: ipara
   CHARACTER(LEN=5)            :: comment
   CHARACTER(LEN=19)           :: comdate,tool_sectodat
   CHARACTER(LEN=lchain)       :: charline
   REAL(KIND=rsh)              :: tool_latlon2i,tool_latlon2j
   REAL(KIND=riosh)            :: coordmin,coordmax,coord1,coord2
   REAL(KIND=rlg)              :: gwext,geext,fisext,finext,tdebtmp,tfintmp
   REAL(KIND=rlg),DIMENSION(4) :: dpara
!MD 08/06/2018
   INTEGER                     :: ncws
   INTEGER                     :: ncvarint
!END MD 08/06/2018

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

   !PRINT*, 'enter OUTPUT_INIT routine'

! save into simu.log
!-------------------
   IF (MASTER) THEN
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) '**************** OUTPUT_INIT.F90 *****************'
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) ' '
   ENDIF

! Initialize counters

   nb_out=0 ! number of output data files
   nbatm =0 ! output file number that will contain cumulated atmospheric fluxes
   nbave =0 ! output file number that will contain averaged dynamical and thermodynamical variables

   nbdp  =0 ! counter for output with Demerliac filter or precise date
   nboff =0 ! counter for saving conservative velocities for offline purpose (only once)

! test on description file existence
   INQUIRE(file=file_output,exist=ex)
   IF (.NOT.ex) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'ERROR: The output description file ',TRIM(file_output)
     WRITE(ierrorlog,*) 'does not exist.'
     WRITE(ierrorlog,*) 'No output file(s) can be written.'
     WRITE(ierrorlog,*) 'Simulation stopped.'
     STOP
   ENDIF

! read total number of diagnostic variables
   eof=0
   OPEN(49,file = file_output,form='formatted')

   comment='debut'
   DO WHILE (comment /= '*****')
     READ(49,'(a)',iostat=eof) comment
   ENDDO
   DO WHILE (eof==0)
     READ(49,'(a)',iostat=eof) comment
     IF ((comment == '*****').OR.(eof==1)) THEN
       nb_out=nb_out+1
       DO WHILE (((comment == '*****').OR.(comment == '     ')).AND.(eof==0))
         READ(49,'(a)',iostat=eof) comment
       END DO
     ENDIF
   ENDDO
   IF (MASTER) &
   WRITE(iscreenlog,*) 'Number of output files = ',nb_out
   IF (nb_out==0) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'ERROR: No defined output file in description file: ',TRIM(file_output)
     WRITE(ierrorlog,*) 'Simulation stopped.'
     STOP
   ENDIF

! allocate arrays for diagnostic variables
   ALLOCATE(fileresu(nb_out),iextmin(nb_out),iextmax(nb_out))

   ALLOCATE(fileresu_verif(nb_out),jextmin(nb_out),jextmax(nb_out))

   ALLOCATE(igextmin(nb_out),igextmax(nb_out))
   ALLOCATE(jgextmin(nb_out),jgextmax(nb_out))
   ALLOCATE(kextmin(nb_out),kextmax(nb_out))
   ALLOCATE(stepi(nb_out),stepj(nb_out),stepk(nb_out))
   ALLOCATE(l_out_ssh(nb_out),l_out_sal(nb_out),l_out_temp(nb_out))
   ALLOCATE(l_out_ubt(nb_out),l_out_u3d(nb_out),l_out_wind(nb_out))
   ALLOCATE(l_out_wave(nb_out),l_out_zaxis_w(nb_out))
   ALLOCATE(l_out_visc(nb_out),l_out_turb(nb_out),l_out_buoy(nb_out),l_out_atmsflx(nb_out))
   ALLOCATE(l_out_timeave(nb_out))
   ALLOCATE(l_out_demer(nb_out),l_out_precdate(nb_out),l_filebydate(nb_out),l_out_off(nb_out))
   ALLOCATE(tdebsor(nb_out),tfinsor(nb_out),tinter(nb_out))
   ALLOCATE(tsort(nb_out),pasor(nb_out),dt_sum(nb_out))
!MD 08/06/2018
   ALLOCATE(l_out_ws(nb_out))
   ALLOCATE(l_out_varint(nb_out))
!END MD 08/06/2018

! initialize cumulative time step
   dt_sum(:)=0.0_rlg
! initialize real parameters
   tdebsor(:)=0.0_rlg;tfinsor(:)=0.0_rlg;pasor(:)=0.0_rlg
! initialize integer parameters
   iextmin(:)=-1;iextmax(:)=imax+1
   jextmin(:)=-1;jextmax(:)=jmax+1
   kextmin(:)=0 ;kextmax(:)=kmax+1
   stepi(:)=0;stepj(:)=0;stepk(:)=0
! initialize boolean parameters
   l_out_ssh(:)  =.FALSE.;l_out_sal(:)     =.FALSE.;l_out_temp(:)   =.FALSE.
   l_out_ubt(:)  =.FALSE.;l_out_u3d(:)     =.FALSE.;l_out_wind(:)   =.FALSE.
   l_out_visc(:) =.FALSE.;l_out_turb(:)    =.FALSE.;l_out_atmsflx(:)=.FALSE.
   l_out_timeave(:)=.FALSE.
   l_out_demer(:)=.FALSE.;l_out_precdate(:)=.FALSE.;l_filebydate(:) =.FALSE.
   l_out_off(:)  =.FALSE.;l_out_buoy(:)    =.FALSE.;l_out_wave(:)   =.FALSE.
   l_out_zaxis_w(:)=.FALSE.
!MD 08/06/2018
   l_out_ws(:)=.FALSE.
   l_out_varint(:)=.FALSE.
!END MD 06/08/2018



! skip comment lines
   REWIND(49)
   comment='debut'
   DO WHILE (comment /= '*****')
     READ(49,'(a)',iostat=eof) comment
   END DO

! beginning of loop on the output file number
   DO nb=1,nb_out

! Initialize the error counter.
     nberror=0

! Initialize the horizontal extension.
     l_posit=.TRUE.
     gwext=-999.9_rlg;geext=-999.9_rlg
     fisext=-999.9_rlg;finext=-999.9_rlg

! Initialize the boolean counters.
     ncpos=0;ncdat=0;ncoff=0;ncpre=0;ncdem=0;ncssh=0;ncubt=0;ncwav=0
     ncu3d=0;ncsal=0;nctem=0;ncwin=0;ncvis=0;nctur=0;ncbuo=0;ncflx=0
     ncave=0
!MD 08/06/2018
     ncws=0
     ncvarint=0
!END MD 08/06/2018

! output file name
     READ(49,'(a)',iostat=eof) charline
     charline=ADJUSTL(charline)
     ispace=INDEX(charline,' ')
     fileresu(nb)=charline(1:ispace-1)

! read parameters for the current output file
     DO WHILE ((charline(1:5) /= '*****').AND.(eof==0))
       READ(49,'(a)',iostat=eof) charline

! start date for saving data
       IF (INDEX(charline,'date_startout')>1) THEN
         CALL date_check(tdebsor(nb),charline,nberror,file_output,nb)

! end date for saving data
       ELSE IF (INDEX(charline,'date_endout')>1) THEN
         CALL date_check(tfinsor(nb),charline,nberror,file_output,nb)

! time step for saving data
       ELSE IF (INDEX(charline,'pasor')>1) THEN
         CALL para_check(dpara,1,charline,nberror,file_output,nb)
         pasor(nb)=dpara(1)*3600.0_rlg

! vertical extension limits
       ELSE IF (INDEX(charline,'kextmin')>1) THEN
         CALL para_check(ipara,2,charline,nberror,file_output,nb)
         kextmin(nb)=ipara(1);kextmax(nb)=ipara(2);

! horizontal extension with indices or positions
       ELSE IF (INDEX(charline,'posit')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_posit=lpara(1);ncpos=ncpos+1

! horizontal extension limits
       ELSE IF (INDEX(charline,'iextmin')>1) THEN
         IF (l_posit) THEN
           CALL para_check(ipara,4,charline,nberror,file_output,nb)
           iextmin(nb)=ipara(1);iextmax(nb)=ipara(2)
           jextmin(nb)=ipara(3);jextmax(nb)=ipara(4)
         ELSE
           CALL para_check(dpara,4,charline,nberror,file_output,nb)
           gwext=dpara(1);geext=dpara(2);fisext=dpara(3);finext=dpara(4)
         ENDIF

! x-step, y-step, k-step for saving data
       ELSE IF (INDEX(charline,'step')>1) THEN
         CALL para_check(ipara,3,charline,nberror,file_output,nb)
         stepi(nb)=ipara(1);stepj(nb)=ipara(2);stepk(nb)=ipara(3)

! read logical parameters (booleans)
       ELSE IF (INDEX(charline,'filebydate')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_filebydate(nb)=lpara(1);ncdat=ncdat+1
       ELSE IF (INDEX(charline,'off')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_off(nb)=lpara(1);ncoff=ncoff+1
         IF (l_out_off(nb)) nboff=nboff+1
       ELSE IF (INDEX(charline,'precdate')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_precdate(nb)=lpara(1);ncpre=ncpre+1
         IF (l_out_precdate(nb)) nbdp = nbdp+1
       ELSE IF (INDEX(charline,'demer')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_demer(nb)=lpara(1);ncdem=ncdem+1
         IF (l_out_demer(nb)) nbdp=nbdp+1
       ELSE IF (INDEX(charline,'ssh')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_ssh(nb)=lpara(1);ncssh=ncssh+1
       ELSE IF (INDEX(charline,'ubt')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_ubt(nb)=lpara(1);ncubt=ncubt+1
       ELSE IF (INDEX(charline,'u3d')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_u3d(nb)=lpara(1);ncu3d=ncu3d+1
!MD 08/06/2018
       ELSE IF (INDEX(charline,'ws')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_ws(nb)=lpara(1);ncws=ncws+1
       ELSE IF (INDEX(charline,'varint')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_varint(nb)=lpara(1);ncvarint=ncvarint+1
!END MD 08/06/2018
       ELSE IF (INDEX(charline,'sal')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_sal(nb)=lpara(1);ncsal=ncsal+1
       ELSE IF (INDEX(charline,'temp')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_temp(nb)=lpara(1);nctem=nctem+1
       ELSE IF (INDEX(charline,'wind')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_wind(nb)=lpara(1);ncwin=ncwin+1
       ELSE IF (INDEX(charline,'wave')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_wave(nb)=lpara(1);ncwav=ncwav+1
       ELSE IF (INDEX(charline,'visc')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_visc(nb)=lpara(1);ncvis=ncvis+1
       ELSE IF (INDEX(charline,'turb')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_turb(nb)=lpara(1);nctur=nctur+1
       ELSE IF (INDEX(charline,'buoy')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_buoy(nb)=lpara(1);ncbuo=ncbuo+1
       ELSE IF (INDEX(charline,'flx')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_atmsflx(nb)=lpara(1);ncflx=ncflx+1
         IF (l_out_atmsflx(nb)) nbatm=nb
       ELSE IF (INDEX(charline,'ave')>1) THEN
         CALL para_check(lpara,1,charline,nberror,file_output,nb)
         l_out_timeave(nb)=lpara(1);ncave=ncave+1
         IF (l_out_timeave(nb)) nbave=nb

       ELSE IF ((charline(1:5)/='*****').AND.(charline(1:5)/='     ')) THEN
         WRITE(iwarnlog,*) ' '
         WRITE(iwarnlog,*) 'WARNING: In file ',TRIM(file_output),', file number : ',nb
         WRITE(iwarnlog,*) 'unrecognized descriptor at the end of this line :'
         WRITE(iwarnlog,*) TRIM(charline)
         WRITE(iwarnlog,*) 'This line is ignored.'
       ENDIF

! end of readding parameters for the current output file
     ENDDO



     IF ( l_out_off(nb) ) THEN
       WRITE(ierrorlog,*) ' '
       WRITE(ierrorlog,*) 'ERROR: In ',TRIM(file_output),' file, file number : ',nb
       WRITE(ierrorlog,*) '"l_out_off" cannot be set to TRUE, because the CPP key '
       WRITE(ierrorlog,*) '"key_offline_outuvflux" has not been selected.'
       nberror=nberror+1
     END IF


     IF (nberror>0) THEN
       WRITE(ierrorlog,*) 'Compare your "output.dat" parameter file to the examples'
       WRITE(ierrorlog,*) 'in the directory $HOMEMARS/../EXAMPLES/'
       WRITE(ierrorlog,*) 'Simulation stopped.'
       STOP
     ENDIF

! Concatenate the ouput filename.
     IF ( .NOT. l_filebydate(nb) ) THEN
       fileresu(nb)=TRIM(out_dir)//TRIM(fileresu(nb))//'_'//TRIM(suffix)//'.nc'

       IF (.NOT. l_out_nc4par) CALL ADD_TAG(fileresu(nb),LEN(fileresu(nb)))

       fileresu_verif(nb)=fileresu(nb)
     END IF

! Check that dates of data saving are defined.
     IF ((tdebsor(nb)==0.0_rlg).OR.(tfinsor(nb)==0.0_rlg).OR.(pasor(nb)==0.0_rlg)) THEN
       WRITE(ierrorlog,*) ' '
       WRITE(ierrorlog,*) 'ERROR: In file ',TRIM(file_output) ,', file number : ',nb
       WRITE(ierrorlog,*) 'the beginning date, ending date and time-step for saving data'
       WRITE(ierrorlog,*) 'must be defined.'
       nberror=nberror+1
     ENDIF

! Check that the output start and stop dates are correct.
     tdebtmp=tdebsor(nb);tfintmp=tfinsor(nb)
     IF (.NOT.l_initfromfile) THEN
       IF (l_modele2d) THEN
         CALL range_check(tdeb,tdebtmp,tfintmp,tfin,'out',file_output,nb)
       ELSE
         CALL range_check(tdeb,tdebtmp,tfintmp,tfin,'out',file_output,nb)
       ENDIF
     ELSE
       CALL range_check(tdebtmp,tdebtmp,tfintmp,tfin,'out',file_output,nb)
     ENDIF
     tdebsor(nb)=tdebtmp;tfinsor(nb)=tfintmp

! Check the output time-step.
     IF (pasor(nb)>tfinsor(nb)-tdebsor(nb)) THEN
       WRITE(iwarnlog,*) ' '
       WRITE(iwarnlog,*) 'WARNING: In file ',TRIM(file_output),', file number : ',nb
       WRITE(iwarnlog,*) 'the output time-step is longer than the output interval:'
       WRITE(iwarnlog,*) 'pasor = ',pasor(nb)/3600.0_rlg,&
         ' hours < tfinsor-tdebsor = ',(tfinsor(nb)-tdebsor(nb))/3600.0_rlg,' hours'
       WRITE(iwarnlog,*) 'The time-step is set to the interval duration.'
       pasor(nb)=tfinsor(nb)-tdebsor(nb)
     ENDIF

! Check vertical extension limits
     ntmpmin=kextmin(nb);ntmpmax=kextmax(nb)
     CALL range_check(1,ntmpmin,ntmpmax,kmax,'kext',file_output,nb)
     kextmin(nb)=ntmpmin;kextmax(nb)=ntmpmax

! Check horizontal extension limits
     IF (l_posit) THEN
! horizontal extent defined by i/j indices
       ntmpmin=iextmin(nb);ntmpmax=iextmax(nb)
       CALL range_check(imin,ntmpmin,ntmpmax,imax,'iext',file_output,nb)
       iextmin(nb)=ntmpmin;iextmax(nb)=ntmpmax

       ntmpmin=jextmin(nb);ntmpmax=jextmax(nb)
       CALL range_check(jmin,ntmpmin,ntmpmax,jmax,'jext',file_output,nb)
       jextmin(nb)=ntmpmin;jextmax(nb)=ntmpmax

     ELSE
       IF ((gwext ==-999.9_rlg).OR.(geext ==-999.9_rlg).OR. &
           (fisext==-999.9_rlg).OR.(finext==-999.9_rlg)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: In file ',TRIM(file_output) ,', file number : ',nb
         WRITE(ierrorlog,*) 'with l_posit=FALSE, the longitudes and latitudes of the saving domain'
         WRITE(ierrorlog,*) 'must be defined.'
         nberror=nberror+1
       ELSEIF ((ABS(gwext)>180.0_rlg).OR.(ABS(geext)>180.0_rlg).OR. &
         (ABS(fisext)>90.0_rlg).OR.(ABS(finext)>90.0_rlg).OR. &
         (gwext>geext).OR.(fisext>finext)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: In file ',TRIM(file_output) ,', file number : ',nb
         WRITE(ierrorlog,*) 'the longitudes and latitudes of the saving domain are out-of-bounds:'
         WRITE(ierrorlog,*) 'gwext  = ',gwext,', geext  = ',geext
         WRITE(ierrorlog,*) 'fisext = ',fisext,', finext = ',finext
         nberror=nberror+1
       ELSE
! horizontal extent defined by longitudes/latitudes
         coordmin=REAL(lon2d_g(imin,jmin),riosh)
         coordmax=REAL(lon2d_g(imax,jmax),riosh)
         coord1=REAL(gwext,riosh);coord2=REAL(geext,riosh)
         CALL range_check(coordmin,coord1,coord2,coordmax,'lonext',file_output,nb)
         gwext=REAL(coord1,rlg);geext=REAL(coord2,rlg)

         coordmin=REAL(lat2d_g(imin,jmin),riosh)
         coordmax=REAL(lat2d_g(imax,jmax),riosh)
         coord1=REAL(fisext,riosh);coord2=REAL(finext,riosh)
         CALL range_check(coordmin,coord1,coord2,coordmax,'latext',file_output,nb)
         fisext=REAL(coord1,rlg);finext=REAL(coord2,rlg)

         iextmin(nb)=NINT(tool_latlon2i(gwext,fisext))
         iextmax(nb)=NINT(tool_latlon2i(geext,finext))
         jextmin(nb)=NINT(tool_latlon2j(gwext,fisext))
         jextmax(nb)=NINT(tool_latlon2j(geext,finext))
       ENDIF
     ENDIF

! Check i-step, j-step, k-step for saving
     steptmp=stepi(nb);stepmax=iextmax(nb)-iextmin(nb)+1
     CALL step_check(1,steptmp,stepmax,'istep',file_output,nb)
     stepi(nb)=steptmp

     steptmp=stepj(nb);stepmax=jextmax(nb)-jextmin(nb)+1
     CALL step_check(1,steptmp,stepmax,'jstep',file_output,nb)
     stepj(nb)=steptmp

     steptmp=stepk(nb);stepmax=kextmax(nb)-kextmin(nb)+1
     CALL step_check(1,steptmp,stepmax,'kstep',file_output,nb)
     stepk(nb)=steptmp

! global spatial extension
     igextmin(nb)=iextmin(nb);igextmax(nb)=iextmax(nb)
     jgextmin(nb)=jextmin(nb);jgextmax(nb)=jextmax(nb)


! correct j indices
     DO WHILE ( jextmin(nb) < ljmin )
       jextmin(nb)=jextmin(nb)+stepj(nb)
     END DO
     jextmax(nb)=MIN(ljmax,jextmax(nb))



! correct i indices
     DO WHILE ( iextmin(nb) < limin )
       iextmin(nb)=iextmin(nb)+stepi(nb)
     END DO
     iextmax(nb)=MIN(limax,iextmax(nb))




! Check boolean parameters
     CALL para_warn(l_posit           ,'l_posit'       ,ncpos,file_output,nb)
     CALL para_warn(l_filebydate(nb)  ,'l_filebydate'  ,ncdat,file_output,nb)
     CALL para_warn(l_out_off(nb)     ,'l_out_off'     ,ncoff,file_output,nb)
     CALL para_warn(l_out_precdate(nb),'l_out_precdate',ncpre,file_output,nb)
     CALL para_warn(l_out_demer(nb)   ,'l_out_demer'   ,ncdem,file_output,nb)
     CALL para_warn(l_out_ssh(nb)     ,'l_out_ssh'     ,ncssh,file_output,nb)
     CALL para_warn(l_out_ubt(nb)     ,'l_out_ubt'     ,ncubt,file_output,nb)
     CALL para_warn(l_out_u3d(nb)     ,'l_out_u3d'     ,ncu3d,file_output,nb)
     CALL para_warn(l_out_sal(nb)     ,'l_out_sal'     ,ncsal,file_output,nb)
     CALL para_warn(l_out_temp(nb)    ,'l_out_temp'    ,nctem,file_output,nb)
     CALL para_warn(l_out_wind(nb)    ,'l_out_wind'    ,ncwin,file_output,nb)
     CALL para_warn(l_out_wave(nb)    ,'l_out_wave'    ,ncwav,file_output,nb)
     CALL para_warn(l_out_visc(nb)    ,'l_out_visc'    ,ncvis,file_output,nb)
     CALL para_warn(l_out_turb(nb)    ,'l_out_turb'    ,nctur,file_output,nb)
     CALL para_warn(l_out_buoy(nb)    ,'l_out_buoy'    ,ncbuo,file_output,nb)
     CALL para_warn(l_out_atmsflx(nb) ,'l_out_atmsflx' ,ncflx,file_output,nb)
     CALL para_warn(l_out_timeave(nb) ,'l_out_timeave' ,ncave,file_output,nb)
!MD 08/06/2018
!     write(109,*) 'FAIT LE CALL DU OUTPUT'
     CALL para_warn(l_out_ws(nb)      ,'l_out_ws'      ,ncws,file_output,nb)
     CALL para_warn(l_out_varint(nb)  ,'l_out_varint'  ,ncvarint,file_output,nb)
!END MD 08/06/2018

! Check dependencies between several boolean parameters.

! Check barotropic and/or baroclinic velocities saving.
     IF ((.NOT.l_modele2d).AND.(l_out_ubt(nb).AND.(.NOT.l_out_u3d(nb)))) THEN
       WRITE(iwarnlog,*) ' '
       WRITE(iwarnlog,*) 'WARNING: In file ',TRIM(file_output),', file number : ',nb
       WRITE(iwarnlog,*) 'the 3D velocity is not saved while the barotropic velocity is saved.'
     ENDIF

     IF ((l_modele2d).AND.(l_out_u3d(nb))) THEN
       l_out_u3d(nb)=.FALSE.
       WRITE(iwarnlog,*) 'WARNING: In file ',TRIM(file_output),', file number : ',nb
       WRITE(iwarnlog,*) 'the 3D velocity can not be saved in a 2D simulation.'
       WRITE(iwarnlog,*) 'l_out_u3d = ',l_out_u3d(nb)
     ENDIF

! Save 3D velocities in conservative flux form (offline purpose).
     IF ((l_modele2d).AND.(l_out_off(nb))) THEN
       l_out_off(nb)=.FALSE.
       WRITE(iwarnlog,*) 'WARNING: In file ',TRIM(file_output),', file number : ',nb
       WRITE(iwarnlog,*) 'conservative fluxes can only be computed in a 3D simulation.'
       WRITE(iwarnlog,*) 'l_out_off = ',l_out_off(nb)
     ENDIF

     IF (l_out_off(nb)) THEN
       IF (.NOT.l_out_u3d(nb)) THEN
         l_out_u3d(nb)=.TRUE.
         WRITE(iwarnlog,*) 'WARNING: In file ',TRIM(file_output),', file number : ',nb
         WRITE(iwarnlog,*) 'the 3D velocity must be saved to compute conservative fluxes.'
         WRITE(iwarnlog,*) 'l_out_u3d = ',l_out_u3d(nb)
       ENDIF
       IF (.NOT.l_filebydate(nb)) THEN
         l_filebydate(nb)=.TRUE.    ! do not change
         WRITE(iwarnlog,*) 'WARNING: In file ',TRIM(file_output),', file number : ',nb
         WRITE(iwarnlog,*) 'the outputs must be written in a file for each date'
         WRITE(iwarnlog,*) 'to save 3D velocities in conservative flux form.'
         WRITE(iwarnlog,*) 'l_filebydate = ',l_filebydate(nb)
       ENDIF
       IF ((l_out_demer(nb)).OR.(l_out_precdate(nb))) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: In file ',TRIM(file_output),', file number : ',nb
         WRITE(ierrorlog,*) 'Saving 3D velocities in conservative flux form (l_out_off=TRUE)'
         WRITE(ierrorlog,*) 'is not compatible with l_out_demer=TRUE or l_out_precdate=TRUE.'
         nberror=nberror+1
       ENDIF
     ENDIF

     IF (l_out_atmsflx(nb)) THEN
       IF ((l_out_demer(nb)).OR.(l_out_precdate(nb))) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: In file ',TRIM(file_output),', file number : ',nb
         WRITE(ierrorlog,*) 'Saving cumulative atmospherical fluxes'
         WRITE(ierrorlog,*) 'is not compatible with l_out_demer=TRUE or l_out_precdate=TRUE.'
         nberror=nberror+1
       ENDIF
     ENDIF

! Prevent saving the mixing coefficients when the Demerliac filter is activated.
     IF (((l_out_visc(nb)).OR.(l_out_turb(nb))).AND.(l_out_demer(nb))) THEN
       l_out_visc(nb)=.FALSE.
       l_out_turb(nb)=.FALSE.
       WRITE(iwarnlog,*) ' '
       WRITE(iwarnlog,*) 'WARNING: In ',TRIM(file_output),', file number : ',nb
       WRITE(iwarnlog,*) 'Filtering the viscosity or turbulence coefficients would make them'
       WRITE(iwarnlog,*) 'useless for physical interpretation, because large variations'
       WRITE(iwarnlog,*) 'in mixing  have an irreversible influence on stratification'
       WRITE(iwarnlog,*) 'Thus, in ',TRIM(file_output),', file number : ',nb
       WRITE(iwarnlog,*) 'l_out_visc = ',l_out_visc(nb)
       WRITE(iwarnlog,*) 'l_out_turb = ',l_out_turb(nb)
     END IF

     IF (nberror>0) THEN
       WRITE(ierrorlog,*) 'Compare your "output.dat" parameter file to the examples'
       WRITE(ierrorlog,*) 'in the directory $HOMEMARS/../EXAMPLES/'
       WRITE(ierrorlog,*) 'Simulation stopped.'
       STOP
     ENDIF

! end of loop on the output file number
   END DO

! Make sure user saves 3D velocities in conservative flux form (offline purpose) only once.
   IF (nboff>1) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'ERROR: In file : ',TRIM(file_output)
     WRITE(ierrorlog,*) 'the parameter l_out_off is set to TRUE ', nboff,' times,'
     WRITE(ierrorlog,*) 'but only one time is allowed.'
     WRITE(ierrorlog,*) 'Simulation stopped.'
     STOP
   END IF

! Make sure user calls demerliac or sorprec only once.
   IF (nbdp>1) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'ERROR: In file : ',TRIM(file_output)
     WRITE(ierrorlog,*) 'the parameter l_out_precdate or l_out_demer is set to TRUE ', nbdp,' times,'
     WRITE(ierrorlog,*) 'but only one time is allowed.'
     WRITE(ierrorlog,*) 'Simulation stopped.'
     STOP
   END IF

! save into simu.log
!-------------------
   IF (MASTER) THEN
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) 'OUTPUTS '
   ENDIF
   DO nb = 1, nb_out
     IF (MASTER) THEN
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) 'output file number : ',nb
     WRITE(iscreenlog,*) 'file name : ', TRIM(fileresu(nb))
     ENDIF
     comdate=tool_sectodat(tdebsor(nb))
     IF (MASTER) &
     WRITE(iscreenlog,*) 'starting date : ', comdate
     comdate=tool_sectodat(tfinsor(nb))
     IF (MASTER) THEN
     WRITE(iscreenlog,*) 'ending date : ', comdate
     WRITE(iscreenlog,*) 'time step (hours) : ', pasor(nb)/3600.0_rlg
     WRITE(iscreenlog,*) 'save from ',igextmin(nb),' to ',igextmax(nb),'every ',stepi(nb),' i'
     WRITE(iscreenlog,*) 'save from ',jgextmin(nb),' to ',jgextmax(nb),'every ',stepj(nb),' j'
     WRITE(iscreenlog,*) 'save from ',kextmin(nb),' to ',kextmax(nb),'every ',stepk(nb),' k'
     WRITE(iscreenlog,*) 'a single file created at each date ? ',l_filebydate(nb)
     WRITE(iscreenlog,*) 'save 3D velocities in conservative flux form (offline purpose) ? ',&
                        l_out_off(nb)
     WRITE(iscreenlog,*) 'save data at precise times ? ',l_out_precdate(nb)
     WRITE(iscreenlog,*) 'filter data outputs ? ',l_out_demer(nb)
     WRITE(iscreenlog,*) 'save sea surface height ? ',l_out_ssh(nb)
     WRITE(iscreenlog,*) 'save 2D barotropic velocities ? ', l_out_ubt(nb)
     WRITE(iscreenlog,*) 'save 3D baroclinic velocities ? ',l_out_u3d(nb)
     WRITE(iscreenlog,*) 'save salinity ? ',l_out_sal(nb)
     WRITE(iscreenlog,*) 'save temperature ? ',l_out_temp(nb)
     WRITE(iscreenlog,*) 'save wind stresses ? ',l_out_wind(nb)
     WRITE(iscreenlog,*) 'save wave forcing ? ',l_out_wave(nb)
     WRITE(iscreenlog,*) 'save viscosity coefficients ? ',l_out_visc(nb)
     WRITE(iscreenlog,*) 'save turbulence coefficients ? ',l_out_turb(nb)
     WRITE(iscreenlog,*) 'save buoyancy ? ',l_out_buoy(nb)
     WRITE(iscreenlog,*) 'save cumulative atmospheric fluxes ? ',l_out_atmsflx(nb)
     WRITE(iscreenlog,*) 'save averaged dynamical and thermodynamical variables ? ',l_out_timeave(nb)
!MD 08/06/2018
     WRITE(iscreenlog,*) 'save 3d settling velocity ? ',l_out_ws(nb)
     WRITE(iscreenlog,*) 'save depth integrated concentrations ? ',l_out_varint(nb)
!END MD 08/06/2018
     WRITE(iscreenlog,*) ' '
     ENDIF


     IF (l_out_turb(nb)) l_out_zaxis_w(nb)=.true.


   END DO

   !PRINT*, 'exit OUTPUT_INIT routine'

  END SUBROUTINE output_init

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

  SUBROUTINE output_mng

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE output_mng  ***
!&E
!&E ** Purpose : gestion des sorties de resultats, filtrage demerliac
!&E              sorties precises ...
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !                             Original code
!&E       !  2004-08 (F. Dumas)
!&E       !  2006-03 (V. Garnier) Dates management according to init_fromfile...
!&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 (flux meanning) variable uz and vz (for offline purpose)
!&E       !  2009-05 (V. Garnier) Save cumulative atmospherical fluxes
!&E       !  2009-09 (S. Michel, HOCER) add the arguments dimensions in the call to out_demerliac
!&E       !  2009-09 (V. Garnier) Shortening of the module
!&E       !  2009-09 (P. Garreau) Saving of viscosity and diffusivity in output files
!&E       !  2010-03 (V. Garnier) tsort(nb) = MAX(t+halfdt,t+pasor(nb)) to
!&E       !                       allow correct cumulative flux even for small output time steps
!&E       !  2014-07 (F. Ganthy)  Add outputs for obstructions module
!&E       !  2017-09 (F. Ganthy)  Add output for turbulence advection (k-epsilon)
!&E
!&E---------------------------------------------------------------------
!! Modules used
   USE comvars2d   , ONLY : nbouc,t,date_startshift,tdeb,dt,dtmin,l_out_subs_diag
   USE comvars2d   , ONLY : l_initfromfile,l_modele2d,halfdt,t
   USE comvars2d   , ONLY : hex,hey,ierrorlog
   USE comvars2d   , ONLY : vis2d,vish_phi,vish_xe
   USE comvars3d   , ONLY : kz,nz

   USE comvarp2d   , ONLY : u,v,ssh,sal,temp
   USE comvarp3d   , ONLY : uz,vz,boz
   USE comdem      , ONLY : xed,ud,uzd,vd,vzd,sald,tempd,nsdem
   USE sflxatm     , ONLY : rad,irlos,latlos,senslos

   USE comfluxdo   , ONLY : evap,rain
   USE cometdyn    , ONLY : tosx,tosy

   USE comsubstance, ONLY : cv_wat,nv_state
!   USE output      , only:cv_wata
   USE comdem      , ONLY : cv_watd


   USE comsubstance, ONLY : tdebsubs



   USE comvarp3d, ONLY : ect
   USE comturb,   ONLY : eps

!! * Arguments

!! * Local declarations
   INTEGER           :: l,nb

   REAL(KIND=rsh),DIMENSION(kmax,limin:limax,ljmin:ljmax) :: c1,c2

   REAL(KIND=rlg)    :: tect


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

   !PRINT*, 'enter OUTPUT_MNG routine'

! INITIALIZATION
! --------------

   IF (nbouc == 0) THEN

! Loop on number of output file
    DO nb = 1,nb_out

! Initialize dates

!     IF (l_modele2d .OR. .NOT.l_initfromfile) THEN
       IF ( (l_out_precdate(nb) .OR. l_out_demer(nb) .OR. l_out_atmsflx(nb)) .AND. tdebsor(nb) <= t) THEN
         tsort(nb)   =MAX(tdebsor(nb),t)+pasor(nb)
         tinter(nb)  =MAX(tdebsor(nb),t)+pasor(nb)
         tdebsor(nb) =MAX(tdebsor(nb),t)+pasor(nb)
       ELSE
         tdebsor(nb) =MAX(tdebsor(nb),t)
         tsort(nb)   =MAX(tdebsor(nb),t)
         tinter(nb)  =MAX(tdebsor(nb),t)
       END IF
!     ELSE
!       tsort(nb)   =MAX(tdebsor(nb),tdeb)
!       tinter(nb)  =MAX(tdebsor(nb),tdeb)
!       tdebsor(nb) =MAX(tdebsor(nb),tdeb)
!     END IF

     IF (l_out_demer(nb)) THEN
       nsdem=0
       pasor(nb)=MAX(pasor(nb),72.0_rlg*3600.0_rlg)
     END IF

     IF (date_startshift == 0.0_rlg) THEN
!       tsort(nb)=tsort(nb)+pasor(nb)
       IF (l_out_precdate(nb) .OR. l_out_demer(nb) ) THEN
!tinter(nb)=tinter(nb)+pasor(nb)
         IF (pasor(nb) < dtmin*1.5_rlg) THEN
           WRITE(ierrorlog,*) 'ERROR: ouput number : ',nb
           WRITE(ierrorlog,*) 'Variable "pasor" ',pasor(nb),' is not large enough'
           WRITE(ierrorlog,*) 'while using l_out_precdate or l_out_demer'
           WRITE(ierrorlog,*) 'Simulation stopped.'
           STOP
         END IF
       END IF
     END IF




! Create obstruction output file
!#ifdef key_obstructions
!     IF (nb==1 .AND. .NOT. l_filebydate(nb)) CALL obst_output(nb,t,ssh(limin:limax,ljmin:ljmax),fileresu(1))
!#endif

! Create sedimentological output file
!cval il y va y avoir sans doute un gros pb de creation de fichier, il faut
! creer fichier si nbouc==0 et l_initfromfile
! puis eliminer les 4 lignes ci-dessous
!#if defined 1 && defined 1
!     IF (nb==1 .AND. .NOT.l_filebydate(nb)) CALL sed_output(nb,t,ssh(limin:limax,ljmin:ljmax),fileresu(1))
!#endif
!#ifdef key_biolo
!       IF (l_out_subs_diag .AND. nb==1 .AND. .NOT.l_filebydate(nb) ) CALL output_diag(t,ssh(limin:limax,ljmin:ljmax),nb,fileresu(1))
!#endif

! Allocate and initialize arrays for cumulative atmospherical fluxes
     IF (nbave /= 0 .AND. l_out_timeave(nb)) THEN
       IF (l_out_ssh(nb)) THEN
         ALLOCATE(ssh_ave(nbave,limin:limax,ljmin:ljmax) )
         ssh_ave(:,:,:)     = 0.0_rsh
       END IF
       IF (l_out_ubt(nb)) THEN
         ALLOCATE(  u_ave(nbave,limin:limax,ljmin:ljmax) )
         ALLOCATE(  v_ave(nbave,limin:limax,ljmin:ljmax) )
         u_ave(:,:,:)     = 0.0_rsh
         v_ave(:,:,:)     = 0.0_rsh
       END IF
       IF (l_out_u3d(nb).AND..NOT.l_modele2d) THEN
         ALLOCATE( uz_ave(nbave,kmax,limin:limax,ljmin:ljmax) )
         ALLOCATE( vz_ave(nbave,kmax,limin:limax,ljmin:ljmax) )
         uz_ave(:,:,:,:)     = 0.0_rsh
         vz_ave(:,:,:,:)     = 0.0_rsh

       END IF
       IF (l_out_temp(nb)) THEN
         ALLOCATE( temp_ave(nbave,kmax,limin:limax,ljmin:ljmax) )
         temp_ave(:,:,:,:)     = 0.0_rsh
       END IF
       IF (l_out_sal(nb)) THEN
         ALLOCATE( sal_ave(nbave,kmax,limin:limax,ljmin:ljmax) )
         sal_ave(:,:,:,:)     = 0.0_rsh
       END IF
       IF (l_out_visc(nb).AND.l_modele2d) THEN
         ALLOCATE(vis2d_ave(nbave,limin:limax,ljmin:ljmax) )
         vis2d_ave(:,:,:)     = 0.0_rsh
       ELSE IF (l_out_visc(nb).AND..NOT.l_modele2d) THEN
         ALLOCATE(vish_xe_ave(nbave,kmax,limin:limax,ljmin:ljmax) )
         ALLOCATE(vish_phi_ave(nbave,kmax,limin:limax,ljmin:ljmax) )
         vish_phi_ave(:,:,:,:)     = 0.0_rsh
         vish_xe_ave(:,:,:,:)     = 0.0_rsh
       END IF
       IF (l_out_turb(nb).AND..NOT.l_modele2d) THEN
         ALLOCATE(nz_ave(nbave,0:kmax,limin:limax,ljmin:ljmax) )
         ALLOCATE(kz_ave(nbave,0:kmax,limin:limax,ljmin:ljmax) )
         nz_ave(:,:,:,:)     = 0.0_rsh
         kz_ave(:,:,:,:)     = 0.0_rsh
       END IF
       IF (l_out_buoy(nb).AND..NOT.l_modele2d) THEN
         ALLOCATE(boz_ave(nbave,kmax,limin:limax,ljmin:ljmax) )
         boz_ave(:,:,:,:)     = 0.0_rsh
       END IF
       ALLOCATE(dt_ave(nbave) )
       dt_ave(:) = 0.0_rlg
     END IF



    END DO    ! End loop on number of output file


! Allocate and initialize arrays for cumulative atmospherical fluxes
     IF (nbatm /= 0) THEN
       ALLOCATE(rad_cum(nbatm,limin:limax,ljmin:ljmax) )
       ALLOCATE(irlos_cum(nbatm,limin:limax,ljmin:ljmax) )
       ALLOCATE(latlos_cum(nbatm,limin:limax,ljmin:ljmax) )
       ALLOCATE(senslos_cum(nbatm,limin:limax,ljmin:ljmax) )
       ALLOCATE(evap_cum(nbatm,limin:limax,ljmin:ljmax) )
       ALLOCATE(rain_cum(nbatm,limin:limax,ljmin:ljmax) )
       ALLOCATE(tosx_cum(nbatm,limin:limax,ljmin:ljmax) )
       ALLOCATE(tosy_cum(nbatm,limin:limax,ljmin:ljmax) )
       ALLOCATE(dt_cum(nbatm),dt_ave(nbatm) )
       rad_cum(:,:,:)     = 0.0_rsh
       irlos_cum(:,:,:)   = 0.0_rsh
       latlos_cum(:,:,:)  = 0.0_rsh
       senslos_cum(:,:,:) = 0.0_rsh
       evap_cum(:,:,:)    = 0.0_rsh
       rain_cum(:,:,:)    = 0.0_rsh
       tosx_cum(:,:,:)    = 0.0_rsh
       tosy_cum(:,:,:)    = 0.0_rsh
       dt_cum(:) = 0.0_rlg

     END IF

   END IF     ! Ends test on nbouc

! AFTER INITIALIZATION IS DONE
! ----------------------------

! GENERAL CASE
! ------------

! Loop on number of output file
   DO nb = 1,nb_out



! averaged nudging fields
! -----------------------




! estimate averaged dynamical and thermodynamical variables
! ---------------------------------------------------------

     IF (nbave /= 0 .AND. l_out_timeave(nb) ) THEN

       IF ( t >= tdebsor(nb) .AND. t < tsort(nb) .AND. t <= tfinsor(nb) ) THEN
! add variables
         IF (l_out_ssh(nb)) &
           ssh_ave(nb,:,:) = ssh_ave(nb,:,:) + halfdt * ssh(limin:limax,ljmin:ljmax)
         IF (l_out_ubt(nb)) THEN
           u_ave(nb,:,:) = u_ave(nb,:,:) + halfdt * u(limin:limax,ljmin:ljmax)
           v_ave(nb,:,:) = v_ave(nb,:,:) + halfdt * v(limin:limax,ljmin:ljmax)
         END IF
         IF (l_out_u3d(nb).AND..NOT.l_modele2d) THEN
           uz_ave(nb,:,:,:) = uz_ave(nb,:,:,:) + halfdt * uz(:,limin:limax,ljmin:ljmax)
           vz_ave(nb,:,:,:) = vz_ave(nb,:,:,:) + halfdt * vz(:,limin:limax,ljmin:ljmax)

         END IF
         IF (l_out_temp(nb)) &
           temp_ave(nb,:,:,:) = temp_ave(nb,:,:,:) + halfdt * temp(:,limin:limax,ljmin:ljmax)
         IF (l_out_sal(nb)) &
           sal_ave(nb,:,:,:) = sal_ave(nb,:,:,:) + halfdt * sal(:,limin:limax,ljmin:ljmax)
         IF (l_out_visc(nb).AND.l_modele2d) THEN
           vis2d_ave(nb,:,:) = vis2d_ave(nb,:,:) + halfdt * vis2d(limin:limax,ljmin:ljmax)
         ELSE IF (l_out_visc(nb).AND..NOT.l_modele2d) THEN
           vish_xe_ave(nb,:,:,:) = vish_xe_ave(nb,:,:,:) + halfdt * vish_xe(:,limin:limax,ljmin:ljmax)
           vish_phi_ave(nb,:,:,:) = vish_phi_ave(nb,:,:,:) + halfdt * vish_phi(:,limin:limax,ljmin:ljmax)
         END IF
         IF (l_out_turb(nb).AND..NOT.l_modele2d) THEN
           nz_ave(nb,:,:,:) = nz_ave(nb,:,:,:) + halfdt * nz(:,limin:limax,ljmin:ljmax)
           kz_ave(nb,:,:,:) = kz_ave(nb,:,:,:) + halfdt * kz(:,limin:limax,ljmin:ljmax)
         END IF
         IF (l_out_buoy(nb).AND..NOT.l_modele2d) &
           boz_ave(nb,:,:,:) = boz_ave(nb,:,:,:) + halfdt * boz(:,limin:limax,ljmin:ljmax)
         dt_ave(nb) = dt_ave(nb) + halfdt
       END IF

       IF ( t >= tsort(nb) .AND. t <= tfinsor(nb) ) THEN
! end adding of fluxes
         dtp = tsort(nb) - (t-halfdt)
         dtn = t - tsort(nb)
         IF (l_out_ssh(nb)) &
         ssh_ave(nb,:,:) = ssh_ave(nb,:,:) + ssh(limin:limax,ljmin:ljmax)*dtp
         IF (l_out_ubt(nb)) THEN
           u_ave(nb,:,:) = u_ave(nb,:,:) + u(limin:limax,ljmin:ljmax)*dtp
           v_ave(nb,:,:) = v_ave(nb,:,:) + v(limin:limax,ljmin:ljmax)*dtp
         END IF
         IF (l_out_u3d(nb).AND..NOT.l_modele2d) THEN
           uz_ave(nb,:,:,:) = uz_ave(nb,:,:,:) + dtp * uz(:,limin:limax,ljmin:ljmax)
           vz_ave(nb,:,:,:) = vz_ave(nb,:,:,:) + dtp * vz(:,limin:limax,ljmin:ljmax)

         END IF
         IF (l_out_temp(nb)) &
           temp_ave(nb,:,:,:) = temp_ave(nb,:,:,:) + dtp * temp(:,limin:limax,ljmin:ljmax)
         IF (l_out_sal(nb)) &
           sal_ave(nb,:,:,:) = sal_ave(nb,:,:,:) + dtp * sal(:,limin:limax,ljmin:ljmax)
         IF (l_out_visc(nb).AND.l_modele2d) THEN
           vis2d_ave(nb,:,:) = vis2d_ave(nb,:,:) + dtp * vis2d(limin:limax,ljmin:ljmax)
         ELSE IF (l_out_visc(nb).AND..NOT.l_modele2d) THEN
           vish_xe_ave(nb,:,:,:) = vish_xe_ave(nb,:,:,:) + dtp * vish_xe(:,limin:limax,ljmin:ljmax)
           vish_phi_ave(nb,:,:,:) = vish_phi_ave(nb,:,:,:) + dtp * vish_phi(:,limin:limax,ljmin:ljmax)
         END IF
         IF (l_out_turb(nb).AND..NOT.l_modele2d) THEN
           nz_ave(nb,:,:,:) = nz_ave(nb,:,:,:) + dtp * nz(:,limin:limax,ljmin:ljmax)
           kz_ave(nb,:,:,:) = kz_ave(nb,:,:,:) + dtp * kz(:,limin:limax,ljmin:ljmax)
         END IF
         IF (l_out_buoy(nb).AND..NOT.l_modele2d) &
           boz_ave(nb,:,:,:) = boz_ave(nb,:,:,:) + dtp * boz(:,limin:limax,ljmin:ljmax)
         dt_ave(nb) =  dt_ave(nb) + dtp

! estimate averaged variables
         IF (dt_ave(nb) /= 0.0_rlg) THEN
           ssh_ave(nb,:,:) = ssh_ave(nb,:,:) / dt_ave(nb)
           IF (l_out_ubt(nb)) THEN
             u_ave(nb,:,:) = u_ave(nb,:,:) / dt_ave(nb)
             v_ave(nb,:,:) = v_ave(nb,:,:) / dt_ave(nb)
           END IF
           IF (l_out_u3d(nb).AND..NOT.l_modele2d) THEN
             uz_ave(nb,:,:,:) = uz_ave(nb,:,:,:) / dt_ave(nb)
             vz_ave(nb,:,:,:) = vz_ave(nb,:,:,:) / dt_ave(nb)

           END IF
           IF (l_out_temp(nb)) &
             temp_ave(nb,:,:,:) = temp_ave(nb,:,:,:) / dt_ave(nb)
           IF (l_out_sal(nb)) &
             sal_ave(nb,:,:,:) = sal_ave(nb,:,:,:) / dt_ave(nb)
           IF (l_out_visc(nb).AND.l_modele2d) THEN
             vis2d_ave(nb,:,:) = vis2d_ave(nb,:,:) / dt_ave(nb)
           ELSE IF (l_out_visc(nb).AND..NOT.l_modele2d) THEN
             vish_xe_ave(nb,:,:,:) = vish_xe_ave(nb,:,:,:) / dt_ave(nb)
             vish_phi_ave(nb,:,:,:) = vish_phi_ave(nb,:,:,:) / dt_ave(nb)
           END IF
           IF (l_out_turb(nb).AND..NOT.l_modele2d) THEN
             nz_ave(nb,:,:,:) = nz_ave(nb,:,:,:) / dt_ave(nb)
             kz_ave(nb,:,:,:) = kz_ave(nb,:,:,:) / dt_ave(nb)
           END IF
           IF (l_out_buoy(nb).AND..NOT.l_modele2d) &
             boz_ave(nb,:,:,:) = boz_ave(nb,:,:,:) / dt_ave(nb)
         END IF
       END IF

     END IF


! estimate cumulative atmospherical surface fluxes
! ------------------------------------------------

     IF (nbatm /= 0 .AND. l_out_atmsflx(nb) ) THEN

       IF ( t/=tdeb .AND. t >= tdebsor(nb)-pasor(nb) .AND. t < tsort(nb) .AND. t <= tfinsor(nb) ) THEN
! add fluxes
         rad_cum(nb,:,:)     = rad_cum(nb,:,:)     + halfdt * rad(limin:limax,ljmin:ljmax)
         irlos_cum(nb,:,:)   = irlos_cum(nb,:,:)   + halfdt * irlos(limin:limax,ljmin:ljmax)
         latlos_cum(nb,:,:)  = latlos_cum(nb,:,:)  + halfdt * latlos(limin:limax,ljmin:ljmax)
         senslos_cum(nb,:,:) = senslos_cum(nb,:,:) + halfdt * senslos(limin:limax,ljmin:ljmax)
         evap_cum(nb,:,:)    = evap_cum(nb,:,:)    + halfdt * evap(limin:limax,ljmin:ljmax)
         rain_cum(nb,:,:)    = rain_cum(nb,:,:)    + halfdt * rain(limin:limax,ljmin:ljmax)
         tosx_cum(nb,:,:)    = tosx_cum(nb,:,:)    + halfdt * tosx(limin:limax,ljmin:ljmax)
         tosy_cum(nb,:,:)    = tosy_cum(nb,:,:)    + halfdt * tosy(limin:limax,ljmin:ljmax)

         dt_cum(nb) =  dt_cum(nb) + halfdt
       END IF

       IF ( t >= tsort(nb) .AND. t <= tfinsor(nb) ) THEN
! end adding of fluxes
         dtp = tsort(nb) - (t-halfdt)
         dtn = t - tsort(nb)
         rad_cum(nb,:,:)     = rad_cum(nb,:,:)     + rad(limin:limax,ljmin:ljmax)*dtp
         irlos_cum(nb,:,:)   = irlos_cum(nb,:,:)   + irlos(limin:limax,ljmin:ljmax)*dtp
         latlos_cum(nb,:,:)  = latlos_cum(nb,:,:)  + latlos(limin:limax,ljmin:ljmax)*dtp
         senslos_cum(nb,:,:) = senslos_cum(nb,:,:) + senslos(limin:limax,ljmin:ljmax)*dtp
         evap_cum(nb,:,:)    = evap_cum(nb,:,:)    + evap(limin:limax,ljmin:ljmax)*dtp
         rain_cum(nb,:,:)    = rain_cum(nb,:,:)    + rain(limin:limax,ljmin:ljmax)*dtp
         tosx_cum(nb,:,:)    = tosx_cum(nb,:,:)    + tosx(limin:limax,ljmin:ljmax)*dtp
         tosy_cum(nb,:,:)    = tosy_cum(nb,:,:)    + tosy(limin:limax,ljmin:ljmax)*dtp

         dt_cum(nb) =  dt_cum(nb) + dtp
! estimate cumulative fluxes
         IF (dt_cum(nb) /= 0.0_rlg) THEN
           rad_cum(nb,:,:)     = rad_cum(nb,:,:)     / dt_cum(nb)    ! units
           irlos_cum(nb,:,:)   = irlos_cum(nb,:,:)   / dt_cum(nb)
           latlos_cum(nb,:,:)  = latlos_cum(nb,:,:)  / dt_cum(nb)
           senslos_cum(nb,:,:) = senslos_cum(nb,:,:) / dt_cum(nb)
           evap_cum(nb,:,:)    = evap_cum(nb,:,:) * 0.001_rsh        ! evap(:,:) and rain(:,:) are in mm/s
           rain_cum(nb,:,:)    = rain_cum(nb,:,:) * 0.001_rsh        ! expected output is in m
           tosx_cum(nb,:,:)    = tosx_cum(nb,:,:)    / dt_cum(nb)
           tosy_cum(nb,:,:)    = tosy_cum(nb,:,:)    / dt_cum(nb)

         END IF
       END IF

     END IF


! save fields without time interpolation or filtering
! ---------------------------------------------------

   IF ( (.NOT.l_out_precdate(nb) .AND. .NOT.l_out_demer(nb)) .OR. nbouc==0 ) THEN   ! cas des sorties brutes
     IF ( (t >= tsort(nb).AND.t <= tfinsor(nb)) .OR. nbouc==0 ) THEN

       IF ( l_out_off(nb) ) THEN



       ELSE


         CALL output_save &
                (nb,t,u(limin:limax,ljmin:ljmax),      uz(:,limin:limax,ljmin:ljmax), &
                    ssh(limin:limax,ljmin:ljmax),                                     &
                      v(limin:limax,ljmin:ljmax),      vz(:,limin:limax,ljmin:ljmax), &
                  sal(:,limin:limax,ljmin:ljmax),    temp(:,limin:limax,ljmin:ljmax), &
                  vis2d(limin:limax,ljmin:ljmax),                                     &
              vish_xe(:,limin:limax,ljmin:ljmax),vish_phi(:,limin:limax,ljmin:ljmax), &
                   nz(:,limin:limax,ljmin:ljmax),      kz(:,limin:limax,ljmin:ljmax), &
                   boz(:,limin:limax,ljmin:ljmax)                                     &


            ,cv_wat(:,:,limin:limax,ljmin:ljmax)                                      &

                )
       END IF






       IF (nb==1) CALL sed_output(nb,t,ssh(limin:limax,ljmin:ljmax),fileresu(1))



!       IF (t >= tsort(nb).AND.t <= tfinsor(nb)) tsort(nb)=tsort(nb)+pasor(nb)
!to allow correct cumulative flux even with a very small time step
       IF (t >= tsort(nb).AND.t <= tfinsor(nb)) tsort(nb) = MAX(t+halfdt,t+pasor(nb))
!to get the same time steps in results.nc and diags.nc         tsort(nb) = tsort(nb)+pasor(nb)



     END IF
   END IF      ! fin du cas des sorties brutes

! SAVE OUTPUTS AT A PRECISE TIME
! ------------------------------

   IF (l_out_precdate(nb)) THEN                   ! cas des sorties a instant precis

     IF ((t < tinter(nb) .AND. t+0.5_rlg*dt >=tinter(nb) ) .OR. &
         (t ==tinter(nb) .AND. l_initfromfile)) THEN
       vis2da(limin:limax,ljmin:ljmax)=vis2d(limin:limax,ljmin:ljmax)
       sala(:,limin:limax,ljmin:ljmax)=sal(:,limin:limax,ljmin:ljmax)
       tempa(:,limin:limax,ljmin:ljmax)=temp(:,limin:limax,ljmin:ljmax)
       vish_xea(:,limin:limax,ljmin:ljmax)=vish_xe(:,limin:limax,ljmin:ljmax)
       vish_phia(:,limin:limax,ljmin:ljmax)=vish_phi(:,limin:limax,ljmin:ljmax)
       nza(:,limin:limax,ljmin:ljmax)=nz(:,limin:limax,ljmin:ljmax)
       kza(:,limin:limax,ljmin:ljmax)=kz(:,limin:limax,ljmin:ljmax)
       boza(:,limin:limax,ljmin:ljmax)=boz(:,limin:limax,ljmin:ljmax)
       xea(limin:limax,ljmin:ljmax)=ssh(limin:limax,ljmin:ljmax)
       ua(limin:limax,ljmin:ljmax)=u(limin:limax,ljmin:ljmax)
       IF (.NOT.l_modele2d) uza(:,limin:limax,ljmin:ljmax)=uz(:,limin:limax,ljmin:ljmax)
       va(limin:limax,ljmin:ljmax)=v(limin:limax,ljmin:ljmax)
       IF (.NOT.l_modele2d) vza(:,limin:limax,ljmin:ljmax)=vz(:,limin:limax,ljmin:ljmax)

       cv_wata(:,:,limin:limax,ljmin:ljmax)=cv_wat(:,:,limin:limax,ljmin:ljmax)

       ta=t

     ELSEIF (t >= tinter(nb)) THEN
       CALL tool_timeinterp_2dvar(tinter(nb),ta,t,vis2da(limin:limax,ljmin:ljmax),vis2d(limin:limax,ljmin:ljmax))
       CALL tool_timeinterp_3dvar(tinter(nb),ta,t,sala(:,limin:limax,ljmin:ljmax),sal(:,limin:limax,ljmin:ljmax))
       CALL tool_timeinterp_3dvar(tinter(nb),ta,t,tempa(:,limin:limax,ljmin:ljmax),temp(:,limin:limax,ljmin:ljmax))
       CALL tool_timeinterp_3dvar(tinter(nb),ta,t,vish_xea(:,limin:limax,ljmin:ljmax),vish_xe(:,limin:limax,ljmin:ljmax))
       CALL tool_timeinterp_3dvar(tinter(nb),ta,t,vish_phia(:,limin:limax,ljmin:ljmax),vish_phi(:,limin:limax,ljmin:ljmax))
       CALL tool_timeinterp_3dvar(tinter(nb),ta,t,nza(1:kmax,limin:limax,ljmin:ljmax),nz(1:kmax,limin:limax,ljmin:ljmax))
       CALL tool_timeinterp_3dvar(tinter(nb),ta,t,kza(1:kmax,limin:limax,ljmin:ljmax),kz(1:kmax,limin:limax,ljmin:ljmax))
       CALL tool_timeinterp_3dvar(tinter(nb),ta,t,boza(:,limin:limax,ljmin:ljmax),boz(:,limin:limax,ljmin:ljmax))

       DO l=1,nv_state
         c1(:,limin:limax,ljmin:ljmax)=cv_wata(l,:,limin:limax,ljmin:ljmax)
         c2(:,limin:limax,ljmin:ljmax)=cv_wat(l,:,limin:limax,ljmin:ljmax)
         CALL tool_timeinterp_3dvar(tinter(nb),ta,t,c1(:,limin:limax,ljmin:ljmax),c2(:,limin:limax,ljmin:ljmax))
         cv_wata(l,:,limin:limax,ljmin:ljmax)=c1(:,limin:limax,ljmin:ljmax)
       END DO

       CALL tool_timeinterp_2dvar(tinter(nb),ta,t,xea(limin:limax,ljmin:ljmax),ssh(limin:limax,ljmin:ljmax))
       CALL tool_timeinterp_2dvar(tinter(nb),ta,t,ua(limin:limax,ljmin:ljmax),u(limin:limax,ljmin:ljmax))
       IF (.NOT.l_modele2d) &
           CALL tool_timeinterp_3dvar(tinter(nb),ta,t,uza(:,limin:limax,ljmin:ljmax),uz(:,limin:limax,ljmin:ljmax))
       CALL tool_timeinterp_2dvar(tinter(nb),ta,t,va(limin:limax,ljmin:ljmax),v(limin:limax,ljmin:ljmax))
       IF (.NOT.l_modele2d) &
           CALL tool_timeinterp_3dvar(tinter(nb),ta,t,vza(:,limin:limax,ljmin:ljmax),vz(:,limin:limax,ljmin:ljmax))
       CALL output_save &
            (nb,tinter(nb),ua(limin:limax,ljmin:ljmax),  uza(:,limin:limax,ljmin:ljmax),     &
                          xea(limin:limax,ljmin:ljmax),                                      &
                           va(limin:limax,ljmin:ljmax),  vza(:,limin:limax,ljmin:ljmax),     &
                       sala(:,limin:limax,ljmin:ljmax),tempa(:,limin:limax,ljmin:ljmax),     &
                       vis2da(limin:limax,ljmin:ljmax),                                      &
                   vish_xea(:,limin:limax,ljmin:ljmax),vish_phia(:,limin:limax,ljmin:ljmax), &
                        nza(:,limin:limax,ljmin:ljmax),      kza(:,limin:limax,ljmin:ljmax), &
                        boza(:,limin:limax,ljmin:ljmax)                                       &

                 ,cv_wata(:,:,limin:limax,ljmin:ljmax)                                       &

            )
! taking account output for diag and sedim if precdate = true


       IF (nb==1) CALL sed_output(nb,t,ssh(limin:limax,ljmin:ljmax),fileresu(1))


       tinter(nb)=tinter(nb)+pasor(nb)
     END IF
   END IF                      ! fin du traitement des sorties precises

! SAVE FILTERED OUTPUTS
! ---------------------

   IF (l_out_demer(nb)) THEN             ! traitement des sorties moyennees

     IF (t < tinter(nb) .AND. t+0.5*dt >=tinter(nb) ) THEN
       sala(:,limin:limax,ljmin:ljmax)=sal(:,limin:limax,ljmin:ljmax)
       tempa(:,limin:limax,ljmin:ljmax)=temp(:,limin:limax,ljmin:ljmax)
       xea(limin:limax,ljmin:ljmax)=ssh(limin:limax,ljmin:ljmax)
       ua(limin:limax,ljmin:ljmax)=u(limin:limax,ljmin:ljmax)
       IF (.NOT.l_modele2d) uza(:,limin:limax,ljmin:ljmax)=uz(:,limin:limax,ljmin:ljmax)
       va(limin:limax,ljmin:ljmax)=v(limin:limax,ljmin:ljmax)
       IF (.NOT.l_modele2d) vza(:,limin:limax,ljmin:ljmax)=vz(:,limin:limax,ljmin:ljmax)

     cv_wata(:,:,limin:limax,ljmin:ljmax)=cv_wat(:,:,limin:limax,ljmin:ljmax)

     ta=t

   ELSE IF (t >= tinter(nb)) THEN
     CALL tool_timeinterp_3dvar(tinter(nb),ta,t,sala(:,limin:limax,ljmin:ljmax),sal(:,limin:limax,ljmin:ljmax))
     CALL tool_timeinterp_3dvar(tinter(nb),ta,t,tempa(:,limin:limax,ljmin:ljmax),temp(:,limin:limax,ljmin:ljmax))

     DO l=1,nv_state
       c1(:,limin:limax,ljmin:ljmax)=cv_wata(l,:,limin:limax,ljmin:ljmax)
       c2(:,limin:limax,ljmin:ljmax)=cv_wat(l,:,limin:limax,ljmin:ljmax)
       CALL tool_timeinterp_3dvar(tinter(nb),ta,t,c1(:,limin:limax,ljmin:ljmax),c2(:,limin:limax,ljmin:ljmax))
       cv_wata(l,:,limin:limax,ljmin:ljmax)=c1(:,limin:limax,ljmin:ljmax)
     END DO

     CALL tool_timeinterp_2dvar(tinter(nb),ta,t,xea(limin:limax,ljmin:ljmax),ssh(limin:limax,ljmin:ljmax))
     CALL tool_timeinterp_2dvar(tinter(nb),ta,t,ua(limin:limax,ljmin:ljmax),u(limin:limax,ljmin:ljmax))
     IF (.NOT.l_modele2d) &
          CALL tool_timeinterp_3dvar(tinter(nb),ta,t,uza(:,limin:limax,ljmin:ljmax),uz(:,limin:limax,ljmin:ljmax))
     CALL tool_timeinterp_2dvar(tinter(nb),ta,t,va(limin:limax,ljmin:ljmax),v(limin:limax,ljmin:ljmax))
     IF (.NOT.l_modele2d) &
          CALL tool_timeinterp_3dvar(tinter(nb),ta,t,vza(:,limin:limax,ljmin:ljmax),vz(:,limin:limax,ljmin:ljmax))
     CALL out_demerliac(xea(  limin:limax,ljmin:ljmax),                                  &
                         ua(  limin:limax,ljmin:ljmax),  uza(:,limin:limax,ljmin:ljmax), &
                         va(  limin:limax,ljmin:ljmax),  vza(:,limin:limax,ljmin:ljmax), &
                       sala(:,limin:limax,ljmin:ljmax),tempa(:,limin:limax,ljmin:ljmax), &
                        xed(  limin:limax,ljmin:ljmax),                                  &
                         ud(  limin:limax,ljmin:ljmax),  uzd(:,limin:limax,ljmin:ljmax), &
                         vd(  limin:limax,ljmin:ljmax),  vzd(:,limin:limax,ljmin:ljmax), &
                       sald(:,limin:limax,ljmin:ljmax),tempd(:,limin:limax,ljmin:ljmax), &

                    cv_wata(:,:,limin:limax,ljmin:ljmax),                                &
                    cv_watd(:,:,limin:limax,ljmin:ljmax),                                &

                        nsdem)

       IF (nsdem == 71) THEN
         tect=tinter(nb)-35.0_rlg*3600.0_rlg

         CALL output_save &
              (nb,tect,ud(limin:limax,ljmin:ljmax),  uzd(:,limin:limax,ljmin:ljmax),    &
                      xed(limin:limax,ljmin:ljmax),                                     &
                       vd(limin:limax,ljmin:ljmax),  vzd(:,limin:limax,ljmin:ljmax),    &
                   sald(:,limin:limax,ljmin:ljmax),tempd(:,limin:limax,ljmin:ljmax),    &
                    vis2d(limin:limax,ljmin:ljmax),                                     &
                vish_xe(:,limin:limax,ljmin:ljmax),vish_phi(:,limin:limax,ljmin:ljmax), &
                     nz(:,limin:limax,ljmin:ljmax),      kz(:,limin:limax,ljmin:ljmax), &
                     boz(:,limin:limax,ljmin:ljmax)                                     &

             ,cv_watd(:,:,limin:limax,ljmin:ljmax)                                      &

              )
         tinter(nb)=tect+pasor(nb)-35.0_rlg*3600.0_rlg
         nsdem=0
       ELSE
         tinter(nb)=tinter(nb)+3600.0_rlg
       END IF
     END IF
   END IF ! fin du traitement des sorties moyennees

   END DO ! End loop on number of output files

   !PRINT*, 'exit OUTPUT_MNG routine'

  END SUBROUTINE output_mng

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


 END MODULE output
