  SUBROUTINE step2dxy

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE step2dxy  ***
!&E
!&E ** Purpose : Estimate the dynamics and temperature, salinity
!&E              equations at n+1/2
!&E
!&E ** Description : dyn2dxy estimates xe^n+1/2 and v^n+1/2
!&E                  plus an intermediate variable u_i^n+1/2
!&E                  tssub2d_eq solves temperature, salinity
!&E                  and substance equations
!&E
!&E ** Called by : step.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : u,v,ssh,sal,temp,cv_wat
!&E
!&E ** Modified variables : uv,xeuv,uv2,u,v,sal,temp,cv_wat
!&E                         phitemp,phieau,phisal
!&E
!&E ** Reference :
!&E                T. Duhaut, M. Honnorat, L. Debreu, 2008 : Rapport
!&E                contrat PREVIMER - Ref : 06/2 210 290 (pages 100-125)
!&E
!&E ** History :
!&E       !           (P. Lazure)  Original code
!&E       !  2007-03  (F. Vandermeirsch) introduction of characteristics method
!&E       !  2008-07  (L. Debreu, T. Duhaut, M. Honnorat) time scheme
!&E       !  2008-07  (T. Odaka) MPI
!&E       !  2009-02  (A. Mano, V. Garnier, T. Odaka) MPI2D
!&E       !  2009-04  (V. Garnier) separation of obc origin (ogcm/tide)
!&E                                removing of xeuv argument for tide_harmcompo and obc2d_rank
!&E       !  2009-05  (F. Dumas)  AGRIF to new time scheme
!&E       !  2009-05  (F. Vandermeirsch)  introduction l_obc_tide,l_obc_ogcm
!&E       !  2009-07  (F. Vandermeirsch)  MC 2nd order scheme
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d,     ONLY : l_modele2d,ig,id,jb,jh,l_obc2drank_save, &
                             nbouc,l_champroche,icon,h0fond,l_tide_surge
   USE comvarp2d,     ONLY : u,v,ssh,sal,temp
   USE traobc2dapply, ONLY : ts_obc2d_apply
   USE obcapplyssh,   ONLY : obc_apply_ssh

   USE comsubstance,  ONLY : nbrej,nv_state,cv_wat
   USE subobc2dapply, ONLY : sub_obc2d_apply
   USE obcogcm,       ONLY : obc_subs
   USE obccombine,    ONLY : l_obc_ogcm_rc

   USE comfluxdo,     ONLY : xeadv,nbfl
   USE comaree,       ONLY : l_tide_harmcompo

   USE obcogcm,       ONLY : obc_ogcm

   USE obccombine,    ONLY : l_obc_ogcm_rssh,obc_combine,l_obc_tide,l_obc_ogcm,l_obc_char

   USE obcapplyu,     ONLY : obc_apply_u

   USE toolmpi,       ONLY : ex_i_rsh,ex_j_rsh


   IMPLICIT NONE

!! * Arguments

!! * Local declarations
   INTEGER                                        :: i,j
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:,:)    :: xeuv
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2) :: uv,uv2

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

! ----------------------------------------
! Initialisation of intermediate variables
! ----------------------------------------

! xeuv = estimated sea surface height
! uv = estimated barotropic current

   ALLOCATE(xeuv(liminm2:limaxp2,ljminm2:ljmaxp2))
   xeuv(liminm2:limaxp2,ljminm2:ljmaxp2)=ssh(liminm2:limaxp2,ljminm2:ljmaxp2)

   uv(liminm2:limaxp2,ljminm2:ljmaxp2)=0.0_rsh
   uv2(liminm2:limaxp2,ljminm2:ljmaxp2)=0.0_rsh


! --------------------------------------------
! change the time step if possible or required
! --------------------------------------------

   CALL dyn2d_dt(u,v)

