 MODULE subreaddat
!!======================================================================
!!                   ***  MODULE  subreaddat  ***
!! Ocean dynamics Bio Sed:  Read files *.dat (rivers, discharges, bio variables...)
!&E
!&E ** History :
!&E       !  2006-09  (V. Garnier)  Original code
!&E       !  2008-11  (V. Garnier)  Remove var_read_obc
!&E       !  2010-06  (V.M. Faure, N. Bousserez, M. Baklouti) coupling with eco3m
!&E
!!======================================================================





























































!! * Modules used
   USE parameters,ONLY : rsh,lchain,imin,imax,jmin,jmax,kmax,limin,limax,ljmin,ljmax
   USE comvars2d, ONLY : iscreenlog,hm
   USE comvars3d, ONLY : dsigu


   USE toolmpi, ONLY : MASTER

   IMPLICIT NONE
   PRIVATE

!! * Accessibility
   PUBLIC sub_read_var,sub_read_vardiag           ! routine called by main.F90

!! * Shared module variables

!! * Private variables

 CONTAINS

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

  SUBROUTINE sub_read_var



!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sub_read_var  ***
!&E
!&E ** Purpose : read "variables.dat" file describing substances
!&E              other than temperature and salinity
!&E
!&E ** Description : usefull for dissolved, biologic, sedimentologic... variables
!&E
!&E ** Called by : main.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables : name_var, long_name_var, standard_name_var, unit_var,
!&E                         valid_min_var, valid_max_var, cini_wat,cini_sed,t90,f90,
!&E                         nv_grav,nv_sand,nv_mud,nv_cp,nv_ncp,nv_sorb,nv_dis,nv_fix,nv_dri,nv_int,
!&E                         itypv_r,itypv,nv_adv,nv_state,nv_tot,typart,sub_flx_atm,cv_rain,typdiss
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2006-08-23 (B. Thouvenin) Original code
!&E       !  2007-07-25 (M. Sourisseau) introduction nom des substances
!&E                                     lues dans fichiers OBC et INIT
!&E       !  2008-01    (M. Sourisseau) introduction of MES (module of biology)
!&E       !  2008-03    (C. Struski) introduction of parameter l_out_subs
!&E                                  to save the substance in output files or not
!&E       !  2008-09    (G. Galibert/ALTRAN) attributes of output files that respect
!&E                                          NetCDF OCO v1.1 file format
!&E                                          reads long_name, standard_name, valid_min and
!&E                                          valid_max attributes from file_tracer
!&E       !  2009-02    (P. Le Hir) new sediment characteristics for "new" sand/mud model
!&E       !  2011-02    (B. Thouvenin) correction for key_subs_part_eqsubdt key
!&E       !  2014-05    (B. Thouvenin) modification of variable.dat reading (settling velocities, options for sand)
!&E       !  2016-12    (B.Thouvenin) version sedim_MUSTANG : only mixsed (mudsand not operationnal)
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters
   USE comvars2d
   USE comvars3d
   USE comsubstance
   USE comMUSTANG



!! * Arguments

!! * Local declarations
   LOGICAL                                   :: ex
   LOGICAL                                   :: l_varassoc
   INTEGER                                   :: ivpc

   INTEGER                                   :: isubs,eof,iv,isubs_r,nfrac,ind_white, &
                                                iiv,indgrp,k,IERR_MPI,ifois_ndt
   CHARACTER(LEN=4)                          :: typvar
   CHARACTER(LEN=lchain)                     :: nomvar_r,unitvar_r,  &
                                                long_name_var_r,standard_name_var_r
   CHARACTER(LEN=5)                          :: comment
   LOGICAL, DIMENSION(nb_var)                :: l_out_subs_r 

   LOGICAL, DIMENSION(nb_var)                :: l_outsandrouse_r,l_sand2D_r

   REAL(KIND=rsh), DIMENSION(4,nb_var)       :: ws_free_para_r
   REAL(KIND=rsh), DIMENSION(nb_var)         :: ws_free_min_r,ws_free_max_r
   REAL(KIND=rsh), DIMENSION(2,nb_var)       :: ws_hind_para_r 
   INTEGER       , DIMENSION(nb_var)         :: ws_hind_opt_r,ws_free_opt_r 
   REAL(KIND=rsh), DIMENSION(nb_var)         :: tocd_r,t90_r,flx_atm_r
   REAL(KIND=rsh), DIMENSION(nb_var)         :: cini_wat_r,cv_rain_r,cini_air_r
   CHARACTER(LEN=lchain), DIMENSION(nb_var)  :: obc_cv_name_r,init_cv_name_r

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

