!&E---------------------------------------------------------------------


  PROGRAM mars2d3d


!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE main  ***
!&E
!&E ** Purpose : main program
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls : ionc4_init, init, bathy, grid_mask,grid_metric
!&E                   : tool_decompdate, meteo, init_xe, init_uvz, init_uv, init_viscoef
!&E                   : init_turbcoef, traceur, river, outflow, init_tssub, ts_stateeq_linear
!&E                   : ts_stateeq_blumberg, sortie, diagnostic, track_points, step
!&E                   : Agrif_Step
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E        !  1987-02  (J.-C. Salomon, P. Lazure) original code
!&E        !  2007-07  (M. Sourisseau) elimination var_read_obc
!&E        !  2009-02  (V. Garnier) MPI2D
!&E        !  2009-02  (P. Le Hir, F. Cayocca) introduction "new" sand/mud module
!&E        !  2009-07  (N. Bruneau, F. Dumas) coupling with swan
!&E        !  2009-09  (B. Thouvenin) coupling with MET&OR (contaminant)
!&E        !  2009-10  (V. Garnier) reproductibility when using restart file
!&E        !  2010-03  (A.-C. Bennis, T. Odaka) use of palm coupling
!&E        !  2010-06  (V.M. Faure, N. Bousserez, M. Baklouti) coupling with eco3m
!&E        !  2010-06  (F. Cayocca) coupling with the module of sedimentology
!&E        !  2010-10  (B. Thouvenin) Distribution of the river discharge throughout a thickness
!&E        !  2011-01 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
!&E        !  2010-10  (B. Thouvenin) Distribution of the river discharge throughout a thickness
!&E        !  2011-03  (M. Huret) Introduction of Individual Based Model
!&E        !  2013-07 (A. Thevenin - CERFACS) coupling using OASIS3-MCT
!&E        !  2014-07  (F. Ganthy) Introduction of obstruction module
!&E        !  2016-12  (B.Thouvenin) version sedim_MUSTANG : only mixsed (mudsand not operationnal)
!&E---------------------------------------------------------------------





























































!!! definitions de dimensions, variables, parametres pour MUSTANG
!!! et couplage avec autre modele




!&E==========================================================================
!&E                   ***    coupleur_MUSTANG  ***
!&E==========================================================================

!/* Spatial MPI limits grid
!*/





!/* Spatial Grid limits definition  of loops
!*/







!/* dimensions table definition
!*/







!/* general variable hydro , bathy, time ... defined in hydro model but using by MUSTANG
!*/


























!/* number of advected and total substances
!*/
! nombre de variables transportees

! nombre de variables totales (variable transportees,+variables fixees, +variables forcantes etc..)


!/* name of variables in MUSTANG (dimensions order could changed)
!*/
! layers thickness (ordre des dimensions differentes)

! vertical eddy viscosity  (ordre des dimensions differentes)


!/* name of fluxes exchange between MUSTANG and hydro model
!*/











!/* dimensions of variables in hydro modele
!*/!! dimensions des variables 3D ou 4D dans le modele hydro !! (si pas MARS)

!!!#define DIM_MODELE_HYDRO_cw_bottom nv_tot,liminm1:limaxp2,ljminm1:ljmaxp2
!#define DIM_COORD_SIG 0:kmax+1





!/* Laterale Erosion
!*/
! Erosion laterale dans les mailles voisines (pourrait dependre architecture grille)













!/* sliding proces of fluid mud
!*/
! glissement depot : pente du fond dans les mailles voisines  (pourrait dependre architecture grille)









!! * Modules used
   USE parameters
   USE ionc4, ONLY : ionc4_init,ionc4_closeall
   USE comvars2d
   USE comvars3d
   USE comvarp2d
   USE comvarp3d
   USE comcolig2d,    ONLY : igdv_old
   USE comcolig3d,    ONLY : intz_exp_vstar
   USE comdiag,       ONLY : l_diag,file_diag
   USE output,        ONLY : output_init,output_mng,l_out_demer,fileresu,nb_out
   USE comdem
   USE comfluxdo
   USE comturb,       ONLY : eps
   USE obcogcm,       ONLY : obc_tsref


