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

  SUBROUTINE init(narg)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE init  ***
!&E
!&E ** Purpose : Read namelists and initialize most of the parameters
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !          Original code
!&E       !  2008-01 (M. Sourisseau) introduction MES (module de biologie)
!&E       !  2008-09 (P. Karleskind, V. Garnier) simultaneous saving of obc for different child grids
!&E       !  2008-10 (V. Garnier) output managing relative to displacement of init_fromfile from step to main
!&E       !  2008-10 (M. Jouan, V. Garnier) save conservative (flux meaning) variables uz and vz (for offline purpose)
!&E       !  2009-04 (V. Garnier) separation of obc origin (ogcm/tide) and
!&E                                read l_obc_* to choose different obc schemes
!&E       !  2009-04 (P. Le Hir, F. Cayocca) no reading of the sedimentology parametrization
!&E       !  2009-09 (B. Thouvenin) reading of contaminant parametrization
!&E       !  2009-10 (F. Dumas) Introduction of overtide functionality
!&E       !  2009-10 (V. Garnier) reproductibility when using restart file
!&E       !  2009-11 (S. Michel, HOCER) write errors and warnings in "error.log" and "warn.log"
!&E       !  2010-07 (R. Verney) Introduction of orbital velocity and wave period variables
!&E       !  2010-07 (B. Thouvenin) Modification for budgets of any substancee
!&E       !  2011-01 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
!&E       !  2014    (B.Thouvenin) Modification for key_biolo
!&E       !  2014-05 (R. Verney) introduction of module flocmod
!&E       !  2015-09 (F. Ganthy) Introduction of a spatialized bottom roughness length
!&E       !  2016-12 (A. Fiandrino, B.Thouvenin) Introduction of key_lagoon
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d
   USE comvars3d
   USE cometdyn
   USE comoasis
   USE comdiag
   USE cometemp
   USE comfluxdo

   USE comsubstance

   USE comMUSTANG
   USE initMUSTANG, ONLY : MUSTANG_param





   USE comaree
   USE comtraj,    ONLY : file_trajec,ndtz,itypetraj
   USE obcogcm,    ONLY : l_obc_ogcm_relax ,obc_coefrel
   USE obccombine, ONLY : l_obc_south, l_obc_north,l_obc_east,l_obc_west,     &
                          l_obc_ogcm_rssh,l_obc_ogcm_rt,l_obc_ogcm_rs,        &
                          l_obc_ogcm_ruv,l_obc_ogcm_ruvz,l_obc_ogcm_rc,       &
                          l_obc_mars,l_obc_diri,l_obc_char,l_obc_cycl,        &
                          l_obc_tide,l_obc_ogcm
   USE sflxatm,    ONLY : imeteo_dragtype,imeteo_exchtype
   USE trackpoints,ONLY : l_points,file_point
   USE output,     ONLY : file_output,l_out_nc4par,l_out_pack
   USE offline,    ONLY : path_offline,l_bz_offline
   USE comsiggen,  ONLY : b_sig,theta_sig,hc

   USE wave,       ONLY : file_wave,valfill_wave,name_wave_dir,name_wave_hs,  &
                          longitude_name_wave,latitude_name_wave,l_wave_rt0m1,&
                          l_wave_rtauwx,l_wave_rtauwy, l_wave_russ,l_wave_rfbb,&
                          l_wave_rtauox,l_wave_rtauoy,l_wave_rwlm,l_wave_rbhd,&
                          l_wave_rdir,l_wave_rhs,l_wave_rubot,l_wave_rfoc,    &
                          l_wave_rtus,l_wave_rtp,l_wave_rradstr,              &
                          name_wave_t0m1,l_wave_rtbb,  &
                          name_wave_ubot,name_wave_tp,name_wave_foc,          &
                          name_wave_tauwx,name_wave_tauwy,name_wave_tauox,    &
                          name_wave_tbbx,name_wave_tbby,name_wave_fbb,        &
                          name_wave_tauoy, name_wave_wlm,name_wave_bhd,       &
                          name_wave_sx,name_wave_sy,name_wave_sxy,            &
                          l_wave_marsgrid, name_wave_uubot, name_wave_vubot,  &
                          l_wave_ruubot,l_wave_rvubot,name_out_wave_dir,name_out_wave_hs, &
                          name_wave_ussx,name_wave_ussy,name_wave_tusx,name_wave_tusy,   &
                          name_out_wave_ubot,l_wave_rt0m1,name_out_wave_t0m1, &
                          name_out_wave_ubotx,name_out_wave_uboty,            &
                          name_out_wave_tauox,name_out_wave_tauoy,            &
                          name_out_wave_tauwx,name_out_wave_tauwy,            &
                          name_out_wave_wlm,name_out_wave_bhd,                &
                          name_out_wave_foc,name_out_wave_tp,                 &
                          name_out_wave_tbbx,name_out_wave_tbby,name_out_wave_fbb, &
                          name_out_wave_tusx,name_out_wave_tusy,              &
                          name_out_wave_ussx,name_out_wave_ussy,              &
                          name_out_wave_usx,name_out_wave_usy,name_out_wave_usz, &
                          name_out_wave_sx,name_out_wave_sy,                  &
                          riog_valid_min_wave_dir,riog_valid_max_wave_dir,    &
                          riog_valid_min_wave_hs,riog_valid_max_wave_hs,      &
                          riog_valid_min_wave_ubot,riog_valid_max_wave_ubot,  &
                          riog_valid_min_wave_tp,riog_valid_max_wave_tp,      &
                          riog_valid_min_wave_foc,riog_valid_max_wave_foc,    &
                          riog_valid_min_wave_t0m,riog_valid_max_wave_t0m,    &
                          riog_valid_min_wave_two,riog_valid_max_wave_two,    &
                          riog_valid_min_wave_wlm,riog_valid_max_wave_wlm,    &
                          riog_valid_min_wave_bhd,riog_valid_max_wave_bhd,    &
                          riog_valid_min_wave_stokes,riog_valid_max_wave_stokes

   USE toolmpi, ONLY: MASTER
   USE ionc4, ONLY : ionc4_openr,ionc4_read_time,ionc4_close

   IMPLICIT NONE

!! * Arguments
   INTEGER, INTENT(in)   ::   narg    ! number of argument (MARSWEB only)

!! * Local declarations
   CHARACTER(LEN=19)     :: tool_sectodat
   CHARACTER(LEN=4)      :: nraclu,obc2d_filename
   CHARACTER(LEN=19)     :: date_start,date_end,date_start3d,date_startagrif
   CHARACTER(LEN=19)     :: date_startdiag,date_enddiag
   CHARACTER(LEN=lchain) :: file_head,nverss,filepc
   CHARACTER(LEN=lchain) :: path_outputs,modele,logfilename
   INTEGER               :: i,j,l,nolu,eofcl,nlast,IERR_MPI
   INTEGER               :: ie,iw,js,jn,obcrk_num ! i,j indexes of obc for each file number obcrk_num
   REAL(KIND=rsh)        :: dtini,x
   REAL(KIND=rlg)        :: tool_datosec,tdb,tfi
   REAL(KIND=riolg)      :: finr,fisr,ger,gwr,dfir,dgr,dxr,dyr ! variables defining the area
   LOGICAL               :: ex


!! * define namelists

