MODULE waveforcing

  !!======================================================================
  !!                   ***  MODULE waveforcing  ***
  !! Compute the surface forcing associated WITH waves.
  !!======================================================================
  !! * Modules used
  USE parameters

  IMPLICIT NONE
  PRIVATE

  !! * Accessibility
  PUBLIC init_wave_forcing_fric,lect_wave,compute_wave_forcing, div_tens_rad ! routine called by main.F90,
  ! dyn2dxy.F90, dyn2dyx.F90
  PUBLIC lec_para_swan

  !! * Shared or PUBLIC MODULE variables
  LOGICAL,PUBLIC :: l_wave_brgm_phillips   ! Phillips theory to account for wave forcing  on circulation
  LOGICAL,PUBLIC :: l_wave_brgm_smith      ! Smith theory to account for wave forcing on circulation
  LOGICAL,PUBLIC :: l_wave_brgm_roller     ! account for roller forcing on circulation
  LOGICAL,PUBLIC :: l_wave_brgm_t2c        ! something to DO WITH computation of radiation stress !!!
  LOGICAL,PUBLIC :: l_wave_brgm_asym       ! account for wave asymetry
  LOGICAL,PUBLIC :: l_wave_brgm_frot_homog ! homogeneous bottom drag
  LOGICAL,PUBLIC :: l_stat_conv       ! homogeneous bottom drag
  INTEGER,PUBLIC :: lect_sens         ! readding mode for SWAN bathymetry
  INTEGER,PUBLIC :: iter_max          ! maximum number of iteration
  REAL(kind=rsh),PUBLIC :: fact_bathy ! bathymetry correction factor
  REAL(kind=rsh),PUBLIC :: M_battjes  ! coefficient for battjes formulation of
  ! turbulent horizontal mixing
  REAL(kind=rsh),PUBLIC :: gammab     ! wave dissipation
  REAL(kind=rsh),PUBLIC :: hlim_calc  ! minimum depth for computation
  REAL(kind=rsh),PUBLIC :: frot_val   ! drag friction coefficient
  REAL(kind=rsh),PUBLIC :: level      ! level shift
  REAL(kind=rsh),PUBLIC :: err_conv   ! level shift

  INTEGER       , PUBLIC,DIMENSION(imin:imax,jmin:jmax) :: front
  REAL(kind=rsh), PUBLIC,DIMENSION(imin:imax,jmin:jmax) :: Sxx,Syy,Sxy,hs,qb
  REAL(kind=rsh), PUBLIC,DIMENSION(imin:imax,jmin:jmax) :: Qx,Qy,Orb,Fradx,Frady
  REAL(kind=rsh), PUBLIC,DIMENSION(imin:imax,jmin:jmax) :: Depth, Tm01, Dissip,tp
  REAL(kind=rsh), PUBLIC,DIMENSION(imin:imax,jmin:jmax) :: Direc,Wlen,Dissip_u,Dissip_v
  REAL(kind=rsh), PUBLIC,DIMENSION(imin:imax,jmin:jmax) :: hcentre,Uorb,Vorb,Uorbm,Vorbm,urms
  REAL(kind=rsh), PUBLIC,DIMENSION(imin:imax,jmin:jmax) :: Orb_u,Orb_v,Qx_u,Qy_v,Qy_u,Qx_v
  REAL(kind=rsh), PUBLIC,DIMENSION(imin:imax,jmin:jmax) :: TmBot,Hrms,Orbwc,Orbwt,UOrbwc,UOrbwt
  REAL(kind=rsh), PUBLIC,DIMENSION(imin:imax,jmin:jmax) :: VOrbwc,VOrbwt,radu,radv_temp,cw,rhoAsurT
  REAL(kind=rsh), PUBLIC,DIMENSION(imin:imax,jmin:jmax) :: radv,CfH,JSMITH,k_x,k_y,Dsmith
  REAL(kind=rsh), PUBLIC,DIMENSION(imin:imax,jmin:jmax) :: temp_radu,temp_radv
  REAL(kind=rsh), PUBLIC,DIMENSION(imin:imax,jmin:jmax) :: Cf,frot_FW,frot_CD,frot_FC 
  REAL(kind=rsh), PUBLIC,DIMENSION(imin:imax,jmin:jmax) :: TauW,TauC,TauM 
  REAL(kind=rsh), PUBLIC,DIMENSION(imin:imax,jmin:jmax) :: TaumW,TaumCW,TauCW 
  REAL(kind=rsh), PUBLIC,DIMENSION(imin:imax,jmin:jmax) :: uf,vf,phicw,dir_W,dir_C,phiwc 
  REAL(kind=rsh) :: za0,kappa,iCD,ifw