!   USE wave_tools,    ONLY : wave_stokes_fill_tabkm


   USE subreaddat,    ONLY : sub_read_var,sub_read_vardiag
   USE subreadflxatm, ONLY : sub_sflx_read
   USE comsubstance,  ONLY : nv_dri, nv_state, cv_wat,l_cvrain_readfile,l_subflxatm_readfile,l_bilan
   USE usersubstance, ONLY : user_init_vdriv
   USE obcogcm,       ONLY : obc_subsref
   USE submassbalance,ONLY : sub_budget_readdomain





   USE comMUSTANG,        ONLY : MUSTANG_alloc,cv_sed,l_repsed,l_morphocoupl,coef_erolat
   USE initMUSTANG,       ONLY : MUSTANG_compatibility,MUSTANG_param,MUSTANG_sedinit
   USE sed_MUSTANG_MARS,  ONLY : sedinit_fromfile
   USE sedimorpho,        ONLY : sed_morphoinit,sed_alloc_morpho
   USE coupleur_MUSTANG,  ONLY : coupl_conv2MUSTANG

   USE sflxatm,     ONLY : drag,exch,imeteo_dragtype,imeteo_exchtype
   USE sflxturb,    ONLY : sflx_coef_init
   USE trackpoints, ONLY : l_points,track_points






   USE casinit,     ONLY : cas_init
   USE toolmpi, ONLY : MPI_START1,MPI_START2,MPI_START3,MPI_END1,NPROCS
   USE toolmpi, ONLY : ADD_TAG,ex_j_rsh,ex_i_rsh,MASTER

   USE tooltiming,  ONLY : set_iluout_timing, time_stat_ll, time_header_ll, second_mars, &
                           timing_legend, timing_separator, xt_start, xt_diag, &
                           xt_store, xt_bound, xt_forcing, xt_nudging, xt_traj, &
                           xt_dynamics, xt_turb, xt_tsbuoy, xt_subs, xt_points, &
                           xt_misc, xt_oasis, xt_agrif
   USE toolmpi, ONLY : xt_mpi_com,xt_mpi_wait


!! * Arguments

!! * Local declarations

   IMPLICIT NONE



   INTEGER                            :: taille_germe,narg,iargc,nb
   INTEGER, ALLOCATABLE, DIMENSION(:) :: germe
   INTEGER                            :: iluout      ! Logical unit number for the output listing


   CHARACTER(LEN=3)     :: tag

! for computing time analysis
   REAL(KIND=rlg),DIMENSION(2)        :: ztime0,ztime1,ztime2,ztot,zall,xt_init,xt_step

!$   INTEGER  :: OMP_GET_NUM_THREADS

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



   ! Agrif_Parallel_Sisters = .FALSE.

   CALL MPI_START1

! 0.0 Time logger initialization
! ------------------------------
   CALL second_mars(ztime0)
   xt_start     = 0.0_rlg ; xt_store     = 0.0_rlg ; xt_bound     = 0.0_rlg ;
   xt_forcing   = 0.0_rlg ; xt_nudging   = 0.0_rlg ; xt_dynamics  = 0.0_rlg ;
   xt_turb      = 0.0_rlg ; xt_tsbuoy    = 0.0_rlg ; xt_subs      = 0.0_rlg ;
   xt_diag      = 0.0_rlg ; xt_traj      = 0.0_rlg ; xt_points    = 0.0_rlg ;
   xt_misc      = 0.0_rlg ; xt_oasis     = 0.0_rlg ; xt_agrif     = 0.0_rlg ;

   xt_mpi_wait = 0.0_rlg
   xt_mpi_com = 0.0_rlg





! 1.0 Initialization
! ------------------

!  Init des boucles
!------------------------
   nbouc=0
   nbouc_macho=0

! initialisation de la procedure de tirage des nombres aleatoires
!-----------------------------------------------------------------
   root='./'
   CALL RANDOM_SEED(size=taille_germe)
   ALLOCATE(germe(taille_germe))
   germe(:)=3
   CALL RANDOM_SEED(put=germe(:))


   CALL ionc4_init()

  
   narg=iargc()

   IF (l_testcase) THEN

! Initialize test case configuration
     CALL cas_init