! namelists in paramain.txt
!--------------------------

   NAMELIST/nmlmain/ date_ref,comdir,jhemisph
   NAMELIST/nmlrestart/ name_in_dt,name_in_xe,name_in_xev,                &
                        name_in_u,name_in_v,name_in_uz,name_in_vz,        &
                        name_in_wz,name_in_temp,name_in_sal,name_in_psi,  &
                        name_in_ect,name_in_kz,name_in_nz,name_in_vis2d,  &
                        name_in_vis3d_xe,name_in_vis3d_phi ,              &
                        name_in_mxe,name_in_mhex,name_in_mhey,name_in_mkz,&
                        name_in_longitude,name_in_longitude_u,            &
                        name_in_longitude_v,name_in_longitude_f,          &
                        name_in_latitude,name_in_latitude_u,              &
                        name_in_latitude_v,name_in_latitude_f
   NAMELIST/nmloutput/ l_maxu,name_out_h0,name_out_hx,name_out_hy,        &
                       name_out_u,name_out_v,name_out_uz,name_out_vz,     &
                       name_out_vis2d,name_out_vishxe,name_out_vishphi,   &
                       name_out_nz,name_out_kz,name_out_bz,               &
                       name_out_mxe,name_out_mhex,name_out_mhey,          &
                       name_out_mkz,name_out_temp,name_out_sal,           &
                       name_out_xe,name_out_tosx,name_out_tosy,           &
                       name_out_mrad,name_out_mirlos,name_out_mlatlos,    &
                       name_out_msenslos,name_out_mevap,name_out_mrain,   &
                       name_out_mtosx,name_out_mtosy
   NAMELIST/nmlconvention/ riog_valid_min_h0,riog_valid_max_h0,           &
                           riog_valid_min_u,riog_valid_max_u,             &
                           riog_valid_min_temp,riog_valid_max_temp,       &
                           riog_valid_min_sal,riog_valid_max_sal,         &
                           riog_valid_min_xe,riog_valid_max_xe,           &
                           riog_valid_min_wstress,riog_valid_max_wstress, &
                           riog_valid_min_flux,riog_valid_max_flux,       &
                           riog_valid_min_visc,riog_valid_max_visc,       &
                           riog_valid_min_nz,riog_valid_max_nz,           &
                           riog_valid_min_kz,riog_valid_max_kz,           &
                           riog_valid_min_bz,riog_valid_max_bz,           &
                           riog_valid_min_dx,riog_valid_min_dy,           &
                           riog_valid_max_dx,riog_valid_max_dy,           &
                           cg_conv_title,cg_conv_comment,cg_conv_source,  &
                           cg_conv_area,cg_conv_institution,              &
                           cg_conv_institution_references,                &
                           cg_conv_data_centre,cg_conv_forecast_range,    &
                           cg_conv_data_centre_references,                &
                           cg_conv_forecast_type,cg_conv_field_type,      &
                           cg_conv_operational_status
   NAMELIST/nmladvtra/ qmax,qmaxz
   NAMELIST/nmlthermo/ rhoref,saliref,tetaref,l_stateeq_lin,              &
                       saliref_lin,tetaref_lin,rhoair,chp,coext
   NAMELIST/nmlinterp/ perchrono,tobs,l_champroche
   NAMELIST/nmlexpertonly/ sflx_rad_maxdepth,dt_minallow

! namelists in paracom.txt
!--------------------------

   NAMELIST/namdate/ date_start,date_end
   NAMELIST/namhead/ file_head,confname_mother,confname,suffix
   NAMELIST/nammeteoc/ l_meteo_stat,l_meteo_hom,wind_veloc,wind_dir,file_meteo, &
                       longitude_name_meteo,latitude_name_meteo
   NAMELIST/namtide/ l_tide_M2harm,l_tide_M2statcoef,tide_M2coef,           &
                      file_tide_M2coef,file_tide_harmcp,                    &
                      l_tide_admittance,l_genpot



! namelists in paraspec.txt
!--------------------------

   NAMELIST/namdate3d/ date_start3d,date_startagrif,date_startshift
   NAMELIST/nambathy/ file_bathy,l_bathy_meanlev,bathy_nbsmooth,mslshift,                    &
                      name_bathy_hx,name_bathy_hy,name_bathy_h0,name_bathy_meanlev,          &
                      l_bathy_comb,l_bathy_pit,hminim,hrdef,h0fond,dhj,fwetmin,               &
                      hminkxky,l_tssubadv_over_smalldepth,hm,l_bathy_save,file_bathy_child
   NAMELIST/namrestart/ l_initfromfile,file_init,l_init_rtime,l_init_rdt,   &
                        l_init_rssh,l_init_rbtvel,l_init_r3dvel,l_init_rwz, &
                        l_init_rturb,l_init_rsal,l_init_rtemp,l_init_restart
   NAMELIST/namobc/ l_obc_mars,l_obc_diri,l_obc_char,l_obc_cycl,l_obc_cycl_x,             &
                    l_obc_cycl_y,l_obc_tide,l_obc_ogcm,l_tide_harmcompo,l_obc2drank_save, &
                    l_obc_south,file_obc_ogcm_s,l_obc_north,file_obc_ogcm_n,l_obc_west,   &
                    file_obc_ogcm_w,l_obc_east,file_obc_ogcm_e,l_obc_ogcm_relax,obc_coefrel, &
                    l_obc_ogcm_rssh,l_obc_ogcm_rt,l_obc_ogcm_rs,l_obc_ogcm_ruv,           &
                    l_obc_ogcm_ruvz
   NAMELIST/nammeteo/ imeteo_dragtype,imeteo_exchtype,cds,paref,          &
                      l_sflx_rpa,l_sflx_rwind,l_sflx_rsat,l_sflx_rrh,     &
                      l_sflx_rsolar,l_sflx_rir,l_sflx_rrain,l_sflx_rfcc,  &
                      l_sflx_rath2m,l_sflx_rflux,                         &
                      name_sflx_sat,name_sflx_pa,name_sflx_rh,            &
                      name_sflx_windx,name_sflx_windy,name_sflx_solar,    &
                      name_sflx_ir,name_sflx_rain,name_sflx_fcc,          &
                      name_sflx_wtaux,name_sflx_wtauy,name_sflx_heat,     &
                      name_sflx_emp
   NAMELIST/namsaverestart/ l_saverestart_1file,l_saverestart_bydate,file_saverestart,saverestart_step
   NAMELIST/namoutput/ iscreen,iscreenlog,ierrorlog,iwarnlog,out_dir,file_output,l_out_nc4par,l_out_pack
   NAMELIST/namvisc/ sponge_nbcell,vismul,l_spongs,l_spongn,l_sponge,l_spongw, &
                     fvisc,l_smagor,vismin,vismax,cosmag
   NAMELIST/namturb/ z0bot,z0surf,l_z0b_hom,file_z0b,hminfrot,hmaxfrot,turb_nbeq,             &
                     turb_0eq_option,turb_2eq_option,bgdiff,nzinit,kzinit, &
                     l_stability
   NAMELIST/namfrot/ l_incbotstress,botstressmax
   NAMELIST/namdiff/ kx,ky
   NAMELIST/namgridsig/ l_equisig,sig,hc,b_sig,theta_sig,grid_angle,l_corio_cst
   NAMELIST/namriver/ file_river,icon
   NAMELIST/namparanum/ dtini,dtmin,dtmax,l_modele2d,cflcrt
   NAMELIST/namthermo/ l_sflx_radlossbot,l_sflx_solarcst,sflx_solarcst
   NAMELIST/namtraj/ file_trajec,itypetraj,ndtz
   NAMELIST/namdiag/l_diag,date_startdiag,date_enddiag,file_diag,l_points,file_point
   NAMELIST/namoffline/path_offline,l_bz_offline
   NAMELIST/namsurge/ l_tide_surge,file_tide_harmcp_surge,file_surge_meteo,file_surge_nometeo