CONTAINS

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

  SUBROUTINE init_wave_forcing_fric

    !&E---------------------------------------------------------------------
    !&E                 ***  ROUTINE init_wave_forcing_fric  ***
    !&E
    !&E ** Purpose : Initialize the 2D fields entering the computation of surface friction.
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls :
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used
    USE comvars2d, ONLY : ig,id,dx_u,dy_v,jb,jh,htx,hty
    USE comcolig2d,ONLY : igdu,igdv

    !! * Arguments

    !! * Local declarations
    REAL(kind=rsh) :: som_rad1,som_rad2,som_rad3,som_rad4
    INTEGER        :: k,i,j

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

    ! RAZ wave

    Sxx(:,:)   =0.0_rsh
    Sxy(:,:)   =0.0_rsh
    Syy(:,:)   =0.0_rsh
    Qx(:,:)    =0.0_rsh
    Qx_u(:,:)  =0.0_rsh
    Qy(:,:)    =0.0_rsh
    Qy_v(:,:)  =0.0_rsh
    Orb(:,:)   =0.0_rsh
    Orb_u(:,:) =0.0_rsh
    Orb_v(:,:) =0.0_rsh
    Fradx(:,:) =0.0_rsh
    Frady(:,:) =0.0_rsh
    Uorb(:,:)  =0.0_rsh
    Vorb(:,:)  =0.0_rsh
    Uorbm(:,:) =0.0_rsh
    Vorbm(:,:) =0.0_rsh
    UOrbwc(:,:)=0.0_rsh
    UOrbwt(:,:)=0.0_rsh
    VOrbwc(:,:)=0.0_rsh
    VOrbwt(:,:)=0.0_rsh
    Orbwc(:,:) =0.0_rsh
    Orbwt(:,:) =0.0_rsh
    hrms(:,:)  =0.0_rsh
    CfH(:,:)   =0.0_rsh
    cw(:,:)    =0.0_rsh
    direc(:,:) =0.0_rsh
    wlen(:,:)  =0.0_rsh

    ! RAZ friction

    Cf(:,:)     =0.0_rsh
    frot_FW(:,:)=0.0_rsh
    frot_CD(:,:)=0.0_rsh
    frot_FC(:,:)=0.0_rsh
    TauW(:,:)   =0.0_rsh
    TauC(:,:)   =0.0_rsh
    TauM(:,:)   =0.0_rsh
    TaumW(:,:)  =0.0_rsh
    TaumCW(:,:) =0.0_rsh
    TauCW(:,:)  =0.0_rsh
    phicw(:,:)  =0.0_rsh
    dir_W(:,:)  =0.0_rsh

  END SUBROUTINE init_wave_forcing_fric

  !!----------------------------------------------------------------------

  SUBROUTINE lec_para_swan

    !&E---------------------------------------------------------------------
    !&E                 ***  ROUTINE lec_para_swan  ***
    !&E
    !&E ** Purpose : READ the wave parameters file.
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls :
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used

    IMPLICIT NONE

    !! * Arguments

    !! * Local declarations
    CHARACTER(len=lchain) :: fileps

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

    !----------------- Nico --------------------
    NAMELIST/nawbathymetrie/ lect_sens,fact_bathy,level
    NAMELIST/nawhoule/gammab,l_wave_brgm_roller,l_wave_brgm_asym,M_battjes,hlim_calc, &
         l_wave_brgm_phillips,l_wave_brgm_smith,iter_max,err_conv,l_stat_conv,l_wave_brgm_t2c
    NAMELIST/nawfrot/iCD,ifw,Kappa,za0,l_wave_brgm_frot_homog,frot_val
    !NAMELIST/nawmorpho/morpho,iformule,ieff,rhos,Te,d50hom,aphi,choix_ws,ordr_legendre, &
         !                   ischem_morpho,dt_morpho,poro_sedim,f_evolxyz,f_evolbot,mDW,nDW

    !----------------- Fin Nico ----------------

    fileps='./parawave.txt'
    OPEN(50,file=fileps,status='old',form='formatted', access='sequential')
    READ(50,nawbathymetrie)
    READ(50,nawhoule)
    READ(50,nawfrot)
    !READ(50,nawmorpho)
    CLOSE(50)

    CfH(:,:)=frot_val


  END SUBROUTINE lec_para_swan

  !!----------------------------------------------------------------------

  SUBROUTINE lect_wave

    !&E---------------------------------------------------------------------
    !&E                 ***  ROUTINE lect_wave  ***
    !&E
    !&E ** Purpose : READ the DATA files generated from the SWAN wave model.
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls :
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used

    !! * Arguments

    !! * Local declarations
    CHARACTER(LEN=lchain) :: fich_hs, fich_bat, fich_direc, fich_qb, fich_depth,test
    CHARACTER(LEN=lchain) :: fich_wlen, fich_per, fich_flux, fich_rad,nf, fich_tp, fich_tmbot

    INTEGER :: nf_hs, nf_bat, nf_direc, nf_qb, nf_depth,borne
    INTEGER :: i, j , nf_wlen, nf_per,indic, nf_frad, nf_tp, nf_tmbot ,iii,jj,eof

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

    ! READ SWAN FILES

    fich_hs    = './HOULE/hs.dat'
    fich_bat   = './HOULE/bathy.bot'
    fich_direc = './HOULE/dir.dat'
    fich_qb    = './HOULE/qb.dat'
    fich_depth = './HOULE/depth.dat'
    fich_wlen  = './HOULE/wlen.dat'
    fich_per   = './HOULE/tm0.dat'
    fich_rad   = './HOULE/Frad.dat'
    fich_tp    = './HOULE/tp.dat'
    fich_tmbot = './HOULE/tmbot.dat'

    ! Channels for reading
    !-------------------------------------------------------------------------

    nf_hs    = 201
    nf_bat   = 202
    nf_direc = 203
    nf_qb    = 204
    nf_depth = 205
    nf_wlen  = 206
    nf_per   = 207
    nf_frad  = 208
    nf_tp    = 209
    nf_tmbot = 210

    ! READING of frad.dat
    !-------------------------------------------------------------------------

    OPEN(UNIT=nf_frad , FILE=fich_rad, STATUS='old',form='formatted')
    DO j = jmax,jmin,-1
      READ(nf_frad,*,iostat=eof) Fradx(imin:imax,j)
    ENDDO
    DO j = jmax,jmin,-1
      READ(nf_frad,*,iostat=eof) Frady(imin:imax,j)
    ENDDO
    CLOSE(nf_frad)

    ! READING of hs.dat
    !-------------------------------------------------------------------------

    OPEN(UNIT=nf_hs , FILE=fich_hs , STATUS='old',form='formatted')
    DO j = jmax,jmin,-1
      READ(nf_hs,*,iostat=eof) Hs(imin:imax,j)
    ENDDO
    CLOSE(nf_hs)

    ! READING of direc.dat
    !-------------------------------------------------------------------------

    OPEN(UNIT=nf_direc , FILE=fich_direc , STATUS='old')
    DO j = jmax,jmin,-1
      READ(nf_direc,*,iostat=eof) direc(imin:imax,j)
    ENDDO
    CLOSE(nf_direc)

    ! READING of qb.dat
    !-------------------------------------------------------------------------

    OPEN(UNIT=nf_qb , FILE=fich_qb , STATUS='old')
    DO j = jmax,jmin,-1
      READ(nf_qb,*,iostat=eof) Qb(imin:imax,j)
    ENDDO
    CLOSE(nf_qb)

    ! READING of depth.dat
    !-------------------------------------------------------------------------

    OPEN(UNIT=nf_depth , FILE=fich_depth , STATUS='old')
    DO j = jmax,jmin,-1
      READ(nf_depth,*,iostat=eof) Depth(imin:imax,j)
    ENDDO
    CLOSE(nf_depth)

    ! READING of wlen.dat
    !-------------------------------------------------------------------------

    OPEN(UNIT=nf_wlen , FILE=fich_wlen , STATUS='old')
    DO j = jmax,jmin,-1
      READ(nf_wlen,*,iostat=eof) Wlen(imin:imax,j)
    ENDDO
    CLOSE(nf_wlen)

    ! READING of tm0.dat
    !-------------------------------------------------------------------------

    OPEN(UNIT=nf_per , FILE=fich_per , STATUS='old')
    DO j = jmax,jmin,-1
      READ(nf_per,*,iostat=eof) Tm01(imin:imax,j)
    ENDDO
    CLOSE(nf_per)

    ! READING of tp.dat
    !-------------------------------------------------------------------------

    OPEN(UNIT=nf_tp , FILE=fich_tp , STATUS='old')
    DO j = jmax,jmin,-1
      READ(nf_tp,*,iostat=eof) Tp(imin:imax,j)
    ENDDO
    CLOSE(nf_tp)

    ! READING of tmbot.dat
    !-------------------------------------------------------------------------

    OPEN(UNIT=nf_tmbot , FILE=fich_tp , STATUS='old')
    DO j = jmax,jmin,-1
      READ(nf_tmbot,*,iostat=eof) TmBot(imin:imax,j)
    ENDDO
    CLOSE(nf_tmbot)

    !-------------------------------------------------------------------------

    ! WRITE(*,*) 'Reading of  SWAN files              : OK'

  END SUBROUTINE lect_wave

  !!-------------------------------------------------------------------------

  SUBROUTINE compute_wave_forcing

    !&E---------------------------------------------------------------------
    !&E                 ***  ROUTINE compute_wave_forcing  ***
    !&E
    !&E ** Purpose : Compute the ocean forcing associated WITH waves.
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls :
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used
    USE comvars2d, ONLY : dx_u,dy_v,dx,dy,l_obc_cycl_x,h0,nbouc
    USE comvars3d, ONLY : rhoref

    !! * Arguments

    !! * Local declarations
    INTEGER :: i,j
    REAL(KIND=rsh),DIMENSION(imin:imax,jmin:jmax) :: HWATER
    CHARACTER(len=20)                             :: test
    ! DECLARATION DES TABLEAUX DE CALCUL
    REAL(KIND=rsh),DIMENSION(imin:imax,jmin:jmax) :: k            ! Nbre d ondes
    REAL(KIND=rsh),DIMENSION(imin:imax,jmin:jmax) :: cg           ! Vitesse de groupe
    REAL(KIND=rsh),DIMENSION(imin:imax,jmin:jmax) :: Energie      ! Energie des vagues
    REAL(KIND=rsh),DIMENSION(imin:imax,jmin:jmax) :: Qrx,Qry
    REAL(KIND=rsh),DIMENSION(imin:imax,jmin:jmax) :: Rxx,Rxy,Ryy
    ! DECLARATION POUR LE ROULEAU
    REAL(KIND=rsh)            :: sec_mbre,fact,moyFr,limFR
    ! PARAMETRES PHYSIQUES
    REAL(KIND=rsh), PARAMETER :: betad=0.2_rsh              ! Roller (Dally (2001))

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

    ! INITIALISATION VARIABLES INTERNES
    !-------------------------------------------------------------------------

    hwater(:,:)=MAX(h0(:,:),0.0001_rsh)


    ! CALCUL DES DIFFERENTS TERMES DE FORCAGE
    !-------------------------------------------------------------------------

    !              Correction de la hauteur d eau
    !-------------------------------------------------------------------------
    ! WHERE( HWATER(:,:) < hlim_calc ) HWATER(:,:)=0._rsh


    ! Calcul des premières variables caractéristiques
    !-------------------------------------------------------------------------
    Hrms(:,:)       = Hs(:,:)/SQRT(2.0_rsh)

    ! linear increase of waves during first ierations
    ! in order to avoid numerical oscillations (shock at t=0)
    IF (nbouc <= 10000) THEN
      Hrms(:,:) = Hrms(:,:)*REAL(nbouc,rsh)/10000.0_rsh
      PRINT*,'Nbouc :',nbouc,' -- ',REAL(100*nbouc,rsh)/10000.0_rsh,' % de Hrms'
    END IF