! Initialize coefficients of turbulence
     IF (turb_nbeq==1 .OR. turb_nbeq==2) CALL init_turbcoef(ect,psi)

! Initialize xeadv
     xeadv(limin:limax,ljmin:ljmax)=ssh(limin:limax,ljmin:ljmax)

! Restart from a file ( restart file from MARS or climatological file )
     IF (l_initfromfile) THEN
       CALL init_fromfile(u,uz,v,vz,ssh,sal,temp,        &
                               intz_exp_vstar,igdv_old,  &

                                            cv_wat,      &

                                            ect,psi,eps)
     END IF

   ELSE

! Read namelists (paramain.txt, paracom.txt, paraspec.txt)
! Read head.conf
     CALL init(narg)

! Read parameters of XXXXX


! Initialize flux exchange coefficients (to prepare the bilinear interpolation)
     IF (imeteo_exchtype==1) CALL sflx_coef_init(drag,exch,imeteo_dragtype)

! For MARS_WEB
     IF (narg>=2) CALL init_marsweb_clim
     IF (narg==3) CALL init_marsweb_simu

! Read the bathymetry
     CALL bathy

! Initialize ig,id,jb,jh
     CALL grid_mask

! Divide the domain for MPI run
     CALL MPI_START2(1,1)
     CALL MPI_START3
! Toujours sur grille mere uniquement
     !CALL Agrif_Declare_Variables


! Allocation of variables
! -----------------------
! Dynamical variables
     CALL alloc
! Intermediate variables
     CALL alloc_inter
! Stationnary variables (grid metrics)
     CALL alloc_stat

! Initialize metrics (dx,dy,sig,lon,lat...)
     CALL grid_metric



! Read file variables.dat describing substances
! Allocate substance variables
! Make variables in order
     CALL sub_read_var

! Initialize open boundary arrays for substance
     CALL obc_subsref





     CALL MUSTANG_alloc




! Read file vardiag.dat describing diagnoses
! Allocate diagnoses variables




! Initialize sea surface height
     CALL init_xe(ssh)
! Initialize height variables for vertical integration
     CALL init_height(ssh)
! Initialize 3D velocities
     CALL init_uvz(uz,vz)
! Initialize 2D velocities
     CALL init_uv(u,v,uz,vz)
! Initialize viscosity coefficients
     CALL init_viscoef
! Initialize coefficients of turbulence
     CALL init_turbcoef(ect,psi)
! Initialize sal, temp, cv_wat
     CALL init_tssub(sal,temp   &

                        ,cv_wat &

                             )
     CALL obc_tsref
! Initialize obstructions


! Restart from a file ( restart file from MARS or climatological file )
     IF (l_initfromfile) THEN
       CALL init_fromfile(u,uz,v,vz,ssh,sal,temp,        &
                               intz_exp_vstar,igdv_old,  &

                                            cv_wat,      &

                                            ect,psi,eps)
     END IF



! coupleur pour recuperation des concentrations au fond
! pour iappel=0 seuls cv_wat,sal,temp vont etre utilises
    CALL coupl_conv2MUSTANG(limin,limax,ljmin,ljmaxm1,0,  &
                            h0,ssh,   &

                       cv_wat,sal,temp )


     IF (l_repsed) THEN

       CALL sedinit_fromfile(limin,limax,ljmin,ljmaxm1,  &
                                 rhoref,h0)

     ELSE
       CALL MUSTANG_sedinit(limin,limax,ljmin,ljmaxm1,  &
                         rhoref,h0)
     END IF
!CALL sed_settlveloc_ini
     IF(l_morphocoupl) THEN
        CALL sed_alloc_morpho
        IF(.NOT.l_initfromfile)CALL sed_morphoinit
     ENDIF

! Initialize wave fields


! Initialize ecosystem variables


!!!! CAR INITIALISE DANS STEP POUR SEDIM REVOIR ?#if defined key_sedim_MUSTANG_anc || defined key_contaminant




! initialization of driving variables
     IF (nv_dri > 0) &
       CALL user_init_vdriv(cv_wat                       &

                           ,cv_sed                       &

                                )