! save into simu.log
!-------------------
   IF (MASTER) THEN
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) '**************************************************'
     WRITE(iscreenlog,*) '**************** SUBREADVAR.F90 **********************'
     WRITE(iscreenlog,*) '**************************************************'
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) 'file defining simulated substances (other than temperature and salinity) :'
     WRITE(iscreenlog,*) TRIM(name_filesubs)
   ENDIF

! Initialize number of variables depending on their type
! (gravel,sand,mud,no-constitutive,constitutive,dissolved,fixed,driving or intermediate)
   nv_grav=0
   nv_sand=0
   nv_mud=0
   nv_cp=0
   nv_ncp=0
   nv_sorb=0
   nv_dis=0
   nv_fix=0
   nv_dri=0
   nv_int=0
   isubs=0


   name_var(-1)='Temperature'
   name_var(0)='Salinite'
   long_name_var(-1)='temperature'
   long_name_var(0)='salinity'
   standard_name_var(-1)='sea_water_potential_temperature'
   standard_name_var(0)='sea_water_salinity'
   unit_var(-1)='degrees_Celsius'
   unit_var(0)='1e-3'
   valid_min_var(-1)=0.0_rsh
   valid_min_var(0)=0.0_rsh
   valid_max_var(-1)=40.0_rsh
   valid_max_var(0)=40.0_rsh

! initialize concentrations
   cini_wat_r(:)=0.0_rsh
!  cini_sed_r(:)=0.0_rsh
   cini_air_r(:)=0.0_rsh
   flx_atm_r(:)=0.0_rsh
   cv_rain_r(:)=0.0_rsh
   unit_modif_mudbio_N2dw(:)=1.0_rsh
   eof=0
   INQUIRE(file=name_filesubs,exist=ex)
   IF (ex) THEN
     OPEN(49,file = name_filesubs,form='formatted')
     comment='debut'
     DO WHILE (comment /= '*****')
       READ(49,'(a)',iostat=eof) comment
     END DO
     DO WHILE (eof==0)
       READ(49,*,iostat=eof)isubs_r
       isubs=isubs+1
       IF (eof==0) THEN
         READ(49,'(a)',iostat=eof) nomvar_r
         ind_white=INDEX(nomvar_r,' ')
         name_var(isubs)=TRIM(ADJUSTL(ADJUSTR(nomvar_r(1:ind_white))))
         READ(49,'(a)',iostat=eof) long_name_var_r
         ind_white=INDEX(long_name_var_r,' ')
         long_name_var(isubs)=TRIM(ADJUSTL(ADJUSTR(long_name_var_r(1:ind_white))))
         READ(49,'(a)',iostat=eof) standard_name_var_r
         ind_white=INDEX(standard_name_var_r,' ')
         standard_name_var(isubs)=TRIM(ADJUSTL(ADJUSTR(standard_name_var_r(1:ind_white))))
         READ(49,'(a)',iostat=eof) unitvar_r 
         ind_white=INDEX(unitvar_r,' ')
         unit_var(isubs)=TRIM(ADJUSTL(ADJUSTR(unitvar_r(1:ind_white))))
         READ(49,*,iostat=eof) valid_min_var(isubs)
         READ(49,*,iostat=eof) valid_max_var(isubs)
         READ(49,'(a)',iostat=eof) typvar
! pour l instant, on se passe du code, on tente le traitement avec les noms seulement
!    READ(49,*,iostat=eof) icodv_r(isubs)
!     WRITE(12345,*)'son code est = ',icodv_r(isubs)
         IF (typvar=='GRAV') THEN
           itypv_r(isubs)=1
           nv_grav=nv_grav+1
         ELSE IF (typvar=='SAND') THEN
           itypv_r(isubs)=2
           nv_sand=nv_sand+1
         ELSE IF (typvar=='MUDS' .OR. typvar=='MUDB') THEN
           itypv_r(isubs)=3
           nv_mud=nv_mud+1
         ELSE IF (typvar=='PART') THEN
           itypv_r(isubs)=3
           nv_cp=nv_cp+1
         ELSE IF (typvar=='NoCP') THEN
           itypv_r(isubs)=4
           nv_ncp=nv_ncp+1
         ELSE IF (typvar=='SORB') THEN
           itypv_r(isubs)=5
           nv_sorb=nv_sorb+1
         ELSE IF (typvar=='DISS') THEN
           itypv_r(isubs)=6
           nv_dis=nv_dis+1
         ELSE IF (typvar=='FIXE') THEN
           itypv_r(isubs)=7
           nv_fix=nv_fix+1
         ELSE IF (typvar=='DRIV') THEN
           itypv_r(isubs)=8
           nv_dri=nv_dri+1
         ELSE IF (typvar=='INTE') THEN
           itypv_r(isubs)=9
           nv_int=nv_int+1
         END IF
         