!cval Vous voulez vraiment mettre la meme variable (le MAX) dans un tableau ????
    !k(:,:)          = MAX(2.0_rsh*pi/Wlen(:,:),0.00001_rsh)
    ! modif pour gfortran
    k(:,:)          = 2.0_rsh*pi/Wlen(:,:)
    WHERE(k(:,:)<0.00001_rsh) k(:,:)=0.00001_rsh
    cw(:,:)         = SQRT( (grav/k(:,:)) * TANH(k(:,:)*HWATER(:,:)) )
    cg(:,:)         = 0.5_rsh * cw(:,:) * ( 1.0_rsh + 2.0_rsh*k(:,:)*HWATER(:,:)/SINH(2.0_rsh*k(:,:)*HWATER(:,:)) )
    Energie(:,:)    = rhoref * grav * Hrms(:,:)**2 / 8.0_rsh
    Qx(:,:)         = Energie(:,:) / (rhoref*cw(:,:)) * COS(Direc(:,:)*pi/180.0_rsh)
    Qy(:,:)         = Energie(:,:) / (rhoref*cw(:,:)) * SIN(Direc(:,:)*pi/180.0_rsh)

    WHERE (Tm01(:,:) <= 0.0001_rsh)
      cw(:,:)     =0.0_rsh
      cg(:,:)     =0.0_rsh
      k(:,:)      =0.0_rsh
      Qx(:,:)     =0.0_rsh
      Qy(:,:)     =0.0_rsh
      energie(:,:)=0.0_rsh
      Hrms(:,:)   =0.0_rsh
    END WHERE

    WRITE(*,*) 'Calcul des variables caracteristiques     : OK'

    ! Calcul des vitesses orbitales
    !-------------------------------------------------------------------------
    !k(:,:)    = MAX(2.0_rsh*pi/Wlen(:,:),0.00001_rsh)
    ! modif pour gfortran
    k(:,:)          = 2.0_rsh*pi/Wlen(:,:)
    WHERE(k(:,:)<0.00001_rsh) k(:,:)=0.00001_rsh
    Tm01(:,:) = MAX(Tm01(:,:),0.00001_rsh)
    cw(:,:)   = MAX(cw(:,:),0.00001_rsh)
    Uorb(:,:) = pi * Hrms(:,:)/ (Tm01(:,:)*SINH(k(:,:)*HWATER(:,:))) * COS(Direc(:,:)*pi/180.0_rsh)
    Vorb(:,:) = pi * Hrms(:,:)/ (Tm01(:,:)*SINH(k(:,:)*HWATER(:,:))) * SIN(Direc(:,:)*pi/180.0_rsh)
    Orb(:,:)  = SQRT( Uorb(:,:)**2 + Vorb(:,:)**2 )
    Urms(:,:) = SQRT( grav * k(:,:) * TANH(k(:,:)*HWATER(:,:)) )* Hrms(:,:)/(2.0_rsh*SINH(k(:,:)*HWATER(:,:)))

    ! WRITE(*,*) 'Calcul des vitesses orbitales             : OK'

    ! Calcul de la dissipation d energie
    !-------------------------------------------------------------------------
    Dissip(:,:) = Qb(:,:) * rhoref * grav * ( gammab * HWATER(:,:) )**2 / (4.0_rsh*Tp(:,:))
    WHERE (Dissip(:,:) < 0.0_rsh) Dissip(:,:) = 0.0_rsh

    ! WRITE(*,*) 'Calcul de la dissipation          : OK'

    ! Calcul des caracteristiques du Roller
    !-------------------------------------------------------------------------
    IF (l_wave_brgm_roller) THEN
      rhoAsurT(:,:)= 0.0_rsh
      DO i=imin+1,imax
        DO j=jmin+1,jmax
          fact  = cw(i,j)**2*COS(Direc(i,j)*pi/180.0_rsh)/dx(i,j) &
               + cw(i,j)**2*SIN(Direc(i,j)*pi/180.0_rsh)/dy(i,j) &
               + grav * betad

          sec_mbre = Dissip(i,j) &
               + cw(i-1,j)**2*COS(Direc(i-1,j)*pi/180.0_rsh)/dx(i-1,j) * rhoAsurT(i-1,j) &
               + cw(i,j-1)**2*SIN(Direc(i,j-1)*pi/180.0_rsh)/dy(i,j-1) * rhoAsurT(i,j-1)
          rhoAsurT(i,j) = sec_mbre / fact

          IF(rhoAsurT(i,j) <= 0.0_rsh) rhoAsurT(i,j)=0.0_rsh
          IF(Qb(i,j)==0.0_rsh) rhoAsurT(i,j)=0.0_rsh
        ENDDO
      ENDDO

      IF(     l_obc_cycl_x) CALL LIM_PERIO_X(rhoAsurT)
      IF(.NOT.l_obc_cycl_x) CALL LIM_LIBRE_X(rhoAsurT)
      CALL LIM_LIBRE_Y(rhoAsurT)

      Qrx(:,:)= rhoAsurT(:,:) * COS(Direc(:,:)*pi/180.0_rsh)/rhoref
      Qry(:,:)= rhoAsurT(:,:) * SIN(Direc(:,:)*pi/180.0_rsh)/rhoref
      Rxx(:,:)= cw(:,:) * rhoAsurT(:,:) * (COS(Direc(:,:)*pi/180.0_rsh)**2+0.5_rsh)
      Ryy(:,:)= cw(:,:) * rhoAsurT(:,:) * (SIN(Direc(:,:)*pi/180.0_rsh)**2+0.5_rsh)
      Rxy(:,:)= cw(:,:) * rhoAsurT(:,:) * SIN(Direc(:,:)*pi/180.0_rsh)*COS(Direc(:,:)*pi/180.0_rsh)
      Qx(:,:) = Qx(:,:) + Qrx(:,:)
      Qy(:,:) = Qy(:,:) + Qry(:,:)

      ! WRITE(*,*) 'Calcul du    Roller                     : OK'
    ENDIF

    ! Calcul du J, theorie de SMITH 2006
    !-------------------------------------------------------------------------

    IF(.NOT.l_wave_brgm_phillips .AND. l_wave_brgm_smith) THEN
      JSMITH(:,:) = (cg(:,:)/cw(:,:)-0.5_rsh)*Energie(:,:)/HWATER(:,:)/rhoref
      k_x(:,:) = k(:,:) * COS(Direc(:,:)*pi/180.0_rsh)
      k_y(:,:) = k(:,:) * SIN(Direc(:,:)*pi/180.0_rsh)
      Dsmith(:,:) = Dissip(:,:)/k(:,:)/cw(:,:)
      IF (l_wave_brgm_roller) THEN
        JSMITH(:,:) = JSMITH(:,:) + cw(:,:)*rhoAsurT(:,:)/HWATER(:,:)/rhoref/2.0_rsh
        Dsmith(:,:) = grav*betad*rhoAsurT(:,:)/k(:,:)/cw(:,:)
      ENDIF
    ENDIF

    ! Calcul des tensions de radiations
    !-------------------------------------------------------------------------
    IF(.NOT.l_wave_brgm_phillips .AND. .NOT.l_wave_brgm_smith) THEN
      ! WRITE(*,*) 'MODELE DE MEI (1989)'
      Sxx(:,:) = 0.5_rsh * Energie(:,:)  * ( (COS(Direc(:,:)*pi/180.0_rsh)**2 + 1.0_rsh) * 2.0_rsh * cg(:,:) &
                                              / cw(:,:) - 1.0_rsh ) ! - rho * Qx(:,:)*Qx(:,:)/HWATER(:,:) ??
      Syy(:,:) = 0.5_rsh * Energie(:,:)  * ( (SIN(Direc(:,:)*pi/180.0_rsh)**2 + 1.0_rsh) * 2.0_rsh * cg(:,:) &
                                              / cw(:,:) - 1.0_rsh ) ! - rho * Qy(:,:)*Qy(:,:)/HWATER(:,:) ??
      Sxy(:,:) = 0.5_rsh * Energie(:,:)  * SIN(2.0_rsh*Direc(:,:)*pi/180.0_rsh) * cg(:,:) / cw(:,:) ! - rho * Qx(:,:)*Qy(:,:)/HWATER(:,:) ??
      IF (l_wave_brgm_roller) THEN
        Sxx(:,:) = Sxx(:,:) + Rxx(:,:)
        Syy(:,:) = Syy(:,:) + Ryy(:,:)
        Sxy(:,:) = Sxy(:,:) + Rxy(:,:)
      ENDIF
    ENDIF

    IF(l_wave_brgm_phillips .AND. .NOT.l_wave_brgm_smith) THEN
      ! WRITE(*,*) 'MODELE DE PHILLIPS (1977)'
      Sxx(:,:) = Energie(:,:) / 2.0_rsh * ( (COS(Direc(:,:)*pi/180.0_rsh)**2 + 1.0_rsh) * 2.0_rsh * cg(:,:) &
                                             / cw(:,:) - 1.0_rsh ) ! - rho * Qx(:,:)*Qx(:,:)/HWATER(:,:)
      Syy(:,:) = Energie(:,:) / 2.0_rsh * ( (SIN(Direc(:,:)*pi/180.0_rsh)**2 + 1.0_rsh) * 2.0_rsh * cg(:,:) &
                                             / cw(:,:) - 1.0_rsh ) ! - rho * Qy(:,:)*Qy(:,:)/HWATER(:,:)
      Sxy(:,:) = Energie(:,:) / 2.0_rsh * SIN(2.0_rsh*Direc(:,:)*pi/180.0_rsh) * cg(:,:) / cw(:,:) ! - rho * Qx(:,:)*Qy(:,:)/HWATER(:,:)
      IF (l_wave_brgm_roller) THEN
        Sxx(:,:) = Sxx(:,:) + Rxx(:,:)
        Syy(:,:) = Syy(:,:) + Ryy(:,:)
        Sxy(:,:) = Sxy(:,:) + Rxy(:,:)
      ENDIF
    ENDIF

    IF(.NOT.l_wave_brgm_phillips .AND. l_wave_brgm_smith) THEN
      ! WRITE(*,*) 'MODELE DE SMITH (2006)'
      Sxx(:,:) = 0.0_rsh
      Syy(:,:) = 0.0_rsh
      Sxy(:,:) = 0.0_rsh
    ENDIF

    IF(l_wave_brgm_phillips .AND. l_wave_brgm_smith) THEN
      WRITE(*,*) 'ANORMAL : ON ARRETE LE CALCUL'
      STOP
    ENDIF

    ! WRITE(*,*) 'Calcul des tensions de radiation  : OK'

    ! Mise a jour au niveau du trait de cote
    !-------------------------------------------------------------------------

    WHERE (Tm01(:,:) <= 0.0001_rsh)
      Sxx(:,:)   = 0.0_rsh
      Sxy(:,:)   = 0.0_rsh
      Syy(:,:)   = 0.0_rsh
      Orb(:,:)   = 0.0_rsh
      Uorb(:,:)  = 0.0_rsh
      Vorb(:,:)  = 0.0_rsh
      JSMITH(:,:)= 0.0_rsh
      k_x(:,:)   = 0.0_rsh
      k_y(:,:)   = 0.0_rsh
      Dissip(:,:)= 0.0_rsh
      Dsmith(:,:)= 0.0_rsh
    ENDWHERE
    WHERE (Hwater(:,:) <= hlim_calc)
      Sxx(:,:)   = 0.0_rsh
      Sxy(:,:)   = 0.0_rsh
      Syy(:,:)   = 0.0_rsh
      Qx(:,:)    = 0.0_rsh
      Qy(:,:)    = 0.0_rsh
      Orb(:,:)   = 0.0_rsh
      Uorb(:,:)  = 0.0_rsh
      Vorb(:,:)  = 0.0_rsh
      JSMITH(:,:)= 0.0_rsh
      k_x(:,:)   = 0.0_rsh
      k_y(:,:)   = 0.0_rsh
      Dissip(:,:)= 0.0_rsh
      Dsmith(:,:)= 0.0_rsh
    ENDWHERE

    ! Calcul des variables en u et en v
    !-------------------------------------------------------------------------
    Qx_u   = MIDDLEX(Qx)
    Qy_v   = MIDDLEY(Qy)
    Qy_u   = MIDDLEX(Qy)
    Qx_v   = MIDDLEY(Qx)

    Orb_u  = MIDDLEX(Orb)
    Orb_v  = MIDDLEY(Orb)

    Dissip_u = MIDDLEX(Dissip)
    Dissip_v = MIDDLEY(Dissip)

    WHERE (Hwater(:,:) <= hlim_calc)
      Qx_u(:,:)   = 0.0_rsh
      Qy_v(:,:)   = 0.0_rsh
      Qy_u(:,:)   = 0.0_rsh
      Qx_v(:,:)   = 0.0_rsh
      Orb_u(:,:)  = 0.0_rsh
      Orb_v(:,:)  = 0.0_rsh
    ENDWHERE
    ! CALCUL DES FORCES
    !-------------------------------------------------------------------------
    radu(:,:)     = 0.0_rsh
    radv_temp(:,:)= 0.0_rsh
    IF(.NOT.l_wave_brgm_smith)THEN
      DO j=jmin+1,jmax-1
        DO i=imin+1,imax-1
          radu(i,j) = ( Sxx(i+1,j  )-Sxx(i,j    ) ) / dx_u(i,j) + &
               ( ( Sxy(i,j+1  )-Sxy(i,j-1  ) ) /(dy_v(i,j-1  )+dy_v(i,j  )) &
               + ( Sxy(i+1,j+1)-Sxy(i+1,j-1) ) /(dy_v(i+1,j-1)+dy_v(i+1,j)) ) /2.0_rsh

          radv_temp(i,j)= ( Syy(i,j+1  )-Syy(i,j    ) ) / dy_v(i,j) + &
               ( ( Sxy(i+1,j  )-Sxy(i-1,j  ) ) /(dx_u(i-1,j  )+dx_u(i,j  )) &
               + ( Sxy(i+1,j+1)-Sxy(i-1,j+1) ) /(dx_u(i-1,j+1)+dx_u(i,j+1)) ) /2.0_rsh

        ENDDO
      ENDDO
      !---------
      ! DECOMMENTER C 2 LIGNES POUR UTILISER LES TENSIONS DE RADIATION DE SWAN
      ! radu     (:,:) = Fradx(:,:)
      ! radv_temp(:,:) = Frady(:,:)
      !---------

      IF(l_obc_cycl_x) CALL LIM_PERIO_X(radu)
      IF(l_obc_cycl_x) CALL LIM_PERIO_X(radv_temp)
      IF(.NOT.l_obc_cycl_x) CALL LIM_LIBRE_X(radu)
      IF(.NOT.l_obc_cycl_x) CALL LIM_LIBRE_X(radv_temp)
      CALL LIM_LIBRE_Y(radu)
      CALL LIM_LIBRE_Y(radv_temp)

      WHERE (Hrms(:,:) .LE. 0.0_rsh) radu(:,:)=0.0_rsh
      WHERE (Hrms(:,:) .LE. 0.0_rsh) radv_temp(:,:)=0.0_rsh

      radv(:,:)=radv_temp(:,:)
    ENDIF

  END SUBROUTINE compute_wave_forcing

  !!----------------------------------------------------------------------

  SUBROUTINE div_tens_rad

    !&E---------------------------------------------------------------------
    !&E                 ***  ROUTINE div_tens_rad  ***
    !&E
    !&E ** Purpose : 
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls :
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used
    USE comvars2d, ONLY : ig,id,dx_u,dy_v,jb,jh,htx,hty
    USE comcolig2d,ONLY : igdu,igdv

    !! * Arguments

    !! * Local declarations
    REAL(kind=rsh) :: som_rad1,som_rad2,som_rad3,som_rad4
    INTEGER        :: k,i,j

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

    IF(l_wave_brgm_t2c) THEN
      DO j=jmin+2,jmax-1
        DO i = ig(j),id(j)-1
          som_rad1 = 0.0_rsh
          som_rad2 = 0.0_rsh
          som_rad3 = 0.0_rsh
          som_rad4 = 0.0_rsh
          temp_radu(i,j) = ( Sxx(i+1,j)-Sxx(i,j) ) /dx_u(i,j)

          IF (igdu(i,j) == 4 .OR. htx(i,j) <= hlim_calc) THEN  ! a priori ce test est un pleonasme.
            temp_radu(i,j) = 0.0_rsh
          ELSE
            IF (igdu(i,j+1) == 2 .OR. igdu(i,j+1) == 1) som_rad1 = ( Sxy(i,j+1)-Sxy(i,j) ) /dy_v(i,j)
            IF (igdu(i,j-1) == 2 .OR. igdu(i,j-1) == 1) som_rad2 = ( Sxy(i,j)-Sxy(i,j-1) ) /dy_v(i,j-1)
            IF (igdu(i,j+1) == 3 .OR. igdu(i,j+1) == 1) som_rad3 = ( Sxy(i+1,j+1)-Sxy(i+1,j) ) /dy_v(i+1,j)
            IF (igdu(i,j-1) == 3 .OR. igdu(i,j-1) == 1) som_rad4 = ( Sxy(i+1,j)-Sxy(i+1,j-1) ) /dy_v(i+1,j-1)
            IF (igdu(i,j-1) == 1 .AND. igdu(i,j+1) == 1) &
                 temp_radu(i,j) = temp_radu(i,j) + ( (som_rad1+som_rad2)/2.0_rsh &
                                                   + (som_rad3+som_rad4)/2.0_rsh )/2.0_rsh
            IF (igdu(i,j-1) == 1 .AND. igdu(i,j+1) == 2) &
                 temp_radu(i,j) = temp_radu(i,j) + ( (som_rad1+som_rad2)/2.0_rsh + som_rad4 )/2.0_rsh
            IF (igdu(i,j-1) == 1 .AND. igdu(i,j+1) == 3) &
                 temp_radu(i,j) = temp_radu(i,j) + ( som_rad2 + (som_rad3+som_rad4)/2.0_rsh )/2.0_rsh
            IF (igdu(i,j-1) == 1 .AND. igdu(i,j+1) == 4) &
                 temp_radu(i,j) = temp_radu(i,j) + ( som_rad2 + som_rad4)/2.0_rsh
            IF (igdu(i,j-1) == 2 .AND. igdu(i,j+1) == 1) &
                 temp_radu(i,j) = temp_radu(i,j) + ( (som_rad1+ som_rad2)/2.0_rsh + som_rad3 )/2.0_rsh
            IF (igdu(i,j-1) == 2 .AND. igdu(i,j+1) == 2) &
                 temp_radu(i,j) = temp_radu(i,j) + ( som_rad1 + som_rad2 )/4.0_rsh
            IF (igdu(i,j-1) == 2 .AND. igdu(i,j+1) == 3) &
                 temp_radu(i,j) = temp_radu(i,j) + ( som_rad2 + som_rad3 )/2.0_rsh
            IF (igdu(i,j-1) == 2 .AND. igdu(i,j+1) == 4) &
                 temp_radu(i,j) = temp_radu(i,j) + som_rad2 /2.0_rsh
            IF (igdu(i,j-1) == 3 .AND. igdu(i,j+1) == 1) &
                 temp_radu(i,j) = temp_radu(i,j) + ( som_rad1 + (som_rad3+som_rad4)/2.0_rsh )/2.0_rsh
            IF (igdu(i,j-1) == 3 .AND. igdu(i,j+1) == 2) &
                 temp_radu(i,j) = temp_radu(i,j) + ( som_rad1 + som_rad4 )/2.0_rsh
            IF (igdu(i,j-1) == 3 .AND. igdu(i,j+1) == 3) &
                 temp_radu(i,j) = temp_radu(i,j) + ( som_rad3 + som_rad4 )/4.0_rsh
            IF (igdu(i,j-1) == 3 .AND. igdu(i,j+1) == 4) &
                 temp_radu(i,j) = temp_radu(i,j) + som_rad4 /2.0_rsh
            IF (igdu(i,j-1) == 4 .AND. igdu(i,j+1) == 1) &
                 temp_radu(i,j) = temp_radu(i,j) + ( som_rad1 + som_rad3 )/2.0_rsh
            IF (igdu(i,j-1) == 4 .AND. igdu(i,j+1) == 2) &
                 temp_radu(i,j) = temp_radu(i,j) + som_rad1 /2.0_rsh
            IF (igdu(i,j-1) == 4 .AND. igdu(i,j+1) == 3) &
                 temp_radu(i,j) = temp_radu(i,j) + som_rad3 /2.0_rsh
            IF (igdu(i,j-1) == 4 .AND. igdu(i,j+1) == 4) &
                 temp_radu(i,j) = temp_radu(i,j)
          END IF
        END DO
      END DO
      !
      DO j=jmin+1,jmax
        DO i =imin+2,imax-1
          IF(j >= jb(i) .AND.j <= jh(i)-1) THEN
            som_rad1 = 0.0_rsh
            som_rad2 = 0.0_rsh
            som_rad3 = 0.0_rsh
            som_rad4 = 0.0_rsh
            temp_radv(i,j) = ( Syy(i,j+1)-Syy(i,j) ) /dy_v(i,j)

            IF(igdv(i,j) == 4 .OR. hty(i,j) <= hlim_calc) THEN ! idem pour le pleonasme.
              temp_radv(i,j)=0.0_rsh
            ELSE
              IF (igdv(i+1,j) == 2 .OR. igdv(i+1,j) == 1) som_rad1 = ( Sxy(i+1,j)-Sxy(i,j) ) /dx_u(i,j)
              IF (igdv(i-1,j) == 2 .OR. igdv(i-1,j) == 1) som_rad2 = ( Sxy(i,j)-Sxy(i-1,j) ) /dx_u(i-1,j)
              IF (igdv(i+1,j) == 3 .OR. igdv(i+1,j) == 1) som_rad3 = ( Sxy(i+1,j+1)-Sxy(i,j+1) ) /dx_u(i,j+1)
              IF (igdv(i-1,j) == 3 .OR. igdv(i-1,j) == 1) som_rad4 = ( Sxy(i,j+1)-Sxy(i-1,j+1) ) /dx_u(i-1,j+1)

              IF (igdv(i-1,j) == 1 .AND. igdv(i+1,j) == 1) &
                   temp_radv(i,j) = temp_radv(i,j) + ( (som_rad1+som_rad2)/2.0_rsh &
                                                     + (som_rad3+som_rad4)/2.0_rsh )/2.0_rsh
              IF (igdv(i-1,j) == 1 .AND. igdv(i+1,j) == 2) &
                   temp_radv(i,j) = temp_radv(i,j) + ( (som_rad1+som_rad2)/2.0_rsh + som_rad4 )/2.0_rsh
              IF (igdv(i-1,j) == 1 .AND. igdv(i+1,j) == 3) &
                   temp_radv(i,j) = temp_radv(i,j) + ( som_rad2 + (som_rad3+som_rad4)/2.0_rsh )/2.0_rsh
              IF (igdv(i-1,j) == 1 .AND. igdv(i+1,j) == 4) &
                   temp_radv(i,j) = temp_radv(i,j) + ( som_rad2 + som_rad4 )/2.0_rsh
              IF (igdv(i-1,j) == 2 .AND. igdv(i+1,j) == 1) &
                   temp_radv(i,j) = temp_radv(i,j) + ( (som_rad1+som_rad2)/2.0_rsh + som_rad3 )/2.0_rsh
              IF (igdv(i-1,j) == 2 .AND. igdv(i+1,j) == 2) &
                   temp_radv(i,j) = temp_radv(i,j) + ( som_rad1 + som_rad2 )/4.0_rsh
              IF (igdv(i-1,j) == 2 .AND. igdv(i+1,j) == 3) &
                   temp_radv(i,j) = temp_radv(i,j) + ( som_rad2 + som_rad3 )/2.0_rsh
              IF (igdv(i-1,j) == 2 .AND. igdv(i+1,j) == 4) &
                   temp_radv(i,j) = temp_radv(i,j) + som_rad2 /2.0_rsh
              IF (igdv(i-1,j) == 3 .AND. igdv(i+1,j) == 1) &
                   temp_radv(i,j) = temp_radv(i,j) + ( som_rad1 + (som_rad3+som_rad4)/2.0_rsh )/2.0_rsh
              IF (igdv(i-1,j) == 3 .AND. igdv(i+1,j) == 2) &
                   temp_radv(i,j) = temp_radv(i,j) + ( som_rad1 + som_rad4 )/2.0_rsh
              IF (igdv(i-1,j) == 3 .AND. igdv(i+1,j) == 3) &
                   temp_radv(i,j) = temp_radv(i,j) + ( som_rad3 + som_rad4 )/4.0_rsh
              IF (igdv(i-1,j) == 3 .AND. igdv(i+1,j) == 4) &
                   temp_radv(i,j) = temp_radv(i,j) + som_rad4 /2.0_rsh
              IF (igdv(i-1,j) == 4 .AND. igdv(i+1,j) == 1) &
                   temp_radv(i,j) = temp_radv(i,j) + ( som_rad1 + som_rad3 )/2.0_rsh
              IF (igdv(i-1,j) == 4 .AND. igdv(i+1,j) == 2) &
                   temp_radv(i,j) = temp_radv(i,j) + som_rad1 /2.0_rsh
              IF (igdv(i-1,j) == 4 .AND. igdv(i+1,j) == 3) &
                   temp_radv(i,j) = temp_radv(i,j) + som_rad3 /2.0_rsh
              IF (igdv(i-1,j) == 4 .AND. igdv(i+1,j) == 4) &
                   temp_radv(i,j) = temp_radv(i,j)

            END IF

          END IF
        END DO
      END DO

    END IF

  END SUBROUTINE div_tens_rad

  !!----------------------------------------------------------------------

  SUBROUTINE LIM_PERIO_X(Var)

    !&E---------------------------------------------------------------------
    !&E                 ***  ROUTINE LIM_PERIO_X  ***
    !&E
    !&E ** Purpose : 
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls :
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used
    IMPLICIT NONE

    !! * Arguments
    REAL(KIND=rsh),DIMENSION(imin:imax,jmin:jmax),INTENT(inout) :: Var

    !! * Local declarations

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

    Var(imin  ,:) = Var(imax-2,:)
    Var(imin+1,:) = Var(imax-1,:)
    Var(imax  ,:) = Var(imin+2,:)

  END SUBROUTINE LIM_PERIO_X

  !!----------------------------------------------------------------------

  SUBROUTINE LIM_LIBRE_X(Var)

    !&E---------------------------------------------------------------------
    !&E                 ***  ROUTINE LIM_LIBRE_X  ***
    !&E
    !&E ** Purpose : 
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls :
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used
    IMPLICIT NONE

    !! * Arguments
    REAL(KIND=rsh),DIMENSION(imin:imax,jmin:jmax),INTENT(inout) :: Var

    !! * Local declarations

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

    Var(imin+1,:) = Var(imin+2,:)
    Var(imin,:)   = Var(imin+1,:)
    Var(imax,:)   = Var(imax-1,:)

  END SUBROUTINE LIM_LIBRE_X

  !!----------------------------------------------------------------------

  SUBROUTINE LIM_LIBRE_Y(Var)

    !&E---------------------------------------------------------------------
    !&E                 ***  ROUTINE LIM_LIBRE_Y  ***
    !&E
    !&E ** Purpose : 
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls :
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used
    IMPLICIT NONE

    !! * Arguments
    REAL(KIND=rsh),DIMENSION(imin:imax,jmin:jmax),INTENT(inout) :: Var

    !! * Local declarations

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

    Var(:,jmin+1) = Var(:,jmin+2)
    Var(:,jmin)   = Var(:,jmin+1)
    Var(:,jmax)   = Var(:,jmax-1)

  END SUBROUTINE LIM_LIBRE_Y

  !!----------------------------------------------------------------------

  SUBROUTINE LIM_NULLE_X(Var)
    !&E---------------------------------------------------------------------
    !&E                 ***  ROUTINE LIM_NULLE_X  ***
    !&E
    !&E ** Purpose : 
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls :
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used
    USE comvars2d, ONLY : htx
    IMPLICIT NONE

    !! * Arguments
    REAL(KIND=rsh),DIMENSION(imin:imax,jmin:jmax),INTENT(inout) :: Var

    !! * Local declarations

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