! namelists in paraoasis.txt
!---------------------------

   NAMELIST/namoaocn/ l_cpl_ocn_maskt, l_cpl_ocn_masku, l_cpl_ocn_maskv, &
                      l_cpl_ocn_ssh, l_cpl_ocn_sst, l_cpl_ocn_u3d, l_cpl_ocn_v3d, &
                      l_cpl_ocn_u2d, l_cpl_ocn_v2d
   NAMELIST/namoaatm/ l_cpl_atm_fwsu, l_cpl_atm_fwsv, l_cpl_atm_heat, l_cpl_atm_rad, &
                      l_cpl_atm_ir, l_cpl_atm_lat, l_cpl_atm_sens, &
                      l_cpl_atm_evap, l_cpl_atm_rain, l_cpl_atm_pa
   NAMELIST/namoawave/l_cpl_wave_mask, l_cpl_wave_tmn, l_cpl_wave_hs, &
                      l_cpl_wave_dir, l_cpl_wave_bhd, l_cpl_wave_two, &
                      l_cpl_wave_uba, l_cpl_wave_foc, l_cpl_wave_taw, &
                      l_cpl_wave_fbb, l_cpl_wave_tbb, l_cpl_wave_lm, l_cpl_wave_stokes

! namelists in parasubs.txt
!--------------------------


! namelist generale pour traitement des variables de l ecosysteme
   NAMELIST/namecosys/ file_outflow,name_filesubs,date_startsub,l_obc_ogcm_rc,            &
                       l_restart_subs,l_out_log,l_bilan,name_bilfil,dtout_budget,         &
                       filevardiag,l_driv_cst, dt_driv,nb_border,                         &
                       l_cvrain_readfile,l_subflxatm_readfile,sflx_sub_atm_depth,         &
                       file_obc_subs_n,file_obc_subs_s,file_obc_subs_e,file_obc_subs_w,   &
                       l_subflxatm_xyt,file_flxatm_subs,date_start_budget

! namelists in paraconta.txt
!---------------------------



! namelists in parabiolo.txt
!---------------------------
! lues dans bioloinit



! namelists in parawave.txt
!---------------------------


   NAMELIST/nawinput/  file_wave,name_wave_dir,name_wave_hs,      &
                       longitude_name_wave,latitude_name_wave,valfill_wave,    &
                       name_wave_dir,name_wave_hs,name_wave_wlm,name_wave_tp,  &
                       name_wave_t0m1,name_wave_tauwx,name_wave_tauwy,         &
                       name_wave_foc,name_wave_tauox,name_wave_tauoy,          &
                       name_wave_bhd,name_wave_ubot, name_wave_uubot,          &
                       name_wave_vubot,name_wave_tusx, name_wave_tusy,         &
                       name_wave_ussx,name_wave_ussy,name_wave_sx,             &
                       name_wave_sy,name_wave_sxy,  &
                       l_wave_rdir,l_wave_rhs,l_wave_rfoc,l_wave_rtauox,       &
		       l_wave_rtauoy,l_wave_rbhd,l_wave_rubot,l_wave_ruubot,   &
                       l_wave_rvubot,l_wave_rfbb,l_wave_rtbb,l_wave_rtp,       &
                       l_wave_rt0m1,l_wave_rtauwx,l_wave_rtauwy,l_wave_rwlm,   &
                       l_wave_rbhd,l_wave_rtus,l_wave_russ,l_wave_rradstr,     & 
                       l_wave_marsgrid
   NAMELIST/nawoutput/ name_out_wave_dir,name_out_wave_hs,name_out_wave_ubot,  &
                       name_out_wave_ubotx,name_out_wave_uboty,                &
                       name_out_wave_foc,name_out_wave_tp,name_out_wave_t0m1,  &
                       name_out_wave_tauox,name_out_wave_tauoy,              &
                       name_out_wave_tauwx,name_out_wave_tauwy,              &
                       name_out_wave_tbbx,name_out_wave_tbby,name_out_wave_fbb, &
                       name_out_wave_wlm,name_out_wave_bhd,                  &
                       name_out_wave_tusx,name_out_wave_tusy,                &
                       name_out_wave_ussx,name_out_wave_ussy,                &
                       name_out_wave_usx,name_out_wave_usy,name_out_wave_usz, &
                       name_out_wave_sx,name_out_wave_sy,                    &
                       riog_valid_min_wave_dir,riog_valid_max_wave_dir,      &
                       riog_valid_min_wave_hs,riog_valid_max_wave_hs,        &
                       riog_valid_min_wave_ubot,riog_valid_max_wave_ubot,    &
                       riog_valid_min_wave_tp,riog_valid_max_wave_tp,        &
                       riog_valid_min_wave_foc,riog_valid_max_wave_foc,      &
                       riog_valid_min_wave_t0m,riog_valid_max_wave_t0m,      &
                       riog_valid_min_wave_two,riog_valid_max_wave_two,      &
                       riog_valid_min_wave_wlm,riog_valid_max_wave_wlm,      &
                       riog_valid_min_wave_bhd,riog_valid_max_wave_bhd,      &
                       riog_valid_min_wave_stokes,riog_valid_max_wave_stokes


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

   !PRINT*, 'enter INIT routine'

! --------------
! Read namelists
! --------------

! read namelist of general parameters
!------------------------------------

   filepc='./paramain.txt'
!  if(nracin/='zzzz') filepc='../'//nracin//'.paramain.txt'
   OPEN(50,file=filepc,status='old',form='formatted',access='sequential')
   READ(50,nmlmain)
   READ(50,nmlrestart)
   READ(50,nmloutput)
   READ(50,nmlconvention)
   READ(50,nmladvtra)
   READ(50,nmlthermo)
   READ(50,nmlinterp)
   READ(50,nmlexpertonly)
   CLOSE(50)

! read namelist common to each rank
!----------------------------------

   filepc=TRIM(comdir)//'paracom.txt'
   !IF (.NOT. .TRUE.) THEN
     !  filepc=TRIM(comdir)//'paracom'//TRIM(Agrif_Cfixed())//'.txt'
     !  INQUIRE(file=filepc,exist=ex)
     !  IF (.NOT. ex) filepc=TRIM(comdir)//'paracom.txt'
   !ENDIF

   OPEN(50,file=filepc,status='old',form='formatted',access='sequential')
   READ(50,namdate)
   READ(50,namhead)
   READ(50,nammeteoc)
   READ(50,namtide)

   CLOSE(50)

! read namelist specific to each rank
!------------------------------------

   filepc='./paraspec.txt'
   !IF (.NOT. .TRUE.) filepc='./paraspec'//TRIM(Agrif_Cfixed())//'.txt'
   OPEN(50,file=filepc,status='old',form='formatted', access='sequential')
   READ(50,namdate3d)
   READ(50,nambathy)
   READ(50,namrestart)
   READ(50,namobc)
   READ(50,nammeteo)
   READ(50,namsaverestart)
   READ(50,namoutput)
   READ(50,namvisc)
   READ(50,namturb)
   READ(50,namfrot)
   READ(50,namdiff)
   READ(50,namgridsig)
   READ(50,namriver)
   READ(50,namparanum)
   READ(50,namthermo)
   READ(50,namtraj)
   READ(50,namdiag)
   READ(50,namoffline)
   READ(50,namsurge)

   CLOSE(50)