! modification of unit if typvar=MUDB (mudbio)
! biological variable which is usually expressed in micromole/L but which must be constitutive particulate variable
!  if 1 because of its own large settling velocities
         IF (typvar=='MUDB') THEN
            unit_modif_mudbio_N2dw(isubs)=12.0e-04_rsh
         ENDIF
         
         READ(49,'(a)',iostat=eof)

! only for variable SAND
         IF (itypv_r(isubs) ==2) THEN

           READ(49,*,iostat=eof) l_sand2D_r(isubs)
           READ(49,*,iostat=eof) l_outsandrouse_r(isubs)

         ENDIF
         READ(49,'(a)',iostat=eof)
         IF (itypv_r(isubs)==3 .OR. itypv_r(isubs)==4 ) THEN
!   reading settling velocities ONLY for MUDS, PART and NoCP
           READ(49,*,iostat=eof) ws_free_opt_r(isubs)  ! option choosed for formulations of free settling velocity
           READ(49,*,iostat=eof) ws_free_min_r(isubs),ws_free_max_r(isubs),ws_free_para_r(1:4,isubs)  ! parameters for free settling velocity
           READ(49,*,iostat=eof) ws_hind_opt_r(isubs)  ! option choosed for formulations of hindered settling velocity
           READ(49,*,iostat=eof) ws_hind_para_r(1:2,isubs) ! parameters for hindered settling velocity
         ELSE IF (itypv_r(isubs) == 5) THEN
!   reading only the associated constitutive variable for SORB variable
           READ(49,*,iostat=eof) nomvar_r
           ind_white=INDEX(nomvar_r,' ')
           name_varpc_assoc(isubs)=TRIM(ADJUSTL(ADJUSTR(nomvar_r(1:ind_white))))
         END IF 
         IF (itypv_r(isubs) < 5) THEN
!   reading diameter of particles, critical stress of deposition and density of particle for all particles EXCEPT SORB
           READ(49,*,iostat=eof) diam_r(isubs)      !diametre des grains
           READ(49,*,iostat=eof) tocd_r(isubs)      !tension critique de depot
           READ(49,*,iostat=eof) ros_r(isubs)       !densite des grains
         ENDIF

         READ(49,'(a)',iostat=eof)
         IF (itypv_r(isubs) < 8) THEN
           READ(49,*,iostat=eof) t90_r(isubs)   ! utilise seulement pour usage de T90
           READ(49,*,iostat=eof) flx_atm_r(isubs)
           READ(49,*,iostat=eof) cv_rain_r(isubs)
           READ(49,*,iostat=eof) cini_wat_r(isubs)
           READ(49,*,iostat=eof) cini_sed_r(isubs)
           READ(49,*,iostat=eof) cini_air_r(isubs)
         END IF
         READ(49,*,iostat=eof) l_out_subs_r(isubs)
         READ(49,'(a)',iostat=eof) nomvar_r
         ind_white=INDEX(nomvar_r,' ')
         init_cv_name_r(isubs)=TRIM(ADJUSTL(ADJUSTR(nomvar_r(1:ind_white))))
         READ(49,'(a)',iostat=eof) nomvar_r
         ind_white=INDEX(nomvar_r,' ')
         obc_cv_name_r(isubs)=TRIM(ADJUSTL(ADJUSTR(nomvar_r(1:ind_white))))
         READ(49,'(a)',err=11,END=11)
       END IF
     END DO
11 CONTINUE
     CLOSE(49)
   END IF     ! find du test sur l existence du fichier