!cval attention la variable htx a ete allouee
    WHERE( htx(:,:) <= hlim_calc ) Var(:,:)=0.0_rsh

  END SUBROUTINE LIM_NULLE_X

  !!----------------------------------------------------------------------

  SUBROUTINE LIM_NULLE_Y(Var)

    !&E---------------------------------------------------------------------
    !&E                 ***  ROUTINE LIM_NULLE_Y  ***
    !&E
    !&E ** Purpose : 
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls :
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used
    USE comvars2d, ONLY : hty
    IMPLICIT NONE

    !! * Arguments
    REAL(KIND=rsh),DIMENSION(imin:imax,jmin:jmax),INTENT(inout) :: Var

    !! * Local declarations

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

!cval attention la variable hty a ete allouee
    WHERE( hty(:,:) <= hlim_calc ) Var(:,:)=0.0_rsh

  END SUBROUTINE LIM_NULLE_Y

  !!----------------------------------------------------------------------

  FUNCTION MIDDLEX(Var)

    !&E---------------------------------------------------------------------
    !&E                 ***  FUNCTION MIDDLEX  ***
    !&E
    !&E ** Purpose : 
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls :
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used

    USE comvars2d,ONLY : l_obc_cycl_x
    IMPLICIT NONE

    !! * Arguments
    REAL(kind=rsh),DIMENSION(imin:imax,jmin:jmax):: Var, MIDDLEX

    !! * Local declarations
    INTEGER :: i,j

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

    ! milieu en u de points sur les mailles scalaires
    !------------------------------

    DO j=jmin,jmax
      DO i=imin,imax-1
        MIDDLEX(i,j) = ( Var(i+1,j)+Var(i,j) ) / 2.0_rsh
      ENDDO
    ENDDO

    IF (l_obc_cycl_x) THEN
      MIDDLEX(imax,:)=MIDDLEX(imin+2,:)
    ELSE
      MIDDLEX(imax,:)=MIDDLEX(imax-1,:)
    ENDIF

  END FUNCTION MIDDLEX

  !!----------------------------------------------------------------------

  FUNCTION MIDDLEY(Var)

    !&E---------------------------------------------------------------------
    !&E                 ***  FUNCTION MIDDLEY  ***
    !&E
    !&E ** Purpose : 
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls :
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used
    IMPLICIT NONE

    !! * Arguments
    REAL(kind=rsh),DIMENSION(imin:imax,jmin:jmax):: Var, MIDDLEY

    !! * Local declarations
    INTEGER :: i,j

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

    ! milieu en v de points sur les mailles scalaires
    !------------------------------

    DO i=imin,imax
      DO j=jmin,jmax-1
        MIDDLEY(i,j) = ( Var(i,j+1)+Var(i,j) ) / 2.0_rsh
      ENDDO
    ENDDO

    MIDDLEY(:,jmax)=MIDDLEY(:,jmax-1)

  END FUNCTION MIDDLEY

  !!----------------------------------------------------------------------

  FUNCTION MIDDLEU(Var)

    !&E---------------------------------------------------------------------
    !&E                 ***  FUNCTION MIDDLEU  ***
    !&E
    !&E ** Purpose : 
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls :
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used
    USE comvars2d,ONLY : l_obc_cycl_x
    IMPLICIT NONE

    !! * Arguments
    REAL(kind=rsh),DIMENSION(imin:imax,jmin:jmax):: Var, MIDDLEU

    !! * Local declarations
    INTEGER :: i,j

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

    ! milieu au noeuds de points sur les demi-mailles u
    !-------------------------------

    DO i=imin+1,imax
      DO j=jmin,jmax
        MIDDLEU(i,j) = ( Var(i-1,j)+Var(i,j) ) / 2.0_rsh
      ENDDO
    ENDDO

    IF (l_obc_cycl_x) THEN
      MIDDLEU(imin,:)=MIDDLEU(imax-2,:)
    ELSE
      MIDDLEU(imin,:)=MIDDLEU(imin+1,:)
    ENDIF

  END FUNCTION MIDDLEU

  !!----------------------------------------------------------------------

  FUNCTION MIDDLEV(Var)

    !&E---------------------------------------------------------------------
    !&E                 ***  FUNCTION MIDDLEV  ***
    !&E
    !&E ** Purpose : 
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls :
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used
    IMPLICIT NONE

    !! * Arguments
    REAL(kind=rsh),DIMENSION(imin:imax,jmin:jmax):: Var, MIDDLEV

    !! * Local declarations
    INTEGER :: i,j

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

    ! milieu au noeuds de points sur les demi-mailles v
    !-------------------------------

    DO i=imin,imax
      DO j=jmin+1,jmax
        MIDDLEV(i,j) = ( Var(i,j-1)+Var(i,j) ) / 2.0_rsh
      ENDDO
    ENDDO

    MIDDLEV(:,jmin)=MIDDLEV(:,jmin+1)

  END FUNCTION MIDDLEV

  !!----------------------------------------------------------------------

  FUNCTION MIDDLEN(VarX,VarY)

    !&E---------------------------------------------------------------------
    !&E                 ***  FUNCTION MIDDLEN  ***
    !&E
    !&E ** Purpose : 
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls :
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used
    USE comvars2d,ONLY : l_obc_cycl_x
    IMPLICIT NONE

    !! * Arguments
    REAL(kind=rsh),DIMENSION(imin:imax,jmin:jmax):: VarX, VarY, MIDDLEN

    !! * Local declarations
    INTEGER :: i,j

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

    ! milieu au noeuds N de points sur les demi-mailles v et u
    !---------------------------------

    DO i=imin+1,imax
      DO j=jmin+1,jmax
        MIDDLEN(i,j) = ( VarY(i,j-1)+VarY(i,j)+VarX(i-1,j)+VarX(i,j) ) / 4.0_rsh
      ENDDO
    ENDDO

    MIDDLEN(:,jmin)=MIDDLEN(:,jmin+1)
    IF (l_obc_cycl_x) THEN
      MIDDLEN(imin,:)=MIDDLEN(imax-2,:)
    ELSE
      MIDDLEN(imin,:)=MIDDLEN(imin+1,:)
    ENDIF

  END FUNCTION MIDDLEN

  !!----------------------------------------------------------------------

  FUNCTION DNDX2(Var)

    !&E---------------------------------------------------------------------
    !&E                 ***  FUNCTION DNDX2  ***
    !&E
    !&E ** Purpose : Compute 1st order derivative in X direction.
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls : LIM_PERIO_X, LIM_LIBRE_X
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used
    USE comvars2d,ONLY : l_obc_cycl_x,dx_u
    IMPLICIT NONE

    !! * Arguments
    REAL(kind=rsh),DIMENSION(imin:imax,jmin:jmax):: Var, DNDX2

    !! * Local declarations
    INTEGER :: i,j

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

    ! Derivee d une variable sur un noeud scalaire dans la direction x
    ! schema centre ordre 2
    !------------------------------------

    DO j=jmin,jmax
      DO i=imin+1,imax-1
        DNDX2(i,j) = ( Var(i+1,j)-Var(i-1,j) ) / ( dx_u(i,j)+dx_u(i-1,j))
      ENDDO
    ENDDO

    IF (l_obc_cycl_x) THEN
      CALL LIM_PERIO_X(DNDX2)
    ELSE
      CALL LIM_LIBRE_X(DNDX2)
    ENDIF

  END FUNCTION DNDX2

  !!----------------------------------------------------------------------

  FUNCTION DNDY2(Var)

    !&E---------------------------------------------------------------------
    !&E                 ***  FUNCTION DNDY2  ***
    !&E
    !&E ** Purpose : Compute 1st order derivative in Y direction.
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls : LIM_LIBRE_Y
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used
    USE comvars2d,ONLY : dy
    IMPLICIT NONE

    !! * Arguments
    REAL(kind=rsh),DIMENSION(imin:imax,jmin:jmax):: Var, DNDY2

    !! * Local declarations
    INTEGER :: i,j

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

    ! Derivee d une variable sur un noeud scalaire dans la direction y
    ! schema centre ordre 2
    !------------------------------------
    DO i=imin,imax
      DO j=jmin+1,jmax-1
        DNDY2(i,j) = ( Var(i,j+1)-Var(i,j-1) ) / (dy(i,j)+dy(i,j-1))
      ENDDO
    ENDDO

    CALL LIM_LIBRE_Y(DNDY2)

  END FUNCTION DNDY2

  !!----------------------------------------------------------------------

  FUNCTION NORMESCAL(VarX,VarY)

    !&E---------------------------------------------------------------------
    !&E                 ***  FUNCTION NORMESCAL  ***
    !&E
    !&E ** Purpose : Compute the norm of a vector.
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by : 
    !&E
    !&E ** EXTERNAL calls : 
    !&E
    !&E ** Used ij-arrays :
    !&E
    !&E ** Modified variables : 
    !&E
    !&E ** Reference :
    !&E
    !&E ** History :
    !&E
    !&E---------------------------------------------------------------------
    !! * Modules used
    IMPLICIT NONE

    !! * Arguments
    REAL(kind=rsh) :: VarX, VarY, NORMESCAL

    !! * Local declarations

    !!----------------------------------------------------------------------
    !! * Executable part
    NORMESCAL = SQRT( VarX**2 + VarY**2 )

  END FUNCTION NORMESCAL

  !!----------------------------------------------------------------------

END MODULE waveforcing