! Read output.dat and allocate variables with the correct number of output files
     CALL output_init

   END IF   ! ends test on l_testcase








!!! introduction of submassbalance
     IF (l_bilan) CALL sub_budget_readdomain(ssh)
     !PRINT*, 'call read_limits'
!!! end submassbalance


   CALL ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,v(liminm2:limaxp2,ljminm2:ljmaxp2))   ! for dyn2d_dt val
   CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2))   ! for dyn2d_dt val

! Save summary of simulation
   IF (MASTER) THEN
   CALL summary

   IF(.NOT. l_testcase)CALL MUSTANG_compatibility

   ENDIF

! Initialize listing file
   IF (narg == 0 .AND. .NOT. l_testcase) THEN
     filename= 'listing'
     CALL ADD_TAG(filename,LEN(filename))
     CALL ADD_TAG(file_diag,LEN(file_diag))
!     CALL edge_initialise_MPI
     IF (MASTER) THEN
     IF (iscreen/=6) OPEN(unit=iscreen,file=filename,status="replace")
     ENDIF
   END IF

! Initialize meteorological variables if l_meteo_stat or l_meteo_hom
   IF (.NOT. l_testcase) CALL meteo

! Read file river.dat, initialize river discharges inputs with time interpolation
   CALL river(ssh)


! Read file outflow.dat, file describing discharges
! initialize discharges inputs with time interpolation
   CALL outflow

! Read files sflx_cvrain and/or sflx_subflxatm
! if temporal variations of cvrain or of atm. fluxes for at least one substance
!  (l_cvrain_readfile=.true. or l_subflxatm_readfile=.true.)
   IF (l_cvrain_readfile .or. l_subflxatm_readfile) CALL sub_sflx_read



! Estimate buoyancy from sal and temp
   IF (l_modele2d) THEN
     IF (icon == 2) CALL ts_stateeq2d(sal)
   ELSE
     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)
   END IF

   CALL second_mars(ztime2)
   xt_start = xt_start + ztime2 - ztime0
   ztime1 = ztime2

! Intialize variables for traking
   IF (l_points) CALL track_points(u,v,uz,vz,ssh,sal,temp, &

                                   cv_wat,                 &

                                   ect,boz)

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

! Initialize diagnoses variables
   IF (l_diag) CALL diagnostic(u,v,uz,vz,ssh,sal,temp, &

                               cv_wat,                 &

                               ect,boz)




!   CALL wave_stokes_fill_tabkm


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

! manage outputs files
   CALL output_mng

! Initialize demerliac filter
   DO nb=1,nb_out

     IF (l_out_demer(nb)) &
        CALL out_demerliac(ssh(limin:limax,ljmin:ljmax),                                     &
                           u(limin:limax,ljmin:ljmax),uz(:,limin:limax,ljmin:ljmax),         &
                           v(limin:limax,ljmin:ljmax),vz(:,limin:limax,ljmin:ljmax),         &
                           sal(:,limin:limax,ljmin:ljmax),temp(:,limin:limax,ljmin:ljmax),   &
                           xed(limin:limax,ljmin:ljmax),                                     &
                           ud(limin:limax,ljmin:ljmax),uzd(:,limin:limax,ljmin:ljmax),       &
                           vd(limin:limax,ljmin:ljmax),vzd(:,limin:limax,ljmin:ljmax),       &
                           sald(:,limin:limax,ljmin:ljmax),tempd(:,limin:limax,ljmin:ljmax), &

                           cv_wat(:,:,limin:limax,ljmin:ljmax),cv_watd(:,:,limin:limax,ljmin:ljmax), &

                           nsdem)
   END DO

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

! Initialize diagnoses variables relative to terms of the momentum equation
! after output_mng to be allowed to complete the main output file ( fileresu(1) )









   xt_init = xt_start + xt_points + xt_diag + xt_store + xt_traj


! Save the repartition of the cpus
! --------------------------------


   IF (MASTER) THEN
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) ' NUMBER OF MPI CPUS' , nprocs
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) ' '
   ENDIF

   iscreenlog=iscreenlog   ! because of bug in AGRIF conv
   IF (MASTER) THEN