! initialize the number of variables according to their type
   nvpc=nv_cp+nv_mud+nv_sand+nv_grav
   nvp=nvpc+nv_ncp+nv_sorb
   nv_adv=nvp+nv_dis
   nv_state=nv_adv+nv_fix
   nv_tot=nv_state+nv_dri+nv_int





   IF (nv_tot /= nb_var) THEN
     WRITE(ierrorlog,*)'WARNING - the total number of variables read from the file'
     WRITE(ierrorlog,*)TRIM(name_filesubs),' is DIFFERENT from the NB_VAR parameter'
     WRITE(ierrorlog,*)'in parameter.F90 routine '
     WRITE(ierrorlog,*)'nb_var in parameter.F90= ',nb_var
     WRITE(ierrorlog,*)'nv_tot read in variable.dat ',nv_tot

     WRITE(ierrorlog,*)' - Make sure there is no mistake in ',TRIM(name_filesubs),'file !'

     WRITE(ierrorlog,*)'The simulation is stopped'
     CALL MPI_FINALIZE(IERR_MPI)
     STOP
   END IF

! allocation des tableaux se rapportant aux variables (declarees dans comsubstance)
! et dont les dimensions dependent des nombres qui viennent d etre lus)
! ---------------------------------------------------------------------------------
   ALLOCATE(cv_wat(nv_tot,kmax,liminm1:limaxp2,ljminm1:ljmaxp2))
! intermediate variables (macho scheme)
   ALLOCATE(c2(nv_tot,kmax,liminm1:limaxp2,ljminm1:ljmaxp2))
   c2(:,:,:,:)=0.0_rsh


! sink-source for substances
   ALLOCATE(phicon(nv_tot,kmax,limin:limax,ljmin:ljmax))
   ALLOCATE(phicon_atm(nv_tot,limin:limax,ljmin:ljmax))
   ALLOCATE(phicon_drycell(nv_tot,kmax,limin:limax,ljmin:ljmax))
   phicon_drycell(:,:,:,:)=0.0_rsh
   phicon(:,:,:,:)=0.0_rsh
   phicon_atm(:,:,:)=0.0_rsh

   ALLOCATE(obc_cv_name(nv_adv))
   ALLOCATE(init_cv_name(nv_tot))
   ALLOCATE(l_out_subs(nv_state+nv_dri))
   ALLOCATE(t90(nv_state))
! plus correct mais genant pour marsWEB   ALLOCATE(f90(nv_state))
   ALLOCATE(f90(nv_tot))
   f90(1:nv_tot) = 1.0_rsh
   ALLOCATE(sub_flx_atm(nv_adv))
   ALLOCATE(cv_rain(nv_adv))
   ALLOCATE(cini_wat(nv_state))
!   ALLOCATE(cini_sed(nv_state))  (P. Le Hir : now in module sedimento)
   ALLOCATE(cini_air(nv_state))
   ALLOCATE(typart(nv_state))
   ALLOCATE(typdiss(nv_adv))
   ALLOCATE(itypv(nv_tot))
   ALLOCATE(obc_cv_wat_e(nv_adv,kmax,imax-obc_width+1:imax,ljmin:ljmax))
   ALLOCATE(obc_cv_wat_w(nv_adv,kmax,imin:imin+obc_width,ljmin:ljmax))
   ALLOCATE(obc_cv_wat_n(nv_adv,kmax,limin:limax,jmax-obc_width+1:jmax))
   ALLOCATE(obc_cv_wat_s(nv_adv,kmax,limin:limax,jmin:jmin+obc_width))
   ALLOCATE(irk_mod(nv_tot))
   ALLOCATE(irk_fil(nv_tot))
   ALLOCATE(ws3(kmax,nv_adv,limin:limax,ljmin:ljmax))
   ws3(:,:,:,:)=0.0_rsh
   IF (nvp > 0) THEN

      ALLOCATE(ws_free_opt(nvp)) 
      ALLOCATE(ws_hind_opt(nvp)) 
      ALLOCATE(ws_free_para(4,nvp)) 
      ALLOCATE(ws_free_min(nvp)) 
      ALLOCATE(ws_free_max(nvp)) 
      ALLOCATE(ws_hind_para(2,nvp)) 
      ALLOCATE(tocd(nvp))
      ws_free_min(:)=0.0_rsh
      ws_free_max(:)=0.0_rsh
      ws_free_para(:,:)=0.0_rsh
      ws_hind_para(:,:)=0.0_rsh
      ws_free_opt(:)=0
      ws_hind_opt(:)=0
      tocd(:)=0.0_rsh
   ENDIF
   CALL alloc_fxfyfz
   ALLOCATE(l_subs2D(nv_state))
   l_subs2D(:)=.false.

   ALLOCATE(l_outsandrouse(nvp))
   l_outsandrouse(:)=.false.