! read namelist relative to the OASIS coupling
! --------------------------------------------


   CALL init_comoasis


! read namelist relative to the ecosystem (biology, sedimentology, contaminants...)
! ---------------------------------------------------------------------------------


   filepc='./parasubs.txt'
   !IF (.NOT. .TRUE.) filepc='./parasubs'//TRIM(Agrif_Cfixed())//'.txt'
   OPEN(50,file=filepc,status='old',form='formatted', access='sequential')
   READ(50,namecosys)

   CLOSE(50)                                  ! end reading of namelists


! read namelist relative to contaminants
! ---------------------------------------------------------------------------------



! read namelist relative to the sedimentology
! ---------------------------------------------------------------------------------


   CALL MUSTANG_param('r')




! read namelist relative to the biology
! ---------------------------------------------------------------------------------



! read the namelist relative to the waves
! ---------------------------------------------------------------------------------


   filepc='./parawave.txt'
   !IF (.NOT. .TRUE.) filepc='./parawave'//TRIM(Agrif_Cfixed())//'.txt'
   OPEN(50,file=filepc,status='old',form='formatted', access='sequential')
   READ(50,nawinput)
   READ(50,nawoutput)
   CLOSE(50)                                  ! end reading of namelist


! ----------------------------------------------------------------------
! save the parametrizations inside a log file : simu.log
! initialize a listing file containing maximum velocities, time steps...
! ----------------------------------------------------------------------

   IF (narg == 0) THEN
     logfilename= 'simu_'//TRIM(suffix)//'.log'
     !IF (.NOT. .TRUE.) THEN
     ! logfilename= 'simu_'//TRIM(suffix)//'_'//TRIM(Agrif_Cfixed())//'.log'
     !ENDIF
     IF (MASTER) THEN
     IF (iscreenlog/=6) OPEN(unit=iscreenlog,file=logfilename,status="replace")
     ENDIF
     logfilename='error_'//TRIM(suffix)//'.log'
     !IF (.NOT. .TRUE.) logfilename= 'error_'//TRIM(suffix)//'_'//TRIM(Agrif_Cfixed())//'.log'
     IF (ierrorlog/=6) OPEN(unit=ierrorlog,file=logfilename,status="replace")
     logfilename='warning_'//TRIM(suffix)//'.log'
     !IF (.NOT. .TRUE.) logfilename= 'warning_'//TRIM(suffix)//'_'//TRIM(Agrif_Cfixed())//'.log'
     IF (iwarnlog/=6) OPEN(unit=iwarnlog,file=logfilename,status="replace")
   ELSE IF (narg == 2) THEN
     CALL getarg(2,path_outputs)
     IF (iscreen/=6) THEN
       modele=TRIM(path_outputs)//'/listing_'//TRIM(suffix)//'.dat'
       OPEN(iscreen,file=modele)
     END IF
     IF (MASTER) THEN
     IF (iscreenlog/=6) THEN
       modele=TRIM(path_outputs)//'/simu_'//TRIM(suffix)//'.log'
       OPEN(iscreenlog,file=modele)
     END IF
     ENDIF
     IF (ierrorlog/=6) THEN
       modele=TRIM(path_outputs)//'/error_'//TRIM(suffix)//'.log'
       OPEN(ierrorlog,file=modele)
     END IF
     IF (iwarnlog/=6) THEN
       modele=TRIM(path_outputs)//'/warning_'//TRIM(suffix)//'.log'
       OPEN(iwarnlog,file=modele)
     END IF
   ELSE IF (narg == 3) THEN
     CALL getarg(3,path_outputs)
     IF (iscreen/=6) THEN
       modele=TRIM(path_outputs)//'/listing_'//TRIM(suffix)//'.dat'
       OPEN(iscreen,file=modele)
     END IF
     IF (MASTER) THEN
     IF (iscreenlog/=6) THEN
       modele=TRIM(path_outputs)//'/simu_'//TRIM(suffix)//'.log'
       OPEN(iscreenlog,file=modele)
     END IF
     ENDIF
     IF (ierrorlog/=6) THEN
       modele=TRIM(path_outputs)//'/error_'//TRIM(suffix)//'.log'
       OPEN(ierrorlog,file=modele)
     END IF
     IF (iwarnlog/=6) THEN
       modele=TRIM(path_outputs)//'/warning_'//TRIM(suffix)//'.log'
       OPEN(iwarnlog,file=modele)
     END IF
   END IF
 

   IF (MASTER) THEN

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) '**************** INIT.F90 ************************'
   WRITE(iscreenlog,*) '**************************************************'

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '!simulation'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,namhead)
   WRITE(iscreenlog,*) 'IMIN = ',imin,' IMAX = ',imax
   WRITE(iscreenlog,*) 'JMIN = ',jmin,' JMAX = ',jmax
   WRITE(iscreenlog,*) 'KMAX = ',kmax
   WRITE(iscreenlog,namgridsig)

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! dates'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,nmlmain)
   WRITE(iscreenlog,namdate)
   WRITE(iscreenlog,namdate3d)

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! numerical parameters'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,namparanum)
   WRITE(iscreenlog,namfrot)
   WRITE(iscreenlog,nmladvtra)

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! bathymetry'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,nambathy)

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! meteo'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,nammeteoc)
   WRITE(iscreenlog,nammeteo)

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! restart'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,nmlrestart)
   WRITE(iscreenlog,namrestart)

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! output'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,nmloutput)
   WRITE(iscreenlog,nmlconvention)
   WRITE(iscreenlog,namoutput)
   WRITE(iscreenlog,namsaverestart)

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! viscosity'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,namvisc)

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! diffusion'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,namdiff)

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! turbulence'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,namturb)

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! thermic parameters'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,nmlthermo)
   WRITE(iscreenlog,namthermo)

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! special parametrization (experts only)'
   WRITE(iscreenlog,*) '****************************************'
   WRITE(iscreenlog,nmlexpertonly)

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! time interpolation of data sets'
   WRITE(iscreenlog,*) '*********************************'
   WRITE(iscreenlog,nmlinterp)

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! diagnoses'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,namdiag)

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! trajectories'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,namtraj)

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! tides'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,namtide)

   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! open boundary conditions'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,namobc)



   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! overtide functionality'
   WRITE(iscreenlog,*) '************************'
   WRITE(iscreenlog,namsurge)



   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! rivers'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,namriver)


   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! substances'
   WRITE(iscreenlog,*) '***********************'
   WRITE(iscreenlog,namecosys)



   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '! sedimentology MUSTANG'
   WRITE(iscreenlog,*) '***********************'
   CALL MUSTANG_param('w')








   ENDIF

   !IF (.NOT. .TRUE.) filename= 'listing_zoom'//TRIM(Agrif_Cfixed())
   !IF (.NOT. .TRUE. .AND. iscreen/=6) OPEN(unit=iscreen,file=filename,status="replace")


! --------------------------------------------------------
! Read head.CONF file
! Initialize variables defining the current areas
! and managing the simultaneous generation of 2d obc files
! --------------------------------------------------------

   !IF (.TRUE.) THEN