!$OMP PARALLEL
!$OMP MASTER
!$   WRITE(iscreenlog,*) ' '
!$   WRITE(iscreenlog,*) '**************************************************'
!$   WRITE(iscreenlog,*) ' NUMBER OF OMP THREADS' , OMP_GET_NUM_THREADS()
!$   WRITE(iscreenlog,*) '**************************************************'
!$   WRITE(iscreenlog,*) ' '
!$OMP END MASTER
!$OMP END PARALLEL
   ENDIF

! 2.0 Time integration
! --------------------

   IF (MASTER) THEN
   WRITE(iscreen,*) ' '
   WRITE(iscreen,*) '**************************************************'
   WRITE(iscreen,*) '**************** Start time integration **********'
   WRITE(iscreen,*) '**************************************************'
   WRITE(iscreen,*) ' '
   ENDIF

   halfdt = dt*0.5_rlg

! Time integration

   DO WHILE (t <= tfin .OR. MOD(nbouc,2)==0)



     CALL step()
   END DO


   CALL second_mars(xt_step)
!xt_step = xt_step - (xt_start + xt_diag + xt_store + xt_traj + xt_points)
   xt_step = xt_step - ztime1


! 3.0 Print statistics
! --------------------

! Set File Timing OUTPUT
   iluout = 1970
   CALL set_iluout_timing(iluout)

! Compute global time
   CALL time_stat_ll(xt_init+xt_step,ztot)  ! reference 100%

   CALL time_header_ll

   CALL time_stat_ll(xt_init,ztot,      'MAIN INITIALIZATION','=')
   CALL time_stat_ll(xt_step,ztot,       'MAIN TIME_STEP','=')
   CALL timing_separator('=')
   CALL timing_separator(' ')

   CALL timing_legend()
   CALL time_stat_ll(xt_start,ztot,      ' INITIALIZATION','=')
   CALL time_stat_ll(xt_bound,ztot,      ' BOUNDARIES','=')
   CALL time_stat_ll(xt_store,ztot,      ' STORE-FIELDS','=')
   CALL time_stat_ll(xt_dynamics,ztot,   ' DYNAMICS (mtm cty)','=')
   CALL time_stat_ll(xt_forcing,ztot,    ' FORCING','=')
   CALL time_stat_ll(xt_tsbuoy,ztot,     ' T_S_BUOY','=')
   CALL time_stat_ll(xt_turb,ztot,       ' TURBULENCE','=')
   CALL time_stat_ll(xt_nudging,ztot,    ' NUDGING','=')
   CALL time_stat_ll(xt_points,ztot,     ' tracking POINTS','=')
   CALL time_stat_ll(xt_diag,ztot,       ' DIAG','=')
   CALL time_stat_ll(xt_traj,ztot,       ' TRAJ','=')
   CALL time_stat_ll(xt_misc,ztot,       ' MISC','=')
   CALL time_stat_ll(xt_oasis,ztot,      ' OASIS','=')
   CALL time_stat_ll(xt_agrif,ztot,      ' AGRIF','=')
   CALL time_stat_ll(xt_subs,ztot,       ' SUBSTANCES','=')

   CALL time_stat_ll(xt_mpi_wait,ztot,   ' MPI_WAIT','=')
   CALL time_stat_ll(xt_mpi_com,ztot,    ' MPI_COM','=')

   CALL timing_separator('=')
   CALL timing_separator(' ')

! sum of call subroutine
   zall = xt_start + xt_bound  + xt_store  + xt_dynamics +  xt_forcing + &
          xt_tsbuoy+ xt_turb + xt_nudging + xt_points + xt_diag + xt_traj + &
          xt_misc + xt_oasis + xt_agrif + xt_subs
  
   CALL timing_legend()
   CALL time_stat_ll(zall,ztot,          ' SUM(CALL)','=')
   CALL second_mars(zall)
   CALL time_stat_ll(zall-ztime0,ztot,   ' TOTAL_RUN','=')
   CALL timing_separator('=')

! 4.0 Deallocations
! -----------------

   CALL dealloc_inter
   CALL dealloc_stat
   CALL ionc4_closeall()

   PRINT*, 'fin normale'

   CALL dealloc





   CALL MPI_END1




  END PROGRAM mars2d3d 