! appel au rangement des variables dans l ordre beni du modele
! conversion de l ordre lu dans l ordre modele : irk_mod(ilu)
! conversion de l ordre modele dans l ordre lu : irk_fil(imod)
!  Exemple : lecture des variables dans fichier dans l odre : MOP,Vas1,O2D,Sabl
!            rangement dans l ordre du modele : Sabl,Vas1,MOP,O2D
!            on obtient : irk_mod(ilu=1)=3 ; irk_mod(2)=2 ; irk_mod(3)=4 ; irk_mod(4)=1
!            et en sens inverse : irk_fil(ivmod=1)=4;irk_fil(2)=2;irk_fil(3)=1;irk_fil(4)=3
! -----------------------------------------------------------------------------------------
   CALL suborder



! reperage de variables particulaires constitutives associees aux variables part. SORB
! ------------------------------------------------------------------------------------
   ALLOCATE(irkm_var_assoc(nvp))
   irkm_var_assoc(:)=0
   DO iv=1,nv_sorb
     isubs=nvpc+nv_ncp+iv
     irkm_var_assoc(isubs)=0
     l_varassoc=.FALSE.
     DO ivpc=1,nvpc
       IF ( TRIM(ADJUSTL(ADJUSTR(name_var(irk_fil(ivpc))))) == TRIM(ADJUSTL(ADJUSTR(name_varpc_assoc(irk_fil(isubs))))) )THEN
         irkm_var_assoc(isubs)=ivpc
         IF (MASTER) THEN
            WRITE(iscreenlog,*) ' '
            WRITE(iscreenlog,*)'constitutive part. variable where is sorbed'
            WRITE(iscreenlog,*)'the variable :',name_var(irk_fil(isubs))
            WRITE(iscreenlog,*)'is the variable ', name_var(irk_fil(ivpc))
         ENDIF
         l_varassoc=.TRUE.
         ws_free_min_r(irk_fil(isubs))=ws_free_min_r(irk_fil(ivpc))
         ws_free_max_r(irk_fil(isubs))=ws_free_max_r(irk_fil(ivpc))
         ws_free_para_r(:,irk_fil(isubs))=ws_free_para_r(:,irk_fil(ivpc))
         ws_hind_para_r(:,irk_fil(isubs))=ws_hind_para_r(:,irk_fil(ivpc))
         ws_free_opt_r(irk_fil(isubs))=ws_free_opt_r(irk_fil(ivpc))
         ws_hind_opt_r(irk_fil(isubs))=ws_hind_opt_r(irk_fil(ivpc))
         tocd_r(irk_fil(isubs))=tocd_r(irk_fil(ivpc))
         ros_r(irk_fil(isubs))=ros_r(irk_fil(ivpc))
       END IF
     END DO
     IF (.NOT.l_varassoc) THEN
       IF (MASTER) THEN
          WRITE(iscreenlog,*)' '
          WRITE(iscreenlog,*)'the SORB variable :',name_var(irk_fil(isubs))
          WRITE(iscreenlog,*)'does not have associated constitutive part. variable'
          WRITE(iscreenlog,*)'See variable.dat to give exactly the name of the constitutive associated variable'
          WRITE(iscreenlog,*)'otherwise, it is not a SORB variable, but a NoCP variable '
       ENDIF
       WRITE(ierrorlog,*)' '
       WRITE(ierrorlog,*)'the SORB variable :',name_var(irk_fil(isubs))
       WRITE(ierrorlog,*)'does not have associated constitutive part. variable'
       WRITE(ierrorlog,*)'See variable.dat to give exactly the name of the constitutive associated variable'
       WRITE(ierrorlog,*)'otherwise, it is not a SORB variable, but a NoCP variable '
       CALL MPI_FINALIZE(IERR_MPI)
       STOP
     END IF
   END DO

! save into simu.log
!-------------------
   IF (MASTER) THEN
    WRITE(iscreenlog,*) ' '
    WRITE(iscreenlog,*) 'SUBSTANCE NUMBER        NAME             UNIT            TYPE       RANK IN MODEL'
    DO isubs=1,nv_tot