! Check all path names end with a "/".
     nlast=LEN(TRIM(comdir))
     IF (comdir(nlast:nlast)/='/') comdir=TRIM(comdir)//'/'

     nlast=LEN(TRIM(out_dir))
     IF (out_dir(nlast:nlast)/='/') out_dir=TRIM(out_dir)//'/'

     nlast=LEN(TRIM(path_offline))
     IF (path_offline(nlast:nlast)/='/') path_offline=TRIM(path_offline)//'/'

   IF (.NOT.l_testcase) THEN

     IF (MASTER) &
     WRITE(iscreenlog,*) 'Number of the current rank of the configuration :', no

! if 2D open boundary conditions must be save only
! read all the lines in head.conf in order to get the number of obc files to save : obc_filemax
! WARNING : the four last numbers on line relative to the rank no correspond to the location of obc
! for each child grid (rank no+1)
! obc2d_filename is the name of the config for which 2D obc are generated
! -------------------------------------------------------------------

     IF (l_obc2drank_save) THEN
       eofcl=0
       OPEN(99,file=file_head,status='old',form='formatted')

       READ(99,3001,iostat=eofcl)nraclu,nolu,x,x,x,x,dfip,dgp,x,x,i,i,i,i,i,i,obc2d_filename
       obcrk_maxnbfile=0                       ! initialize obcrk_maxnbfile to zero
       DO WHILE (nolu<=no .AND. eofcl==0)    ! no need to read ranks larger than no
         IF (nolu == no) THEN
           obcrk_maxnbfile=obcrk_maxnbfile+1
           IF (MASTER) THEN
           WRITE(iscreenlog,*) ' '
           WRITE(iscreenlog,*) obcrk_maxnbfile,'2D open boundary conditions saved for configuration :',&
                               TRIM(obc2d_filename)
           ENDIF
         END IF
         READ(99,3001,iostat=eofcl)nraclu,nolu,x,x,x,x,dfip,dgp,x,x,i,i,i,i,i,i,obc2d_filename ! reads nolu
       END DO
       CLOSE(99)



! allocate dimension of ie, iw, js, jn and obcrk_nameconf
       ALLOCATE(obcrk_ie(obcrk_maxnbfile))
       ALLOCATE(obcrk_iw(obcrk_maxnbfile))
       ALLOCATE(obcrk_js(obcrk_maxnbfile))
       ALLOCATE(obcrk_jn(obcrk_maxnbfile))
       ALLOCATE(obcrk_nameconf(obcrk_maxnbfile))
     END IF

! get informations from head.CONF file
! ------------------------------------

     eofcl=0
     OPEN(99,file=file_head,status='old',form='formatted')    ! open filehead again
     obcrk_num=0
     READ(99,3001,iostat=eofcl)nraclu,nolu,finr,fisr,ger,gwr,dfir,dgr,dxr,dyr,&
                               imaxac,jmaxac,jn,js,ie,iw,obc2d_filename         ! read file_head

! get indexes of open boundaries inside the mother grid (rank no-1)
! -----------------------------------------------------------------

     DO WHILE (nolu<no)
       IF ((nolu .EQ. no-1) .AND. (obc2d_filename .EQ. nrac)) THEN
! read parameters in the line where nrac (name of the current conf no)
! fits with obc2d_filename (name of the conf for which obc file is generated from rank no-1)
         jnac=jn
         jsac=js
         ieac=ie
         iwac=iw
         dfip=dfir
         dgp=dgr
       END IF
       READ(99,3001,iostat=eofcl)nraclu,nolu,finr,fisr,ger,gwr,dfir,dgr,dxr,dyr,&
                                 imaxac,jmaxac,jn,js,ie,iw,obc2d_filename         ! read file_head
     END DO

! reads lines corresponding to rank no
! initializes location of the current grid : fisac,finac,gwac,geac,dfiac,dgac,dxac,dyac
! initializes indexes for which 2D open boundary conditions have to be saved
!---------------------------------------------------------------

     DO WHILE (nolu<=no .AND. eofcl==0)   ! reads lines corresponding to rank no
! verifications
       IF (no /= nolu.OR.nrac /= nraclu) THEN
         WRITE(ierrorlog,*) 'ERROR: Mistake inside the "HEAD.conf" file or in parameters.F90* file'
         WRITE(ierrorlog,*) 'The number of rank, no = ',no,' is different from the number in head.CONF file, nolu = ',nolu
         WRITE(ierrorlog,*) 'OR'
         WRITE(ierrorlog,*) 'The root name of the configuration nrac = ',nrac, &
                            ' is different from the root in head.CONF file, nraclu =',nraclu
         WRITE(ierrorlog,*) 'Simulation stopped'
         CALL MPI_FINALIZE(IERR_MPI)
         STOP
       END IF

       IF ( imaxac/=imax .OR. jmaxac/=jmax ) THEN
         WRITE(ierrorlog,*) 'ERROR: Mistake inside the "HEAD.conf" file for rank ',no
         WRITE(ierrorlog,*) 'imax = ',imax,' imax lu = ',imaxac, ' jmax = ',jmax,' jmax lu = ',jmaxac
         WRITE(ierrorlog,*) 'Simulation stopped'
         CALL MPI_FINALIZE(IERR_MPI)
         STOP
       END IF
       IF ((js>jn).OR.(js<jmin).OR.(jn>jmax)) THEN
         WRITE(ierrorlog,*) 'ERROR: Variable "js" and (or) "jn" is not correct.'
         WRITE(ierrorlog,*) 'Simulation stopped'
         CALL MPI_FINALIZE(IERR_MPI)
         STOP
       END IF
       IF ((iw>ie).OR.(iw<imin).OR.(ie>imax)) THEN
         WRITE(ierrorlog,*) 'ERROR: WARNING : Variable "iw" and (or) "ie" is not correct.'
         WRITE(ierrorlog,*) 'Simulation stopped'
         CALL MPI_FINALIZE(IERR_MPI)
         STOP
       END IF


       finac=finr
       fisac=fisr
       geac=ger
       gwac=gwr
       dfiac=dfir
       dgac=dgr
       dxac=dxr
       dyac=dyr
       IF (l_obc2drank_save) THEN
         obcrk_num=obcrk_num+1
         obcrk_ie(obcrk_num)=ie
         obcrk_iw(obcrk_num)=iw
         obcrk_js(obcrk_num)=js
         obcrk_jn(obcrk_num)=jn
         obcrk_nameconf(obcrk_num)=obc2d_filename
       END IF
       IF (l_obc2drank_save) THEN
         READ(99,3001,iostat=eofcl)nraclu,nolu,finr,fisr,ger,gwr,dfir,dgr,dxr,dyr,&
                                   imaxac,jmaxac,jn,js,ie,iw,obc2d_filename         ! read file_head
       ELSE
         nolu=nolu+1
       END IF
     END DO
     CLOSE(99)

! name of file saving 2D open boundary conditions
!------------------------------------------------

     WRITE(nverss,'(i1)') no+1

     fileclsv =  '../'//TRIM(comdir)//'inputs/obc2d_'//TRIM(confname)//'.rank'//TRIM(nverss)//'_'


! name of file from which 2D open boundary conditions are read
!-------------------------------------------------------------

     WRITE(nverss,'(i1)') no

     filecl =  '../'//TRIM(comdir)//'inputs/obc2d_'//TRIM(confname_mother)//'.rank'

     filecl = TRIM(filecl)//TRIM(nverss)//'_'//TRIM(nrac)//'.'//TRIM(suffix)

