  SUBROUTINE step3dxy

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE step3dxy  ***
   !&E
   !&E ** Purpose : estimate of the dynamics and heat/salt/substance equation
   !&E              from time n to n+1/2
   !&E
   !&E ** Description :  estimate of
   !&E                   - the intermediate zonal velocities : U^n+1/2, UZ^n+1/2
   !&E                   - the final meridional velocities : V^n+1/2, VZ^n+1/2
   !&E                   - the final sea surface height : XE^n+1/2
   !&E                   - temperature, salinity : temp^n+1/2, sal^n+1/2
   !&E                   - substances : cv_wat^n+1/2
   !&E
   !&E ** Called by : step
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays :
   !&E
   !&E ** Modified variables :
   !&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-07  (J. Marin et F. Vandermeirsch) introduction of characteristics method
   !&E       !  2008-07  (L. Debreu, M. Honnorat, T. Duhaut) time scheme
   !&E       !  2008-08  (T. Odaka) MPI
   !&E       !  2009-04  (V. Garnier) separation of obc origin (ogcm/tide)
   !&E                                 removing of xeuv argument for tide_harmcompo and obc2d_rank
   !&E       !  2009-02  (P. Le Hir, F. Cayocca) removing of sedimentology calls
   !&E       !  2009-04  (A. Mano, V. Garnier) MPI2D
   !&E       !  2009-05  (F. Dumas)  AGRIF to new time scheme
   !&E       !  2009-05  (F. Vandermeirsch)  introduction l_obc_tide,l_obc_ogcm,characteristics method
   !&E       !  2009-07  (F. Vandermeirsch)  MC 2nd order scheme
   !&E       !  2009-09  (B. Thouvenin)  coupling with MET&OR (contaminant)
   !&E       !  2011-03  (M. Huret) Introduction of Individual Based Model
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2015-04  (B.Thouvenin) call forcingKtest for some test cases just after initflux
   !&E       !  2017-09  (F. Ganthy) Turbulence advection (k-epsilon)
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d
   USE comvars3d
   USE comvarp2d
   USE comvarp3d ,     ONLY : uz,vz,boz,pn2,ect,psi
   USE cometdyn
   USE comfluxdo
   USE comaree
#ifdef key_lagoon
   USE obcogcm,        ONLY : obc_ogcm_lagoon
#else
   USE obcogcm,        ONLY : obc_ogcm
#endif
   USE obccombine,     ONLY : obc_combine,l_obc_tide,l_obc_ogcm,l_obc_char,l_obc_ogcm_ruvz
#ifdef key_trajec3d
   USE_MPI comtraj,    ONLY : down_give,up_give,right_give,left_give
   USE traject3d
#endif
#ifdef key_ibm
   USE ibm, ONLY : ibm_3d
#endif
   USE obcapplyu
   USE obcapplyuz
   USE traobc3dapply
   USE obcapplyssh
#if defined key_substance
   USE subobc3dapply,  ONLY : sub_obc3d_apply
   USE obcogcm,        ONLY : obc_subs
   USE obccombine,     ONLY : l_obc_ogcm_rc
   USE comsubstance,   ONLY : tdebsubs,nbrej,cv_wat,itypv,iv_ws_max, &
                              sv_part,ndt_part,nb_ws_max,nvp,l_bilan,nv_adv
   USE submassbalance, ONLY : sub_budget_flxcum
#if defined key_biolo
   USE comsubstance,   ONLY : dtbio
   USE biolodynzwat,   ONLY : biolo_dyn_zwat
   USE biolosinksource,ONLY : biolo_sksc_wat
#endif
#if defined key_contaminant
   USE contareaddat,   ONLY : l_flim3_re,nreacmax,l_treat_re_dc
   USE contasinksource,ONLY : conta_sksc_wat
   USE contaradsigz,   ONLY : radz_reduce
#endif
#endif
#ifdef key_offline
   USE offline
#endif
#ifdef key_MPI_2D
   USE toolmpi,        ONLY : ex_i_rsh,ex_j_rsh
