  SUBROUTINE step3dyx

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE step3dyx  ***
!&E
!&E ** Purpose : estimate of the dynamics and heat/salt/substance equation
!&E              from time n+1/2 to n+1
!&E
!&E ** Description : estimate of
!&E                   - the final zonal velocities : U^n+1, UZ^n+1
!&E                   - the intermediate meridional velocities : V^n+1, VZ^n+1
!&E                   - the final sea surface height : XE^n+1
!&E                   - temperature, salinity : temp^n+1, sal^n+1
!&E                   - substances : cv_wat^n+1
!&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---------------------------------------------------------------------





























































!! * 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

   USE obcogcm,        ONLY : obc_ogcm

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


   USE obcapplyv
   USE obcapplyvz
   USE traobc3dapply
   USE obcapplyssh

   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





   USE toolmpi,        ONLY : ex_i_rsh,ex_j_rsh

   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*, 'enter STEP3DYX 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)

   ! 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



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

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

   !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) THEN

       CALL obc_ogcm ! read ogcm open boundary conditions
     END IF


     IF (l_obc_ogcm_rc) CALL obc_subs  ! read "bio" open boundary conditions

     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 of the dynamics
! ------------------------

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

   !IF (.NOT.agrif_grandmother()) THEN
   CALL dyn3dyx(v,u,ssh,uv,uv2,xeuv,vz,uz,boz)
   !ENDIF



! Estimate of the vertical velocity

!toseelater We can remove uz,vz,u,v in dyn3dwz
   !IF (.NOT.agrif_grandmother()) THEN

   CALL dyn3d_wz(ssh,xeuv,uz,vz)

   !ENDIF

! Save intermediate barotropic zonal velocity and the sea surface height
! into corresponding global 2d variables
! inside the interior domain (obc excluded)

   DO j=ljmin,ljmax
     DO i=liminp1,limax
       IF ((j>=jb(i)-jshift_obc_south) .AND. (j<=jh(i)+jshift_obc_north)) &
         ssh(i,j)=xeuv(i,j)
     END DO
   END DO

   DEALLOCATE(xeuv)

   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
     DO i=liminm2,limaxp2
       IF ((j.GE.jb(i)-1-jshift_obc_south) .AND. (j.LE.jh(i)+jshift_obc_north)) v(i,j)=uv(i,j)
       IF ((j.GE.jb(i)  -jshift_obc_south) .AND. (j.LE.jh(i)+jshift_obc_north)) u(i,j)=uv2(i,j)
     END DO
   END DO

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

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

! Initialize phicon and phieau
   CALL init_flux









! Outflow

   IF (nbrej > 0) CALL outflow



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



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

! Estimate of atmospherical surface flux
! Update phitemp, phisal, phieau


   CALL sflx_surf(temp,sal,ssh)


! Solve equation of heat, salt and substances

   IF (icon >= 1) THEN

     CALL tssub_eq(xeadv,ssh,sal,temp  &

                  ,cv_wat              &

                  )



! 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



! Apply open boundary conditions on temperature, salinity and substances

     !IF (Agrif_Root()) THEN


       CALL ts_obc3d_apply(uz,vz,sal,temp)



       CALL sub_obc3d_apply(uz,vz,cv_wat)


     !ELSE

       !CALL ts_obc3d_apply_agrif

       !CALL sub_obc3d_apply_agrif

     !ENDIF

! Estimate the buoyancy


     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)



     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 ! ends test on icon

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



! -------------------------------------------
! 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)

     CALL turbpsi(ssh,uz,vz,pn2,psi,ect)

     CALL turblepskznz(ssh,pn2,psi,ect)
   END IF


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

! key_offline


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



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




   !CALL agrif_update_fields
   CALL second_mars(ztime2)
   xt_agrif = xt_agrif + ztime2 - ztime1
   ztime1 = ztime2


   !PRINT*, 'exit STEP3DYX routine'

  END SUBROUTINE step3dyx