! name of saving file for restart
!--------------------------------

     file_saverestart=TRIM(out_dir)//TRIM(file_saverestart)

   END IF

   !ENDIF

! ---------------------------------
! Initializations and verifications
! ---------------------------------

! name of file saving diagnoses
!------------------------------

   file_diag=TRIM(out_dir)//TRIM(file_diag)//TRIM(suffix)//'.txt'

   IF (l_saverestart_1file.AND.l_saverestart_bydate) THEN
     WRITE(ierrorlog,*) 'ERROR : Only one type of saving must be chosen.'
     WRITE(ierrorlog,*) 'l_saverestart_1file  : a unique saving file (new saving replaces old one)'
     WRITE(ierrorlog,*) 'l_saverestart_bydate : a new file for each date'
     WRITE(ierrorlog,*) 'Simulation stopped'
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   END IF

   IF (l_init_restart.AND.(.NOT.l_initfromfile)) THEN
     WRITE(ierrorlog,*) 'ERROR : An initialisation file must be read if you want the simulation'
     WRITE(ierrorlog,*) 'to restart from exact results (for reproductibility purpose).'
     WRITE(ierrorlog,*) 'You have to set the following parameters in the "paraspec.txt" file :' 
     WRITE(ierrorlog,*) 'l_initfromfile =  .true. '
     WRITE(ierrorlog,*) 'file_init = initial conditions filename'
     WRITE(ierrorlog,*) 'Simulation stopped'
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   END IF

   IF (l_init_restart.AND.(.NOT.l_init_rtime .OR. .NOT.l_init_rdt .OR. .NOT.l_init_rssh .OR. .NOT.l_init_rbtvel      &
                           .OR. .NOT.l_init_r3dvel .OR. .NOT.l_init_rwz .OR. .NOT.l_init_rsal .OR. .NOT.l_init_rtemp &
                           .OR. .NOT.l_init_rturb ) ) THEN
     WRITE(ierrorlog,*) 'ERROR : In order to use l_init_restart functionality'
     WRITE(ierrorlog,*) 'and make sure to get correct results,'
     WRITE(ierrorlog,*) 'You have to set to .true. the following parameters in the "paraspec.txt" file :' 
     WRITE(ierrorlog,*) 'l_init_rtime =  .true. '
     WRITE(ierrorlog,*) 'l_init_rdt =  .true. '
     WRITE(ierrorlog,*) 'l_init_rssh =  .true. '
     WRITE(ierrorlog,*) 'l_init_rbtvel =  .true. '
     WRITE(ierrorlog,*) 'l_init_r3dvel =  .true. '
     WRITE(ierrorlog,*) 'l_init_rwz =  .true. '
     WRITE(ierrorlog,*) 'l_init_rsal =  .true. '
     WRITE(ierrorlog,*) 'l_init_rtemp =  .true. '
     WRITE(ierrorlog,*) 'l_init_rturb =  .true. '
     WRITE(ierrorlog,*) 'Simulation stopped'
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   END IF

! make sure you force the model with the correct atmospheric forcing
!-------------------------------------------------------------------

   IF ( .NOT.l_meteo_stat .AND. .NOT. l_meteo_hom ) THEN
     IF ( l_sflx_rflux .AND. (.NOT.l_sflx_rwind .OR. .NOT.l_sflx_rsolar .OR. .NOT.l_sflx_rrain) ) THEN
       WRITE(ierrorlog,*) ''
       WRITE(ierrorlog,*) 'Since l_meteo_stat and l_meteo_hom are .false. and l_sflx_rflux is true'
       WRITE(ierrorlog,*) 'you want to read the atmospheric forcing in flux form'
       WRITE(ierrorlog,*) 'Specify'
       WRITE(ierrorlog,*) ' - in paracom.txt, the name of the file : file_meteo '
       WRITE(ierrorlog,*) ' - in paraspec.txt, l_sflx_rwind=.true. (for the wind stress actually)'
       WRITE(ierrorlog,*) '                    l_sflx_rsolar=.true.'
       WRITE(ierrorlog,*) '                    l_sflx_rrain=.true. (for EmP field actually)'
       WRITE(ierrorlog,*) '         as well as name_sflx_wtaux,name_sflx_wtauy,name_sflx_heat,name_sflx_solar and name_sflx_emp'
       WRITE(ierrorlog,*) ''
       WRITE(ierrorlog,*) 'Simulation stopped'
       WRITE(ierrorlog,*) ''
       CALL MPI_FINALIZE(IERR_MPI)
       STOP
     ELSE IF ( .NOT.l_sflx_rpa .AND. .NOT.l_sflx_rwind .AND. .NOT.l_sflx_rsat .AND. .NOT.l_sflx_rrh .AND. &
          .NOT.l_sflx_rsolar .AND. .NOT.l_sflx_rir .AND. .NOT.l_sflx_rrain .AND. .NOT.l_sflx_rfcc ) THEN
       WRITE(ierrorlog,*) ''
       WRITE(ierrorlog,*) 'Since l_meteo_stat and l_meteo_hom are .false.'
       WRITE(ierrorlog,*) 'You must read the atmospheric forcing from a netcdf data set'
       WRITE(ierrorlog,*) 'and specify'
       WRITE(ierrorlog,*) ' - in paracom.txt, the name of the file : file_meteo '
       WRITE(ierrorlog,*) ' - in paraspec.txt, l_sflx_rwind=.true. at least'
       WRITE(ierrorlog,*) '                    and the name of the wind variables : name_sflx_windx and name_sflx_windy'
       WRITE(ierrorlog,*) ''
       WRITE(ierrorlog,*) 'Think about the parametrization of the rest of the atmospherical forcing'
       WRITE(ierrorlog,*) '(realistic atmospheric pressure, solar, infrared thermal and turbulent forcings'
       WRITE(ierrorlog,*) 'or constant solar flux l_sflx_solarcst,sflx_solarcst)'
       WRITE(ierrorlog,*) 'see paraspec.txt and the available cpp keys'
       WRITE(ierrorlog,*) ''
       WRITE(ierrorlog,*) 'If you want to run with a homogeneous (uniform in space)'
       WRITE(ierrorlog,*) 'and/or stationar (constant in time)'
       WRITE(ierrorlog,*) ' - in paracom.txt, specify l_meteo_stat, l_meteo_hom, wind_veloc, wind_dir'
       WRITE(ierrorlog,*) ' - in paraspec.txt, set l_sflx_solarcst=.true. and specify the mount of the solar flux (sflx_solarcst)'
       WRITE(ierrorlog,*) ''
       WRITE(ierrorlog,*) 'Simulation stopped'
       WRITE(ierrorlog,*) ''
       CALL MPI_FINALIZE(IERR_MPI)
       STOP
     ELSE IF ( .NOT.l_sflx_rwind ) THEN
       IF (MASTER) THEN
       WRITE(iwarnlog,*) ''
       WRITE(iwarnlog,*) 'Since l_meteo_stat and l_meteo_hom are .false.'
       WRITE(iwarnlog,*) 'You must read a realistic wind field from an atmospheric netcdf data set'
       WRITE(iwarnlog,*) 'and specify'
       WRITE(iwarnlog,*) ' - in paracom.txt, the name of the file : file_meteo '
       WRITE(iwarnlog,*) ' - in paraspec.txt, l_sflx_rwind=.true. '
       WRITE(iwarnlog,*) '                    and the name of the wind variables : name_sflx_windx and name_sflx_windy'
       WRITE(iwarnlog,*) ''
       WRITE(iwarnlog,*) 'Think about the parametrization of the rest of the atmospherical forcing'
       WRITE(iwarnlog,*) '(realistic atmospheric pressure, solar, infrared thermal and turbulent forcings'
       WRITE(iwarnlog,*) 'or constant solar flux l_sflx_solarcst,sflx_solarcst)'
       WRITE(iwarnlog,*) 'see paraspec.txt and the available cpp keys'
       WRITE(iwarnlog,*) ''
       WRITE(iwarnlog,*) 'Simulation stopped'
       WRITE(iwarnlog,*) ''