#endif
   USE tooltiming,     ONLY : second_mars,xt_subs,xt_tsbuoy,xt_agrif,xt_dynamics,xt_traj,xt_turb

   IMPLICIT NONE

   !! * Arguments

   !! * Local declarations
   INTEGER            :: i,j,l,k,iv
   ! intermediate 2d variables of velocity and sea level
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:,:)      :: xeuv
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2)   :: uv,uv2
                                ! for computing time analysis
   REAL(KIND=rlg),DIMENSION(2)         :: ztime1,ztime2

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

   PRINT_DBG*, 'enter STEP3DXY routine'

   CALL second_mars(ztime1)

   !-----------------------------------------------------------------------
   ! fisrt step (half time step incrementation)
   !-----------------------------------------------------------------------

   ! Allocate and initialize 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)

   INST_AGRIF ssh_old(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

#if ! defined key_offline

   ! ------------------------
   ! Open Boundary Conditions
   ! ------------------------

#ifdef key_obc_char_barocline

   ! 3D nesting - characteristics method

!jm Calcul des modes verticaux ...
!! codage de la frontiere par un entier
!! Nord=1,Sud=2,Est=3,Ouest=4;

#ifdef key_obc_char_north
   CALL cal_modes(1,nModes,imax)
#endif

#ifdef key_obc_char_south
   CALL cal_modes(2,nModes,imax)
#endif

#ifdef key_obc_char_east
   CALL cal_modes(3,nModes,jmax)
#endif

#ifdef key_obc_char_west
   CALL cal_modes(4,nModes,jmax)
#endif

#endif

   ! offline : sea surface height is prescribed
   ! agrif : nesting of u and xe

   IF_AGRIF (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

#ifdef key_lagoon
     CALL obc_ogcm_lagoon
#else
     IF (l_obc_ogcm) THEN
#if defined key_siggen || defined key_gencoord
       IF (l_obc_ogcm_ruvz) CALL obc_init_height(xeuv)
#endif
       CALL obc_ogcm ! read ogcm open boundary conditions
     END IF
#endif
#ifdef key_substance
     IF (l_obc_ogcm_rc) CALL obc_subs  ! read "bio" open boundary conditions
#endif
     CALL obc_combine(xeuv)         ! combine ogcm and tide obc
     IF (.NOT. l_obc_char) CALL obc_apply_ssh(xeuv)       ! apply ssh obc
   ELSE_AGRIF
     CALL_AGRIF agrif_obc_ssh(xeuv)
   ENDIF_AGRIF

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

   ! intermediate estimate of the zonal velocity : U^n+1/2 and UZ^n+1/2
   ! final estimate of the meridional velocities : V^n+1/2 and VZ^n+1/2
   ! final estimate of the sea surface height : XE^n+1/2

   IF_AGRIF (.NOT.agrif_grandmother()) THEN
   CALL dyn3dxy(u,v,ssh,uv,uv2,xeuv,uz,vz,boz)
   ENDIF_AGRIF

#elif key_offline
   CALL dyn3dxy_offline(ssh,xeuv,uz,vz,boz,pn2)
   ! key_offline
#endif

   ! Estimate of the vertical velocity

!toseelater We can remove uz,vz,u,v in dyn3dwz
   IF_AGRIF (.NOT.agrif_grandmother()) THEN
#if defined key_siggen || defined key_gencoord
!ATTENTION modif timescheme
   CALL dyn3d_wz(ssh,xeuv,uz,vz,uv,v)
#else
   CALL dyn3d_wz(ssh,xeuv,uz,vz)
#endif
   ENDIF_AGRIF

   ! Save intermediate barotropic zonal velocity and the sea surface height
   ! into corresponding global 2d 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

   DEALLOCATE(xeuv)

   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI 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
   ! If rank nesting, save 2d variables at open boundaries of the child domain
   IF (l_obc2drank_save) CALL tide_saveobc(ssh,u,v)

   CALL second_mars(ztime2)
   xt_dynamics = xt_dynamics + ztime2 - ztime1
   ztime1 = ztime2

   !----------
   ! Ecosystem
   !----------

   ! Initialize phicon and phieau
   CALL init_flux

#if defined key_castest
   IF (num_testcase==5 .OR. num_testcase == 36 .OR. num_testcase == 37 .OR. num_testcase == 40) THEN
       CALL forcing_Ktest
   ENDIF
#endif

#if defined key_substance && defined key_biolo
   IF (t .GE. tdebsubs) THEN
     dtbio=0.5_rlg*dt
     CALL biolo_dyn_zwat(ssh,cv_wat,sal,temp)
   END IF
#endif

#if defined key_substance && defined key_contaminant 
!  attenuation of radiation into water column
   IF (t .GE. tdebsubs .and.(l_flim3_re) ) THEN
     CALL radz_reduce(ssh)
   END IF
#endif

#ifdef key_substance

   ! Outflow

   IF (nbrej > 0) CALL outflow

#ifdef key_biolo
   IF (t .GE. tdebsubs) THEN
     CALL biolo_sksc_wat(cv_wat,sal,temp,ssh)
   END IF
#elif key_contaminant
   IF (t .GE. tdebsubs) THEN
     IF (nreacmax > 0 .AND. l_treat_re_dc) THEN
        CALL conta_sksc_wat( cv_wat,sal,temp,ssh)
     END IF
   END IF
#endif

   CALL second_mars(ztime2)
   xt_subs = xt_subs + ztime2 - ztime1
   ztime1 = ztime2

#endif

   ! -------------------------------------
   ! Equation of heat, salt and substances
   ! -------------------------------------

   ! Estimate of atmospherical surface flux
   ! Update phitemp, phisal, phieau, phicon, (phicon_atm)

#if ! defined key_offline
   CALL sflx_surf(temp,sal,ssh)
#endif

   ! Solve equation of heat, salt and substances

   IF (icon >= 1) THEN

     CALL tssub_eq(xeadv,ssh,sal,temp  &
#ifdef key_substance
                  ,cv_wat              &
#endif
                  )

#ifdef key_substance

   ! Solve vertical equation of variables with settling velocity
     IF (nvp>0 .AND. t >= tdebsubs) THEN
       DO l=1,nb_ws_max
         ! not for gravel
         IF(itypv(iv_ws_max(l,1)) > 1 ) THEN
              CALL parsub_eq(xeadv,ssh,cv_wat,sv_part(l),ndt_part(l),l)
         ENDIF
       END DO
     END IF
     
     ! cumulating substances fluxes for budget estimation
     IF (l_bilan .AND. t >= tdebsubs) THEN
       CALL sub_budget_flxcum(cv_wat,ssh)
     END IF


#endif

   ! Apply open boundary conditions on temperature, salinity and substances

     IF_AGRIF (Agrif_Root()) THEN

#if ! defined key_offline
       CALL ts_obc3d_apply(uz,vz,sal,temp)
#endif

#ifdef key_substance
       CALL sub_obc3d_apply(uz,vz,cv_wat)
#endif

     ELSE_AGRIF
#ifdef key_agrif_obc_sponge_ts
!       CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal(:,liminm1:limaxp2,ljminm1:ljmaxp2))
!       CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal(:,liminm1:limaxp2,ljminm1:ljmaxp2))
!       CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
!       CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
       CALL_MPI ex_i_rsh(-1,1,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal(:,liminm1:limaxp2,ljminm1:ljmaxp2))
       CALL_MPI ex_j_rsh(-1,1,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal(:,liminm1:limaxp2,ljminm1:ljmaxp2))
       CALL_MPI ex_i_rsh(-1,1,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
       CALL_MPI ex_j_rsh(-1,1,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
       CALL_AGRIF ts_obc3d_apply_agrif
#ifdef key_substance
       CALL_AGRIF sub_obc3d_apply_agrif
#endif
     ENDIF_AGRIF

   ! Estimate the buoyancy

#if ! defined key_offline
     IF (icon >= 1 .AND.      l_stateeq_lin) CALL ts_stateeq_linear(sal,temp,boz,pn2)
     IF (icon >= 1 .AND. .NOT.l_stateeq_lin) CALL ts_stateeq_blumberg(sal,temp,boz,pn2)
#endif

#if ! defined key_turbadv
     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
#endif

   END IF ! ends test on icon

   CALL second_mars(ztime2)
   xt_tsbuoy = xt_tsbuoy + ztime2 - ztime1
   ztime1 = ztime2

#if ! defined key_offline

   ! -------------------------------------------
   ! Estimate of vertical turbulent coefficients
   ! -------------------------------------------

   IF (turb_nbeq == 0) &
     CALL turb(ssh,uz,vz,pn2)

   IF (turb_nbeq == 1) &
     CALL turbkl(ssh,uz,vz,boz,pn2,ect)

   IF (turb_nbeq == 2) THEN
     CALL turbk(ssh,uz,vz,pn2,ect)
#ifdef key_turbadv
     CALL turbpsi(xeadv,ssh,uz,vz,pn2,psi,ect)
#else
     CALL turbpsi(ssh,uz,vz,pn2,psi,ect)
#endif
     CALL turblepskznz(ssh,pn2,psi,ect)
   END IF

#ifdef key_turbadv
   IF (icon >= 1) 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)
     ENDDO
   ENDIF
#endif

   CALL second_mars(ztime2)
   xt_turb = xt_turb + ztime2 - ztime1
   ztime1 = ztime2

   ! key_offline
#endif

   ! ------------
   ! Trajectories
   ! ------------

#ifdef key_trajec3d
   CALL traj_3d(ssh,uz,vz)
   CALL second_mars(ztime2)
   xt_traj = xt_traj + ztime2 - ztime1
   ztime1 = ztime2
#endif

   ! ------------
   ! IBM
   ! ------------

#ifdef key_ibm
   CALL ibm_3d(ssh,uz,vz,sal,temp,          &
#ifdef key_substance
                                   cv_wat,  &
#endif
                                   boz       )
   CALL second_mars(ztime2)
   xt_traj = xt_traj + ztime2 - ztime1
   ztime1 = ztime2
#endif

#if ! defined key_agrif_1way
   CALL_AGRIF agrif_update_fields
   CALL second_mars(ztime2)
   xt_agrif = xt_agrif + ztime2 - ztime1
   ztime1 = ztime2
#endif

   PRINT_DBG*, 'exit STEP3DXY routine'

  END SUBROUTINE step3dxy