!     WRITE(iscreenlog,*)  &
     WRITE(iscreenlog,'(4x,i4,6x,a18,2x,a18,5x,i4,9x,i4)')  &
      isubs,TRIM(name_var(isubs)),TRIM(unit_var(isubs)),itypv_r(isubs),irk_mod(isubs)
    END DO
    DO isubs=1,nv_tot
     WRITE(iscreenlog,*)' '
     WRITE(iscreenlog,*)'VARIABLE NAME : ',TRIM(name_var(isubs))
     IF (itypv_r(isubs)==3 .OR. itypv_r(isubs)==4) THEN
       
       WRITE(iscreenlog,*)'Free settling velocity method : ',ws_free_opt_r(isubs)
       WRITE(iscreenlog,*)'Free settling velocity MIN and MAX : ',ws_free_min_r(isubs),ws_free_max_r(isubs)
       WRITE(iscreenlog,*)'Free settling velocity parameters : ',ws_free_para_r(1:4,isubs)
       WRITE(iscreenlog,*)'Hindered settling velocity method : ',ws_hind_opt_r(isubs)
       WRITE(iscreenlog,*)'Hindered settling velocity parameters : ',ws_hind_para_r(1:2,isubs)
     ELSE IF (itypv_r(isubs) == 5) THEN
       WRITE(iscreenlog,*)'Particulate Constitutive associated Variable  : ',name_varpc_assoc(isubs)
     ENDIF
     IF (itypv_r(isubs)< 6) THEN
       WRITE(iscreenlog,*)'critical shear stress for deposit : ',tocd_r(isubs)
       WRITE(iscreenlog,*)'grain density                     : ',ros_r(isubs)
       WRITE(iscreenlog,*)'grain diameter                    : ',diam_r(isubs)
     END IF
     IF (itypv_r(isubs) < 8) THEN
       WRITE(iscreenlog,*)'T90 (inutile si cycle bio complet ou si key_contaminant) : ',t90_r(isubs)
       WRITE(iscreenlog,*)'depot atmospherique (masse/m2/seconde) : ',flx_atm_r(isubs)
       WRITE(iscreenlog,*)'concentration in rain water            : ',cv_rain_r(isubs)
       WRITE(iscreenlog,*)'uniform initial conc. in water column  : ',cini_wat_r(isubs)
       WRITE(iscreenlog,*)'uniform initial conc. in sediment      : ',cini_sed_r(isubs)
       WRITE(iscreenlog,*)'uniform initial conc. in air           : ',cini_air_r(isubs)
     END IF
    END DO

    WRITE(iscreenlog,*)' '
    WRITE(iscreenlog,*)'STOCK OF VARIABLES'
    WRITE(iscreenlog,*)'number of GRAV                        : ',nv_grav
    WRITE(iscreenlog,*)'number of SAND                        : ',nv_sand
    WRITE(iscreenlog,*)'number of MUDS                        : ',nv_mud
    WRITE(iscreenlog,*)'number of PART (other constitutive)   : ',nv_cp
    WRITE(iscreenlog,*)'number of part. var. constitutive     : ',nvpc
    WRITE(iscreenlog,*)'number of part. var. NO constitutive  : ',nv_ncp
    WRITE(iscreenlog,*)'number of part. var. SORB             : ',nv_sorb
    WRITE(iscreenlog,*)'number of dissolved var. (DISS)       : ',nv_dis
    WRITE(iscreenlog,*)'number of FIXE                        : ',nv_fix
    WRITE(iscreenlog,*)'number of DRIV (driving var. )        : ',nv_dri
    WRITE(iscreenlog,*)'number of INTE (intermediate var. )   : ',nv_int
    WRITE(iscreenlog,*)
    WRITE(iscreenlog,*) 'number of state variables            : ',nv_state, '+ salinite et temperature'
    WRITE(iscreenlog,*) 'total number of variables (+fix +driving +inter) : ',nv_tot
   ENDIF

! ---------------------------------------------------------
! SEDIMENTOLOGY
! ---------------------------------------------------------

! identification of igrav1,igrav2,isand1,isand2,imud1,imud2
! ---------------------------------------------------------
   igrav1=1
   igrav2=nv_grav
   isand1=igrav2+1
   isand2=igrav2+nv_sand
   imud1=isand2+1
   imud2=isand2+nv_mud