!CALL MPI_FINALIZE(IERR_MPI)
!STOP
       ENDIF
     END IF
   END IF
   IF ( l_meteo_stat .AND. l_meteo_hom ) THEN
     IF (.NOT. l_sflx_solarcst) THEN
       WRITE(ierrorlog,*) ''
       WRITE(ierrorlog,*) 'To use a constant wind forcing, you must prescribe a constant thermal forcing also'
       WRITE(ierrorlog,*) '- in paraspec.txt, specify l_sflx_solarcst=.true. and the value of sflx_solarcst'
       WRITE(ierrorlog,*) ''
       WRITE(ierrorlog,*) 'Simulation stopped'
       WRITE(ierrorlog,*) ''
       CALL MPI_FINALIZE(IERR_MPI)
       STOP
     END IF
   END IF

   IF ( l_sflx_rflux) THEN
     name_sflx_windx = name_sflx_wtaux
     name_sflx_windy = name_sflx_wtauy
     name_sflx_ir    = name_sflx_heat
     name_sflx_rain  = name_sflx_emp
     l_sflx_rir=.true.   ! to read the non solar component of the flux
   END IF

! make sure the surface drag coefficient is not set to zero
!------------------------------------------
   IF (turb_nbeq==2 .AND. z0surf ==0.0_rsh) THEN
     WRITE(ierrorlog,*) ''
     WRITE(ierrorlog,*) 'z0surf ==0.0_rsh !!!'
     WRITE(ierrorlog,*) 'You must specify a surface drag coefficient (in m)'
     WRITE(ierrorlog,*) 'It is constant (at the moment) and equal to 1.6 * Hs (the significant wave height)'
     WRITE(ierrorlog,*) ''
     WRITE(ierrorlog,*) 'Simulation stopped'
     WRITE(ierrorlog,*) ''
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   END IF
   IF (turb_nbeq>1 .AND. z0bot ==0.0_rsh) THEN
     WRITE(ierrorlog,*) ''
     WRITE(ierrorlog,*) 'z0bot ==0.0_rsh !!!'
     WRITE(ierrorlog,*) 'You must specify a bottom drag coefficient (of the order of a few mm)'
     WRITE(ierrorlog,*) ''
     WRITE(ierrorlog,*) 'Simulation stopped'
     WRITE(ierrorlog,*) ''
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   END IF


! make sure you have no sponge when using smagorinsky viscosity
!--------------------------------------------------------------



! make sure l_modele2d is correctly defined
!------------------------------------------

   IF ((kmax==1).AND.(.NOT.(l_modele2d))) THEN
!     l_modele2d=.true.
!     PRINT*,'variable l_modele2d should be set to .true. as kmax=1'
!     PRINT*,'corrected directly in routine init.F90'
     WRITE(ierrorlog,*) 'ERROR: kmax = ',kmax,', l_modele2d = ',l_modele2d
     WRITE(ierrorlog,*) 'The boolean "l_modele2d" must be set to .true. '
     WRITE(ierrorlog,*) 'as there is only 1 vertical level.'
     WRITE(ierrorlog,*) 'Simulation stopped'
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   END IF



! times and time step
!--------------------

   IF (l_initfromfile .AND. l_init_rtime) THEN
     CALL ionc4_openr(file_init,l_out_nc4par=.true.)
     CALL ionc4_read_time(file_init,1,t)
     tdeb=t
     CALL ionc4_close(file_init)
   ELSE
     tdeb    =tool_datosec(date_start)
   END IF
   tfin    =tool_datosec(date_end)
   IF (l_obc2drank_save) tfin=tfin+2*REAL(6-no,rlg)*dtmax
   tdeb3d  =tool_datosec(date_start3d)

   IF ( .NOT.l_initfromfile .OR. (l_initfromfile .AND. .NOT.l_init_rtime) ) t=tdeb+date_startshift*86400.0_rlg
   dt=dtini
   dtmaxi = dtmax
   saverestart_step=saverestart_step*86400.0_rlg
   tsauv=tdeb + saverestart_step + date_startshift*86400.0_rlg
   tobs=tobs*3600.0_rlg
   tdebobs=t


! initialisation of cdate, ijour, imois...
! for outputs and tide use
!-----------------------------------------

   cdate = tool_sectodat(t)
   CALL tool_decompdate(cdate,ijour,imois,ian,iheure,iminu,isec)

! checking...
!------------

   IF (.NOT.l_modele2d .AND. tdeb3d + 1.0_rlg <tdeb) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'ERROR : Variable "date_start3d" is uncorrect (before "date_start").'
     WRITE(ierrorlog,*) 'dates summary'
     WRITE(ierrorlog,*) 't =',cdate,'( date_start = ',date_start,' but the time departure can be read from the restart file'
     WRITE(ierrorlog,*) ' or shifted throughout date_startshift variable)'
     WRITE(ierrorlog,*) 'date_start3d =',date_start3d
     WRITE(ierrorlog,*) 'If you want to run in 3D, please define "date_start3d" properly'
     WRITE(ierrorlog,*) 'If you want to run in 2D, please set "l_modele2d=.true."'
     WRITE(ierrorlog,*) 'The simulation is stopped'
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
!     WRITE(iwarnlog,*) '3D simulations will start '
!     WRITE(iwarnlog,*) 'at the beginning of the run "date_start": ', date_start
!     WRITE(iwarnlog,*) 'or at the date of the restart file.'
!     tdeb3d=MAX(tdeb,tdeb3d)
   END IF
   IF (.NOT.l_modele2d .AND. tdeb3d>tfin) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'ERROR : Variable "date_start3d" is uncorrect (after "date_end").'
     WRITE(ierrorlog,*) 'If you want to run in 3D, please define "date_start3d" properly'
     WRITE(ierrorlog,*) 'If you want to run in 2D, please set "l_modele2d=.true."'
     WRITE(ierrorlog,*) 'The simulation is stopped'
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   END IF
   tdebagrif=tool_datosec(date_startagrif)
   tdebdiag=tool_datosec(date_startdiag)
   tfindiag=tool_datosec(date_enddiag)
   IF (l_modele2d) THEN
     IF (tdebdiag + 1.0_rlg <tdeb .AND. l_diag) THEN
       IF (MASTER) THEN
       WRITE(iwarnlog,*) ' '
       WRITE(iwarnlog,*) 'WARNING: Variable "date_startdiag" is uncorrect.'
       WRITE(iwarnlog,*) '2D simulations start on ', date_start
       WRITE(iwarnlog,*) 'Consequently, diagnoses are also saved since ', date_start
       ENDIF
       tdebdiag=MAX(tdebdiag,tdeb)
     END IF
     IF (tfindiag<tdeb .AND. l_diag) THEN
       IF (MASTER) THEN
       WRITE(iwarnlog,*) ' '
       WRITE(iwarnlog,*) 'WARNING: Variable "date_enddiag" is uncorrect.'
       WRITE(iwarnlog,*) '2D simulations end on ',date_end
       WRITE(iwarnlog,*) 'Consequently, diagnoses are also saved until ', date_end
       ENDIF
       tfindiag=MAX(tfindiag,tfin)
     END IF
   ELSE
     IF (tdebdiag + 1.0_rlg < tdeb3d .AND. l_diag) THEN
       IF (MASTER) THEN
       WRITE(iwarnlog,*) ' '
       WRITE(iwarnlog,*) 'WARNING: Variable "date_startdiag" is uncorrect.'
       WRITE(iwarnlog,*) '3D simulations start on ', date_start3d
       WRITE(iwarnlog,*) 'Consequently, diagnoses are also saved since ', date_start3d
       ENDIF
       tdebdiag=MAX(tdebdiag,tdeb3d)
     END IF
     IF (tfindiag<tdeb3d .AND. l_diag) THEN
       IF (MASTER) THEN
       WRITE(iwarnlog,*) ' '
       WRITE(iwarnlog,*) 'WARNING: Variable "date_enddiag" is uncorrect.'
       WRITE(iwarnlog,*) '3D simulations end on ', date_end
       WRITE(iwarnlog,*) 'Consequently, diagnoses are also saved until ', date_end
       ENDIF
       tfindiag=MAX(tfindiag,tfin)
     END IF
   END IF

   tdebsubs=tool_datosec(date_startsub)
   IF (tdebsubs + 1.0_rlg <tdeb3d) THEN
     IF (MASTER) THEN
     WRITE(iwarnlog,*) ' '
     WRITE(iwarnlog,*) 'WARNING: Variable "date_startsub" is uncorrect.'
     WRITE(iwarnlog,*) 'Biologic or sedimentologic simulations will start '
     WRITE(iwarnlog,*) 'at the begining of 3D simulations: ', date_start3d  
     ENDIF
   END IF
   tdebsubs=MAX(tdebsubs,tdeb3d)
   IF(l_bilan) THEN
     tdeb_budget=tool_datosec(date_start_budget)
     iout_budget=-1
   ENDIF