! ----------------------------
! Read and apply open boundary
! ----------------------------

   !IF (Agrif_Root()) THEN
     IF (l_obc_tide) THEN
       IF (l_tide_surge) THEN
         CALL tide_harmcompo
         CALL obc2d_rank
       ELSE
         IF (l_tide_harmcompo) THEN
           CALL tide_harmcompo
         ELSE
           IF (l_testcase) THEN
             CALL obc2d_cas(xeuv)
           ELSE
             CALL obc2d_rank
           END IF
         END IF
       END IF
     END IF

     IF (l_obc_ogcm) CALL obc_ogcm


     IF (l_obc_ogcm_rc) CALL obc_subs

     CALL obc_combine(xeuv)         ! combine ogcm and tide obc
     IF (.NOT. l_obc_char) CALL obc_apply_ssh(xeuv)       ! apply ssh obc
   !ELSE
     !CALL agrif_obc_ssh(xeuv)
   !ENDIF


! ---------------------
! Estimate the dynamics
! ---------------------

   CALL dyn2dxy(u,v,ssh,uv,xeuv,uv2)

! Update u,v,ssh variables

   DO j=ljminp1,ljmax
     ssh(MAX0(limin,ig(j)-ishift_obc_west):MIN0(limax,id(j)+ishift_obc_east),j)= &
                                         xeuv(MAX0(limin,ig(j)-ishift_obc_west):MIN0(limax,id(j)+ishift_obc_east),j)
   END DO

   CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh(liminm2:limaxp2,ljminm2:ljmaxp2))

   DO j=ljminm2,ljmaxp2
     u(MAX0(liminm2,ig(j)-1-ishift_obc_west):MIN0(limaxp2,id(j)+ishift_obc_east),j)= &
                                             uv(MAX0(liminm2,ig(j)-1-ishift_obc_west):MIN0(limaxp2,id(j)+ishift_obc_east),j)
     v(MAX0(liminm2,ig(j)  -ishift_obc_west):MIN0(limaxp2,id(j)+ishift_obc_east),j)= &
                                            uv2(MAX0(liminm2,ig(j)  -ishift_obc_west):MIN0(limaxp2,id(j)+ishift_obc_east),j)
   END DO
   DEALLOCATE(xeuv)

! --------------------------------------------------------
! Save open boundary conditions (xe,u,v) for the next rank
! --------------------------------------------------------

   IF (l_obc2drank_save) CALL tide_saveobc(ssh,u,v)

! Initialize phicon and phieau
   CALL init_flux

! -------
! Outflow
! -------


   IF (nbrej > 0) CALL outflow


! -----------------------------------------------------------
! Estimate the temperature, salinity and substances equations
! -----------------------------------------------------------


   IF (l_modele2d .AND. (nv_state >= 1 .OR. nbfl >= 1)) THEN


     CALL sflx_surf(temp,sal,ssh)


     CALL tssub2d_eq(                                             &

                     cv_wat,                                      &

                     sal,temp,xeadv,ssh )

     !IF (Agrif_Root()) THEN

       CALL ts_obc2d_apply(u,v,sal,temp)


       CALL sub_obc2d_apply(u,v,cv_wat)

     !ELSE
       !CALL ts_obc2d_apply_agrif

       !CALL sub_obc2d_apply_agrif

     !ENDIF

     IF (icon == 2) CALL ts_stateeq2d(sal)



     DO j=ljminp2,ljmaxm1
       xeadv(MAX0(limin,ig(j)+1):MIN0(limax,id(j)-1),j)=ssh(MAX0(limin,ig(j)+1):MIN0(limax,id(j)-1),j)
     END DO

   ELSE IF (.NOT.l_modele2d) THEN

     DO j=ljminp2,ljmaxm1
       xeadv(MAX0(limin,ig(j)+1):MIN0(limax,id(j)-1),j)=ssh(MAX0(limin,ig(j)+1):MIN0(limax,id(j)-1),j)
     END DO

   END IF

! -----------------------
! Estimate the trajectory
! -----------------------



  END SUBROUTINE step2dxy