! transfert des tableaux dans le bon ordre pour les donnees sedimento
! et les reactions autres que celles de key_biolo
! -------------------------------------------------------------------

   DO iv=imud1,nvp
     ws_free_opt(iv)=ws_free_opt_r(irk_fil(iv))
     ws_free_min(iv)=ws_free_min_r(irk_fil(iv))
     ws_free_max(iv)=ws_free_max_r(irk_fil(iv))
     ws_free_para(:,iv)=ws_free_para_r(:,irk_fil(iv))
     ws_hind_opt(iv)=ws_hind_opt_r(irk_fil(iv))
     ws_hind_para(:,iv)=ws_hind_para_r(:,irk_fil(iv))
     tocd(iv)=tocd_r(irk_fil(iv))
     ws3(:,iv,:,:)=ws_free_max(iv)
   END DO

   DO iv=1,igrav2
     l_subs2D(iv)=.TRUE.
   ENDDO
   DO iv=isand1,isand2
     l_subs2D(iv)=l_sand2D_r(irk_fil(iv))
     l_outsandrouse(iv)=l_outsandrouse_r(irk_fil(iv))
   ENDDO

   DO iv=1,nv_tot
     itypv(iv)=itypv_r(irk_fil(iv))
   END DO
   DO iv=1,nv_state
     t90(iv)=t90_r(irk_fil(iv))
     cini_wat(iv)=cini_wat_r(irk_fil(iv))*unit_modif_mudbio_N2dw(irk_fil(iv))
     cini_air(iv)=cini_air_r(irk_fil(iv))*unit_modif_mudbio_N2dw(irk_fil(iv))
     init_cv_name(iv)=init_cv_name_r(irk_fil(iv))
     l_out_subs(iv)=l_out_subs_r(irk_fil(iv))


   IF (t90(iv) /= 0.0_rsh) THEN
      WRITE(ierrorlog,*) 'Are you sure you want to use t90 for the ',iv,' th substance ?'
      WRITE(ierrorlog,*) 'Since you are supposed to use conservative substance,'
      WRITE(ierrorlog,*) 'the simulation is stopped'
      CALL MPI_FINALIZE(IERR_MPI)
      STOP
   END IF


   END DO

   DO iv=nv_state+1,nv_state+nv_dri
     l_out_subs(iv)=l_out_subs_r(irk_fil(iv))
   END DO

   DO iv=1,nv_adv
     obc_cv_name(iv)=obc_cv_name_r(irk_fil(iv))
     cv_rain(iv)=cv_rain_r(irk_fil(iv))*unit_modif_mudbio_N2dw(irk_fil(iv))
     sub_flx_atm(iv)=flx_atm_r(irk_fil(iv))*unit_modif_mudbio_N2dw(irk_fil(iv))


   END DO

! Mise en groupe des variables particulaires
   IF (nvp>0) THEN
! Allocation des tableaux se rapportant aux variables non dissoutes
     ALLOCATE(iv_ws_max(nvp,nvp)) 
     ALLOCATE(nv_ws_max(nvp)) 
     ALLOCATE(ndt_part(nvp)) 
     ALLOCATE(sv_part(nvp))
     ALLOCATE(ap_ws_gr(nvp))

! Initialisation des différentes variables et tableaux
     nb_ws_max=0
     iv_ws_max(:,:)=0
     nv_ws_max(:)=0
     ap_ws_gr(:)=0
     ndt_part(:)=0
     sv_part(:)=0.0_rsh


! Equation of substances with settling velocity
! is solved with a sub time step

! estimate of the minimum layer thickness
! ---------------------------------------
     dhmin=hm

     IF(.NOT.l_modele2d) THEN
       DO k=1,kmax
         dhmin=MIN(dhmin,hm*dsigu(k))
       ENDDO
     ENDIF




     nb_ws_max=0
     DO iv=1,nvp
       nb_ws_max=iv     ! number of variables with different settling velocities
! all of them for 1 because ws3 is a 4D variable
       nv_ws_max(iv)=1  ! each particulate variable is relative to a unique group
! because ws3 is a 4D variable
       sv_part(iv)=ws_free_max(iv)   ! important to determinate the number of iterations