! initialize time vector tchrono relative to
! river dischage and outflow events
!-----------------------------------------

   l_chrono=.FALSE.
   perchrono=perchrono*3600.0_rlg
!   IF (l_meteo_hom .and. .not.l_meteo_stat) perchrono=3.0_rlg*3600.0_rlg
   IF (l_initfromfile .AND. l_init_rtime) THEN
     tdb=tdeb-perchrono
   ELSE
     tdb=tdeb+date_startshift*86400.0_rlg-perchrono
   END IF
   tfi=tfin+perchrono
   nvalmax=NINT((tfi-tdb)/perchrono+0.5_rlg)+1
   ALLOCATE(tchrono(nvalmax))
   DO i=1,nvalmax 
     tchrono(i)=tdb+(i-1.0_rlg)*perchrono
   END DO
   indchrono=1
   l_chrono=.TRUE.

! meteorological parameters
!--------------------------

   l_meteoverif=.false.



! initializations
!----------------

   cflumaxmax=0.0_rsh
   cflvmaxmax=0.0_rsh

! Unset the generating potential because not ready yet
!-----------------------------------------------------

   IF (l_genpot) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'ERROR: l_genpot=.true.'
     WRITE(ierrorlog,*) 'Potential generator is not available yet.'
     WRITE(ierrorlog,*) 'Simulation stopped'
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   END IF

! verify CPP keys
!----------------

   IF ( .NOT.l_obc_mars .AND. .NOT.l_obc_cycl .AND.l_obc_char .AND.l_obc_diri ) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'ERROR: You must choose one of the following logicals '
     WRITE(ierrorlog,*) 'to define open boundary conditions schemes:'
     WRITE(ierrorlog,*) ' l_obc_cycl : apply cyclic open boundary conditions'
     WRITE(ierrorlog,*) ' l_obc_mars : apply MARS standard s open boundary conditions'
     WRITE(ierrorlog,*) ' l_obc_diri : apply Dirichlet s open boundary conditions'
     WRITE(ierrorlog,*) ' l_obc_char : apply characteristics method for open boundary conditions'
     WRITE(ierrorlog,*) 'Simulation stopped'
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   END IF

   IF ( .NOT. l_obc_cycl) THEN
   IF ( l_obc_cycl_x .OR. l_obc_cycl_y ) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'ERROR: You do not use cyclic open boundary conditions,'
     WRITE(ierrorlog,*) 'consequently l_obc_cycl_x and l_obc_cycl_y must be equal to .false.'
     WRITE(ierrorlog,*) 'Simulation stopped'
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   END IF
   END IF

   IF (l_obc_cycl) THEN
   IF ( .NOT.l_obc_cycl_x .AND. .NOT.l_obc_cycl_y ) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'ERROR: You use cyclic open boundary conditions,'
     WRITE(ierrorlog,*) 'consequently l_obc_cycl_x or l_obc_cycl_y must be set to .true.'
     WRITE(ierrorlog,*) 'Simulation stopped'
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   END IF
   END IF








! Check the variables in parameters.F90
! -------------------------------------
   IF (obc_width < 2) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'ERROR : Whatever your configuration,'
     WRITE(ierrorlog,*) 'you MUST have obc_width larger than 2 or equal to 2'
     WRITE(ierrorlog,*) 'The simulation is stopped'
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   ENDIF

   IF (valmanq < rg_valmanq_io-epsilon .OR. valmanq > rg_valmanq_io+epsilon) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'ERROR in parameters.F90 : valmanq=',valmanq,' and rg_valmanq_io=',rg_valmanq_io
     WRITE(ierrorlog,*) 'valmanq and rg_valmanq_io must be equal in value (but not for the type)'
     WRITE(ierrorlog,*) 'The simulation is stopped'
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   ENDIF


! Check the variables in MUSTANG_param
! --------------------------------
   IF (l_dzsmaxuni .AND. dzsmaxuni<hseduni .AND. ksdmax==1) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'ERROR in parasedim :  Only one layer allowed (ksdmax=1) and '
     WRITE(ierrorlog,*) 'the maximum sediment layer thickness dzsmaxuni :',dzsmaxuni, ' m'
     WRITE(ierrorlog,*) 'is smaller than the initial sediment thickness hseduni : ',hseduni,' m'
     WRITE(ierrorlog,*) 'The simulation is stopped'
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   ENDIF

   IF (MASTER) THEN
     IF (tstart_sed<tdebsubs .AND. l_dyn_insed) THEN
          WRITE(iwarnlog,*) ' '
          WRITE(iwarnlog,*) 'WARNING : Variable date_start_sed is uncorrect. '
          WRITE(iwarnlog,*) 'WARNING : consolidation/diffusion simulations will start '
          WRITE(iwarnlog,*) 'WARNING : at the begining of substances simulations :' , date_startsub
     ENDIF
   ENDIF
   tstart_sed=MAX(tstart_sed,tdebsubs)
   t_sed=tstart_sed
!!!


! formats of reading
!-------------------

 3001 FORMAT(a4,i1,6f12.7,5x,2f10.2,6i4,1x,a4)

   !PRINT*, 'exit INIT routine'

  END SUBROUTINE init
