  SUBROUTINE out_demerliac(xe,u,uz,v,vz,sal,temp,                       &
                           xei,ui,uzi,vi,vzi,sali,tempi,                &

                                                        cv_wat,cv_wati, &

                                                        nsdem)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE out_demerliac  ***
!&E
!&E ** Purpose : Demerliac filtering for all quantities
!&E
!&E ** Description :
!&E
!&E ** Called by : gessor
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : xei,sali,tempi,cv_wati,ui,vi,uzi,vzi,xe,u,v,uz,vz,sal,temp,cv_wat
!&E
!&E ** Modified variables : xei,sali,tempi,cv_wati,ui,vi,uzi,vzi
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2004-08-17
!&E       !  2009-04    (A. Mano) MPI2D (not tested yet)
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d,  ONLY : nbouc, ig,id,l_modele2d

   USE comsubstance, ONLY : nv_state

   IMPLICIT NONE

!! * Arguments
   REAL(kind=rsh),DIMENSION(     limin:limax,ljmin:ljmax),INTENT(in)    :: xe,u,v
   REAL(kind=rsh),DIMENSION(     limin:limax,ljmin:ljmax),INTENT(inout) :: xei,ui,vi
   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(inout) :: uzi,vzi,sali,tempi

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

   INTEGER,INTENT(inout) :: nsdem

!! * Local declarations
   INTEGER :: i,j,k,l
   REAL(kind=rsh),DIMENSION(71) :: cofi

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

   !PRINT*, 'enter OUTDEMERLIAC  routine'

   cofi(36:71)=(/768,766,762,752,738,726,704,678,658,624,586,558,512,465,435,392,    &
                 351,325,288,253,231,200,171,153,128,105,91,72,55,45,32,21,15,8,3,1/)
   cofi(1:35)=cofi(71:37:-1)

   IF(nbouc==0) THEN
     xei(:,:)=-valmanq
     sali(:,:,:)=valmanq
     tempi(:,:,:)=valmanq

     cv_wati(:,:,:,:)=valmanq

   ELSE 
     IF(nsdem == 0) THEN
       xei(:,:)    =0.0_rsh
       ui(:,:)     =0.0_rsh
       vi(:,:)     =0.0_rsh
       uzi(:,:,:)  =0.0_rsh
       vzi(:,:,:)  =0.0_rsh
       sali(:,:,:) =0.0_rsh
       tempi(:,:,:)=0.0_rsh

       cv_wati(:,:,:,:)=0.0_rsh

     ENDIF

     nsdem=nsdem+1

     xei(limin:limax,ljmin:ljmax)    =xei(limin:limax,ljmin:ljmax)    +cofi(nsdem)*xe(limin:limax,ljmin:ljmax)
     ui(limin:limax,ljmin:ljmax)     =ui(limin:limax,ljmin:ljmax)     +cofi(nsdem)*u(limin:limax,ljmin:ljmax)
     vi(limin:limax,ljmin:ljmax)     =vi(limin:limax,ljmin:ljmax)     +cofi(nsdem)*v(limin:limax,ljmin:ljmax)
     IF (.NOT.l_modele2d) THEN
       uzi(:,limin:limax,ljmin:ljmax)  =uzi(:,limin:limax,ljmin:ljmax)  +cofi(nsdem)*uz(:,limin:limax,ljmin:ljmax)
       vzi(:,limin:limax,ljmin:ljmax)  =vzi(:,limin:limax,ljmin:ljmax)  +cofi(nsdem)*vz(:,limin:limax,ljmin:ljmax)
     END IF
     sali(:,limin:limax,ljmin:ljmax) =sali(:,limin:limax,ljmin:ljmax) +cofi(nsdem)*sal(:,limin:limax,ljmin:ljmax)
     tempi(:,limin:limax,ljmin:ljmax)=tempi(:,limin:limax,ljmin:ljmax)+cofi(nsdem)*temp(:,limin:limax,ljmin:ljmax)

     DO l=1,nv_state
       cv_wati(l,:,limin:limax,ljmin:ljmax)=cv_wati(l,:,limin:limax,ljmin:ljmax)+cofi(nsdem)*cv_wat(l,:,limin:limax,ljmin:ljmax)
     ENDDO


     IF(nsdem == 71) THEN
       DO j=ljmin,ljmax
!DO i=MAX0(limin,ig(j)),MIN0(limax,id(j))
         DO i=limin,limax
           xei(i,j)=xei(i,j)/SUM(cofi(:))
           IF(ABS(xei(i,j)-valmanq)/valmanq < epsilon) xei(i,j)=valmanq
           ui(i,j)=ui(i,j)/SUM(cofi(:))
           IF(ABS(ui(i,j)) > 20.0_rsh) ui(i,j)=0.0_rsh
           vi(i,j)=vi(i,j)/SUM(cofi(:))
           IF(ABS(vi(i,j)) > 20.0_rsh) vi(i,j)=0.0_rsh
           DO k=1,kmax
             sali(k,i,j)=sali(k,i,j)/SUM(cofi(:))
             IF(ABS(sali(k,i,j)-valmanq)/valmanq < epsilon) sali(k,i,j)=valmanq
             tempi(k,i,j)=tempi(k,i,j)/SUM(cofi(:))
             IF(ABS(tempi(k,i,j)-valmanq)/valmanq < epsilon) tempi(k,i,j)=valmanq
             uzi(k,i,j)=uzi(k,i,j)/SUM(cofi(:))
             IF(ABS(uzi(k,i,j)) > 20.0_rsh) uzi(k,i,j)=0.0_rsh
             vzi(k,i,j)=vzi(k,i,j)/SUM(cofi(:))
             IF(ABS(vzi(k,i,j)) > 20.0_rsh) vzi(k,i,j)=0.0_rsh

             DO l=1,nv_state
               cv_wati(l,k,i,j)=cv_wati(l,k,i,j)/SUM(cofi(:))
               IF(ABS(cv_wati(l,k,i,j)-valmanq)/valmanq < epsilon) cv_wati(l,k,i,j)=valmanq
             ENDDO

           ENDDO
         ENDDO
       ENDDO
     ENDIF
   ENDIF

   !PRINT*, 'exit OUTDEMERLIAC  routine'

  END SUBROUTINE out_demerliac