! index of variable associated with a settling velocity
       iv_ws_max(iv,1)=iv
       IF(iv < imud1 ) THEN
           ndt_part(nb_ws_max)=1  ! will be calculated after estimation of settling velocity
       ELSE
         IF (ws_free_max(iv) == 0.0_rsh) THEN
           WRITE(ierrorlog,*) ''
           WRITE(ierrorlog,*) ' ******* Error detected from SUB_READ_VAR from SUBREADDAT.F90 ******* '
           WRITE(ierrorlog,*) 'the variable :',name_var(irk_fil(iv))
           WRITE(ierrorlog,*) 'is a particulate one'
           WRITE(ierrorlog,*) 'but its settling velocity is set to 0.0 in ',TRIM(name_filesubs), ' file.'
           WRITE(ierrorlog,*) 'Consequently, the sub time step used to solve substance equation cannot be estimated'
           WRITE(ierrorlog,*) 'The simulation is stopped'
           CALL MPI_FINALIZE(IERR_MPI)
           STOP
         END IF
         IF(l_modele2d) THEN
           ndt_part(nb_ws_max)=1
         ELSE
! ajuster si besoin en augmentant ifois_ndt mais attention, c est cher en temps calcul
           ifois_ndt=1
           ndt_part(nb_ws_max)=ifois_ndt*NINT(dt/MIN((dhmin/ws_free_max(iv)),dt))
         ENDIF
       END IF
     ENDDO





   ENDIF

!!!!#ifdef key_subs_part_eqsubdt
! les flux horizontaux sont recalcules de toute facon a chaque sous pas de temps
   nvpi=nvp
!!!!#else
!!!!   nvpi=0
!!!!#endif

   IF (MASTER) THEN

     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) 'SUBSTANCES WITH A VARIABLE SETTLING VELOCITY '

     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*)'NUMBER GROUP OF SUBSTANCES            : ',nb_ws_max
     DO iiv=1,nb_ws_max
      WRITE(iscreenlog,*) ' '
      WRITE(iscreenlog,*)'group number                   : ',iiv
      WRITE(iscreenlog,*)'number of substances           : ',nv_ws_max(iiv)
      WRITE(iscreenlog,*)'group settling velocity        : ',sv_part(iiv)
      WRITE(iscreenlog,*)'group number of iterations     : ',ndt_part(iiv)
      DO iv=1,nv_ws_max(iiv)
       WRITE(iscreenlog,*) 'substance number        name' 
       WRITE(iscreenlog,'(4x,i4,6x,a18,2x)')  &
       iv,TRIM(name_var(irk_fil(iv_ws_max(iiv,iv))))
      ENDDO      
     ENDDO      
     IF (nvp/=0 .AND. nvpi==0) THEN
!       WRITE(iscreenlog,*) 'The horizontal fluxes of particulate variables'
!       WRITE(iscreenlog,*) 'are estimated at each hydrodynamics time step'
!       WRITE(iscreenlog,*) 'and remain constant.'
!       WRITE(iscreenlog,*) 'The sub time step (if any) is only used to solve'
!       WRITE(iscreenlog,*) 'the vertical equation of the particulate variables'
!       WRITE(iscreenlog,*) '(pas fractionnaire)'
       WRITE(iscreenlog,*) '  **** SUBREADDAT.F90 (sub_read_var) ****'
       WRITE(ierrorlog,*) 'You should not have nvpi=0'
       WRITE(ierrorlog,*) 'The simulation is stopped'
       WRITE(ierrorlog,*) ''
         CALL MPI_FINALIZE(IERR_MPI)
         STOP
     ELSE IF (nvp /=0) THEN
       WRITE(iscreenlog,*) 'The complete subtances equation is solved from a sub time step : '
       WRITE(iscreenlog,*) 'The horizontal fluxes of particulate variables'
       WRITE(iscreenlog,*) 'are estimated at each sub time step'
     END IF
     WRITE(iscreenlog,*)' nvpi=',nvpi
   ENDIF
   typart(1:nv_state)=0.
   typart(1:nvpc)=1.
   typdiss(nvp+1:nv_adv)=1.
   typdiss(1:nvp)=0.

!  initialization  of l_sflxsubatm
   l_subflxatm=.false.
   IF(l_subflxatm_xyt) THEN
     l_subflxatm=.true.
   ELSE
     DO iv=1,nv_adv
       IF (cv_rain(iv) /= 0.0_rsh .OR. sub_flx_atm(iv) /= 0.0_rsh) THEN
         l_subflxatm=.true.
       ENDIF
     END DO
   ENDIF



  END SUBROUTINE sub_read_var

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

  SUBROUTINE sub_read_vardiag


  WRITE(*,*) 'sub_read_vardiag : You should not have seen this print! error?'


  END SUBROUTINE sub_read_vardiag


 END MODULE

