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

#include "toolcpp.h"

   !! * 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
#ifdef key_biolo
   USE bioloinit, ONLY : ndiag_1d,ndiag_2d,ndiag_3d,ndiag_3d_wat,ndiag_3d_sed, &
                         diag_1d,diag_2d,diag_3d_wat,                          &
                         ndiag_tot,idimv_r,l_diag_sortie,name_vardiag,         &
                         long_name_vardiag,standard_name_vardiag,unit_vardiag, &
                         valid_min_vardiag,valid_max_vardiag,irk_diag,         &
                         name_vardiag,ihour_previous
   USE bioloinit, ONLY : biolo_init_id,biolo_init_iv
#endif
#if defined key_siggen || defined key_gencoord
   USE comsiggen, ONLY : hc,dcusds
#endif
   USE_MPI 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

#ifdef key_substance

   !&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
#ifdef key_contaminant
   USE contareaddat 
#endif
#if defined key_biolo && (defined key_N_tracer || defined key_P_tracer)
   USE bioloinit, ONLY : biolo_create_vartracer
#endif

   !! * Arguments

   !! * Local declarations
   LOGICAL                                   :: ex
   LOGICAL                                   :: l_varassoc
   INTEGER                                   :: ivpc
#ifdef key_contaminant
   INTEGER                                   :: ifracv
   INTEGER,DIMENSION(nb_var)                 :: nb_frac_var_r,icon_var_r
   INTEGER,DIMENSION(5)                      :: ivdl_conta_tmp 
#if defined key_species
   REAL(KIND=rsh),DIMENSION(nb_var,nfracvmax):: frac_cwat_r,frac_csed_r
   INTEGER,DIMENSION(nb_var,nfracvmax)       :: ispc_var_r
#else
   INTEGER,DIMENSION(nb_var)                 :: ispc_var_r
#endif
#endif
   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 
#ifdef key_sand2D
   LOGICAL, DIMENSION(nb_var)                :: l_outsandrouse_r,l_sand2D_r
#endif
   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_MPI (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_MPI

   ! 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
#ifdef key_contaminant
  l_conta=.false.
  nb_conta=0
#endif

   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
         
#ifdef key_contaminant
   ! typvar='CONTA' : new and specific intermediate type in order 
   ! to introduce chemical contaminant species  
         IF (typvar=='CONTA') THEN
            nb_conta=nb_conta+1
            l_conta=.true.
            itypv_r(isubs)=6
            nv_dis=nv_dis+1
            IF (nb_conta > 5) THEN
             IF_MPI (MASTER) THEN
               WRITE(iscreenlog,*)'number of contaminants in ',trim(name_filesubs)
               WRITE(iscreenlog,*)'is greater than 5 '
               WRITE(iscreenlog,*)'MARS/MET&OR is limited to 5 contaminants'
               CALL_MPI MPI_FINALIZE(IERR_MPI)
               STOP
             ENDIF_MPI
            END IF
!           first species = free (or true) dissolved  species
            ivdl_conta_tmp(nb_conta)=isubs
         END IF
#endif
    ! modification of unit if typvar=MUDB (mudbio)
    ! biological variable which is usually expressed in micromole/L but which must be constitutive particulate variable
    !  if key_sedim_MUSTANG 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
#ifdef key_sand2D
           READ(49,*,iostat=eof) l_sand2D_r(isubs)
           READ(49,*,iostat=eof) l_outsandrouse_r(isubs)
#else
           READ(49,'(a)',iostat=eof)
           READ(49,'(a)',iostat=eof)
#endif
         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 defined key_biolo && (defined key_N_tracer || defined key_P_tracer)
   CALL biolo_create_vartracer(tocd_r,t90_r,flx_atm_r,cv_rain_r,cini_wat_r,     &
                          cini_air_r,l_out_subs_r,init_cv_name_r,obc_cv_name_r, &
                          ws_free_opt_r,ws_free_min_r,ws_free_max_r,            &
                          ws_free_para_r,ws_hind_opt_r,ws_hind_para_r)
#endif

#ifdef key_contaminant
   ! for the variable 'CONTA' : introduction of contaminant species  
   ! and then possible adition of dissolved or sorbed variables 
   IF (l_conta) THEN
        ALLOCATE(ivdl_conta(nb_conta))
        ivdl_conta(1:nb_conta)=ivdl_conta_tmp(1:nb_conta)
        CALL conta_read_data(obc_cv_name_r,nb_frac_var_r,                         &
                             icon_var_r, init_cv_name_r,ispc_var_r,               &
                             flx_atm_r,cv_rain_r,cini_wat_r,cini_sed_r,           &
                             cini_air_r,l_out_subs_r,tocd_r,ros_r,diam_r,         &
                             ws_free_opt_r,ws_free_min_r,ws_free_max_r,           &
                             ws_free_para_r,ws_hind_opt_r,ws_hind_para_r          &
#ifdef key_species
                            ,frac_cwat_r,frac_csed_r                              &
#endif
                                            )
    END IF   
#endif

   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
#if defined key_biolo || defined key_contaminant || defined key_sedim_MUSTANG || defined key_eco3m
     WRITE(ierrorlog,*)' - Make sure there is no mistake in ',TRIM(name_filesubs),'file !'
#endif
     WRITE(ierrorlog,*)'The simulation is stopped'
     CALL_MPI 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
#if defined key_tssub_macho3d || defined key_parsub_macho3d
   ALLOCATE(c3(nv_tot,kmax,liminm1:limaxp2,ljminm1:ljmaxp2))
   c3(:,:,:,:)=0.0_rsh
#endif
#if defined key_tssub_adv_fifthorder || defined key_tssub_adv_fourthorder
   ALLOCATE(c5(nv_adv,kmax,liminm2:limaxp3,ljminm2:ljmaxp3))
   c5(:,:,:,:)=0.0_rsh
#endif
                                  ! 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
#if defined key_parsub_newdt
      ALLOCATE(ws3max(nvp))     
      ws3max(:)=0.0_rsh
#endif
      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.
#if defined key_sand2D
   ALLOCATE(l_outsandrouse(nvp))
   l_outsandrouse(:)=.false.
#endif

#if defined key_biolo || defined key_contaminant || defined key_eco3m
   ALLOCATE( dcdt(nv_state,kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(attenuation(kmax,limin:limax,ljmin:ljmax))
   ALLOCATE(extinction(kmax,limin:limax,ljmin:ljmax))
   ALLOCATE(alumplafond(kmax,limin:limax,ljmin:ljmax))
   dcdt(:,:,:,:)=0.0_rsh
   attenuation(:,:,:)=0.0_rsh
   extinction(:,:,:)=0.0_rsh
   alumplafond(:,:,:)=0.0_rsh
#endif

#ifdef key_biolo
   ALLOCATE(extinction_ave4d(kmax,limin:limax,ljmin:ljmax))
   ALLOCATE(extinction_aveh(kmax,limin:limax,ljmin:ljmax))
   ALLOCATE(extinction_tab(4,24,kmax,limin:limax,ljmin:ljmax))
   extinction_ave4d(:,:,:)=0.0_rsh
   extinction_aveh(:,:,:)=0.0_rsh
   extinction_tab(:,:,:,:,:)=0.0_rsh
   t_cum_extinctionh=0.0_rsh
   ihour_previous=0
   !ALLOCATE(phytomoy(limin:limax,ljmin:ljmax))
   !ALLOCATE(phytocarre(limin:limax,ljmin:ljmax))
   !phytomoy(:,:)=0.0_rsh
   !phytocarre(:,:)=0.0_rsh
   ALLOCATE(cmes_3dmgl(kmax,limin:limax,ljmin:ljmax))
   cmes_3dmgl(:,:,:)=0.0_rsh
#ifdef key_messat
   ALLOCATE(messat(limin:limax,ljmin:ljmax,51))
   messat(:,:,:)=0.0_rsh
#endif
#ifdef key_oyster_SFG
   ALLOCATE(tpostpontecoq(limin:limax,ljmin:ljmax))
   ALLOCATE(tpostpontegam(limin:limax,ljmin:ljmax))
#endif
#endif

#ifdef key_contaminant
   ALLOCATE(icon_var(nv_adv))
  
#ifdef key_species
   ! allocation des tableaux se referant aux especes de contaminant dans chaque variable
   nfrac=MAXVAL(nb_frac_var_r)
   ALLOCATE(nb_frac_var(nv_adv))
   ALLOCATE(frac_cwat(nv_adv,nfrac,kmax,limin:limax,ljmin:ljmax))
   ALLOCATE(frac_csed(nv_adv,nfrac,kmax,limin:limax,ljmin:ljmax))
   ALLOCATE(dc_spc_wat(nv_adv,nfrac,kmax,limin:limax,ljmin:ljmax))
   ALLOCATE(ispc_var(nv_adv,nfrac))
   ALLOCATE(ispc_var(nv_adv,nfrac))
#else
   ALLOCATE(ispc_var(nv_adv))
#endif
#endif

   ! 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

#ifdef key_biolo
   DO isubs=1,nv_tot
     CALL biolo_init_iv(isubs,standard_name_var(isubs))
   END DO
#endif

   ! 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_MPI (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_MPI
         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_MPI (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_MPI
       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 MPI_FINALIZE(IERR_MPI)
       STOP
     END IF
   END DO

   ! save into simu.log
   !-------------------
   IF_MPI (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_MPI

   ! ---------------------------------------------------------
   ! 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
#ifdef key_sand2D
   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
#endif
   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 defined key_sedim_MUSTANG || defined key_biolo || defined key_conta
   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 MPI_FINALIZE(IERR_MPI)
      STOP
   END IF
#endif

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

#ifdef key_contaminant
     icon_var(iv)=icon_var_r(irk_fil(iv))
#ifdef key_species     
     nb_frac_var(iv)=nb_frac_var_r(irk_fil(iv))
     ispc_var(iv,:)=ispc_var_r(irk_fil(iv),:)
     DO ifracv=1,nb_frac_var(iv)
        frac_cwat(iv,ifracv,:,:,:)=frac_cwat_r(irk_fil(iv),ifracv)
        frac_csed(iv,ifracv,:,:,:)=frac_csed_r(irk_fil(iv),ifracv)
     END DO
#else
     ispc_var(iv)=ispc_var_r(irk_fil(iv))
#endif
#endif
   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 defined key_siggen || defined key_gencoord
     IF (hc > hm) THEN
        DO k=1,kmax
           dhmin=MIN(dhmin,hm*dsigu(k))
        ENDDO
     ELSE
       DO k=1,kmax
         dhmin=MIN(dhmin,(hc+(hm-hc)*dcusds(k))*dsigu(k))
       ENDDO
     END IF
#else
     IF(.NOT.l_modele2d) THEN
       DO k=1,kmax
         dhmin=MIN(dhmin,hm*dsigu(k))
       ENDDO
     ENDIF
#endif

#ifdef key_subs_part_eqsubdt
     DO iv=1,nvp                      ! loop on particulate variables
       IF (ap_ws_gr(iv)==0) THEN

     ! estimate the number of different settling velocities 
         nb_ws_max=nb_ws_max+1
     ! number of variables with the same settling velocity
         nv_ws_max(nb_ws_max)=1
         indgrp=1
     ! index of variable associated with a settling velocity
         iv_ws_max(nb_ws_max,indgrp)=iv
         ap_ws_gr(iv)=nb_ws_max
     ! number of sub time step
         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
             IF_MPI (MASTER) THEN
             WRITE(ierrorlog,*) ''
             WRITE(ierrorlog,*) ' ******* Error detected from SUB_READ_VAR from SUBREADDAT.F90 ******* '
             WRITE(ierrorlog,*) 'the variable :',name_var(irk_fil(isubs))
             WRITE(ierrorlog,*) 'is a particulate one'
             WRITE(ierrorlog,*) 'but its maximum 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 MPI_FINALIZE(IERR_MPI)
             STOP
             ENDIF_MPI
           ELSE
             IF(l_modele2d) THEN
               ndt_part(nb_ws_max)=1
             ELSE
               ndt_part(nb_ws_max)=NINT(dt/MIN((dhmin/ws_free_max(iv)),dt))
             END IF
           END IF
         ENDIF
         sv_part(nb_ws_max)=ws_free_max(iv)
         IF (iv<nvp) THEN
             DO iiv=iv+1,nvp
               IF ((ap_ws_gr(iiv)==0) .AND. (ws_free_max(iiv)==sv_part(nb_ws_max))) then
                 ! number of variables with the same settling velocity
                 nv_ws_max(nb_ws_max)=nv_ws_max(nb_ws_max)+1
                 indgrp=indgrp+1
                 ! index of variable associated with a settling velocity
                 iv_ws_max(nb_ws_max,indgrp)=iiv
                 ap_ws_gr(iiv)=nb_ws_max
               ENDIF
             ENDDO
         ENDIF
       ENDIF
     ENDDO

#else

     nb_ws_max=0
     DO iv=1,nvp
       nb_ws_max=iv     ! number of variables with different settling velocities
                        ! all of them for key_sedim_MUSTANG 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 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

#if defined key_parsub_macho3d
      ALLOCATE(nbouc_machop(nb_ws_max))
      nbouc_machop(:)=0
#endif

   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_MPI (MASTER) THEN
#if defined key_sedim_MUSTANG 
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) 'SUBSTANCES WITH A VARIABLE SETTLING VELOCITY '
#else
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) 'SUBSTANCES WITH A TIME CONSTANT SETTLING VELOCITY '
#endif
     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 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_MPI
   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

#else
   RETURN
#endif

  END SUBROUTINE sub_read_var

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

  SUBROUTINE sub_read_vardiag

#if defined key_substance && defined key_biolo

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE sub_read_vardiag  ***
   !&E
   !&E ** Purpose : lecture du fichier descriptif des variables diagnostiques
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2006-09 (V. Garnier) Original code
   !&E       !  2008-01 (C. Penard)  saving of diagnoses variables
   !&E       !  2008-09 (G. Galibert / ALTRAN) attributes that respect NetCDF OCO v1.1 file format
   !&E                                         reads long_name, standard_name, valid_min and valid_max
   !&E                                         attributes from file_vardiag
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comsubstance, ONLY : filevardiag
   USE comvars2d,    ONLY : l_out_subs_diag
#ifdef key_biolo
   USE bioloinit, ONLY : effetlumiere_day_diat,effetlumiere_day_nano,effetlumiere_day_dino
#ifdef key_karenia
   USE bioloinit, ONLY : effetlumiere_day_karenia
#endif
#ifdef key_psnz
   USE bioloinit, ONLY : effetlumiere_day_psnz
#endif
#ifdef key_phaeocystis
   USE bioloinit, ONLY : effetlumiere_day_phaeocystis
#endif
#ifdef key_N_tracer
   USE parameters, ONLY : nb_var_tracerN 
   USE bioloinit, ONLY : ndiag_tracerN,nb_source_tracerN,isubs_tracer_N,id_tracer_signN,name_source_tracerN
   USE comsubstance, ONLY : name_var,long_name_var,standard_name_var,valid_min_var,valid_max_var
   USE bioloinit, ONLY : biolo_create_vardiagtracer
#ifdef key_age_tracer
   USE bioloinit, ONLY : id_tracer_ageN
#endif
#endif
#ifdef key_P_tracer
   USE parameters, ONLY : nb_var_tracerP 
   USE bioloinit, ONLY : ndiag_tracerP,nb_source_tracerP,isubs_tracer_P,id_tracer_signP,name_source_tracerP
   USE comsubstance, ONLY : name_var,long_name_var,standard_name_var,valid_min_var,valid_max_var
   USE bioloinit, ONLY : biolo_create_vardiagtracer
#ifdef key_age_tracer
   USE bioloinit, ONLY : id_tracer_ageP
#endif
#endif
#endif

   !! * Arguments

   !! * Local declarations
   LOGICAL               :: ex,l_diag_wat,l_diag_sed
   INTEGER               :: eof,isubs,isubs_r,dimvar,it,ind_white,IERR_MPI
   CHARACTER(LEN=lchain) :: namvar_r,long_name_var_r,standard_name_var_r,unitvar_r
   CHARACTER(LEN=5)      :: comment
#if defined key_N_tracer || defined key_P_tracer
   INTEGER               :: is,id,ivtra
#endif

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

   ! save into simu.log
   !-------------------
   IF_MPI (MASTER) THEN
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) '**************************************************'
     WRITE(iscreenlog,*) '**************** VAR_READ_DIAG.F90 ***************'
     WRITE(iscreenlog,*) '**************************************************'
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) 'file defining usefull diagnostic variables : ',TRIM(filevardiag)
   ENDIF_MPI

   ! Initialize number of diagnostic variables according to their dimensions
   ndiag_1d = 0
   ndiag_2d = 0
   ndiag_3d = 0
   ndiag_3d_wat = 0
   ndiag_3d_sed = 0
   ndiag_tot = 0
   isubs = 0
   l_out_subs_diag=.false.  ! no saving of diagnoses variable by defaulf

   ! read total number of diagnostic variables
   eof = 0
   INQUIRE(file=filevardiag,exist=ex)
   IF (ex) THEN
     OPEN(49,file = filevardiag,form='formatted')
     comment='debut'
     DO WHILE (comment /= '*****')
       READ(49,'(a)',iostat=eof) comment
     END DO
     DO WHILE(eof==0)
       READ(49,*,iostat=eof)     !num var
       READ(49,'(a)',iostat=eof) !variable name
       READ(49,'(a)',iostat=eof) !variable long_name
       READ(49,'(a)',iostat=eof) !variable standard_name
       READ(49,'(a)',iostat=eof) !unit
       READ(49,*,iostat=eof)     !variable valid_min value
       READ(49,*,iostat=eof)     !variable valid_max value
       READ(49,*,iostat=eof)     !dimension (1=1D ; 2=2D ; 3=3D)
       READ(49,*,iostat=eof)     !variable in water    (k,i,j)
       READ(49,*,iostat=eof)     !variable in sediment (i,j,k)
       READ(49,*,iostat=eof)     !saving in file
       IF (eof==0) ndiag_tot=ndiag_tot+1
       READ(49,*,iostat=eof)
     END DO
     CLOSE(49)


#if defined key_N_tracer
       ndiag_tracerN=0
       DO is=1,nb_source_tracerN
            ndiag_tot=ndiag_tot+nb_var_tracerN+1
            ndiag_tracerN=ndiag_tracerN+nb_var_tracerN+1
#if defined key_age_tracer
            ndiag_tot=ndiag_tot+nb_var_tracerN+1
            ndiag_tracerN=ndiag_tracerN+nb_var_tracerN+1
#endif
       ENDDO
#endif
#if defined key_P_tracer
       ndiag_tracerP=0
       DO is=1,nb_source_tracerP
            ndiag_tot=ndiag_tot+nb_var_tracerP+1
            ndiag_tracerP=ndiag_tracerP+nb_var_tracerP+1
#if defined key_age_tracer
            ndiag_tot=ndiag_tot+nb_var_tracerP+1
            ndiag_tracerP=ndiag_tracerP+nb_var_tracerP+1
#endif
       ENDDO
#endif

     IF_MPI (MASTER) THEN
        WRITE(iscreenlog,*) 'Number of diagnostic variables = ',ndiag_tot
     ENDIF_MPI
   ! allocate arrays for diagnostic variables
     ALLOCATE( idimv_r(ndiag_tot) )
     ALLOCATE( l_diag_sortie(ndiag_tot) )
     ALLOCATE( name_vardiag(ndiag_tot) )
     ALLOCATE( long_name_vardiag(ndiag_tot) )
     ALLOCATE( standard_name_vardiag(ndiag_tot) )
     ALLOCATE( unit_vardiag(ndiag_tot) )
     ALLOCATE( valid_min_vardiag(ndiag_tot) )
     ALLOCATE( valid_max_vardiag(ndiag_tot) )
     ALLOCATE( irk_diag(ndiag_tot) )

   ! read diagnostic variables and order them according to their matrix dimensions
     OPEN(49,file = filevardiag,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) namvar_r
         ind_white=INDEX(namvar_r,' ')
         name_vardiag(isubs)=TRIM(ADJUSTL(ADJUSTR(namvar_r(1:ind_white))))
         READ(49,'(a)',iostat=eof) long_name_var_r
         ind_white=INDEX(long_name_var_r,' ') 
         long_name_vardiag(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_vardiag(isubs)=TRIM(ADJUSTL(ADJUSTR(standard_name_var_r(1:ind_white))))
         READ(49,'(a)',iostat=eof) unitvar_r
         ind_white=INDEX(unitvar_r,' ')
         unit_vardiag(isubs)=TRIM(ADJUSTL(ADJUSTR(unitvar_r(1:ind_white))))
         READ(49,*,iostat=eof) valid_min_vardiag(isubs)
         READ(49,*,iostat=eof) valid_max_vardiag(isubs)
         READ(49,*,iostat=eof) dimvar
         IF (dimvar==1) THEN
           idimv_r(isubs)=1
           ndiag_1d=ndiag_1d+1
         ELSE IF (dimvar==2) THEN
           idimv_r(isubs)=2
           ndiag_2d=ndiag_2d+1
         ELSE IF (dimvar==3) THEN
           idimv_r(isubs)=3
           ndiag_3d=ndiag_3d+1
         END IF
         READ(49,*,iostat=eof) l_diag_wat
         IF (l_diag_wat .and. dimvar==3) ndiag_3d_wat=ndiag_3d_wat+1
         READ(49,*,iostat=eof) l_diag_sed
         IF (l_diag_sed .and. dimvar==3) THEN
           ndiag_3d_sed=ndiag_3d_sed+1
           idimv_r(isubs)=4
         END IF
         READ(49,*,iostat=eof) l_diag_sortie(isubs)
         READ(49,*,iostat=eof)
         IF (l_diag_sortie(isubs)) l_out_subs_diag=.true.
         IF_MPI (MASTER) THEN
           WRITE(iscreenlog,*)' '
           WRITE(iscreenlog,*)' DIAG. VAR. NAME : ',TRIM(name_vardiag(isubs))
           WRITE(iscreenlog,*)' UNIT            : ',TRIM(unit_vardiag(isubs))
           WRITE(iscreenlog,*)' DIMENSION       : ', dimvar,'D'
           WRITE(iscreenlog,*)' DIAG SAVE ?     : ', l_diag_sortie(isubs)
           IF (l_diag_wat) WRITE(iscreenlog,*)' DIAGNOSTIC INTO THE SEA'
           IF(l_diag_sed) WRITE(iscreenlog,*)' DIAGNOSTIC INTO THE SEDIMENT'
         ENDIF_MPI
       END IF

     END DO
     CLOSE(49)
   END IF     ! fin du test sur l existence du fichier

#if defined key_N_tracer || defined key_P_tracer
       call biolo_create_vardiagtracer
#endif

   ! save into simu.log
   !-------------------
   IF_MPI (MASTER) THEN
     WRITE(iscreenlog,*)' '
     WRITE(iscreenlog,*)' STOCK OF DIAGNOSTIC VARIABLES'
     WRITE(iscreenlog,*)' number of diagnostic variables total :',ndiag_tot
     WRITE(iscreenlog,*)' number of diagnostic variables 1 dim :',ndiag_1d
     WRITE(iscreenlog,*)' number of diagnostic variables 2 dim :',ndiag_2d
     WRITE(iscreenlog,*)' number of diagnostic variables 3 dim :',ndiag_3d
     WRITE(iscreenlog,*)' number of diagnostic variables 3d wat:',ndiag_3d_wat
     WRITE(iscreenlog,*)' number of diagnostic variables 3d sed:',ndiag_3d_sed
   ENDIF_MPI

   IF (ndiag_tot /= ndiag_1d+ndiag_2d+ndiag_3d) THEN
      PRINT*,'WARNING number of diagnostic variables is incoherent.'
      PRINT*,'Check in file :',filevardiag
   END IF

#if ! defined key_sedim_MUSTANG
   IF (ndiag_3d_sed /= 0) THEN
      PRINT*,'WARNING no key_sedim_MUSTANG, you should not have any diagnostic variable for sediment'
      PRINT*,'simulation stopped'
      CALL_MPI MPI_FINALIZE(IERR_MPI)
      STOP
   END IF
#endif

   ! allocate diagnostic variables
   ALLOCATE( diag_1d(1:ndiag_1d) )
   ALLOCATE( diag_2d(ndiag_1d+1:ndiag_1d+ndiag_2d,limin:limax,ljmin:ljmax) )
   ALLOCATE( diag_3d_wat(ndiag_2d+1:ndiag_2d+ndiag_3d_wat,kmax,limin:limax,ljmin:ljmax) )
   diag_1d(:)=0.0_rsh
   diag_2d(:,:,:)=0.0_rsh
   diag_3d_wat(:,:,:,:)=0.0_rsh
#ifdef key_sedim_anc
   ALLOCATE( diag_3d_sed(ndiag_3d_wat+1:ndiag_3d_wat+ndiag_3d_sed,limin:limax,ljmin:ljmax,ksdmin:ksdmax) )
   diag_3d_sed(:,:,:)=0.0_rsh
#endif

#ifdef key_biolo
   !allocate effet lumiere integre sur la journee
   ALLOCATE( effetlumiere_day_diat(kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE( effetlumiere_day_dino(kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE( effetlumiere_day_nano(kmax,limin:limax,ljmin:ljmax) )
   effetlumiere_day_diat(:,:,:)=0.0_rsh
   effetlumiere_day_dino(:,:,:)=0.0_rsh
   effetlumiere_day_nano(:,:,:)=0.0_rsh
#ifdef key_karenia
   ALLOCATE( effetlumiere_day_karenia(kmax,limin:limax,ljmin:ljmax) )
   effetlumiere_day_karenia(:,:,:)=0.0_rsh
#endif
#ifdef key_psnz
   ALLOCATE( effetlumiere_day_psnz(kmax,limin:limax,ljmin:ljmax) )
   effetlumiere_day_psnz(:,:,:)=0.0_rsh
#endif

#ifdef key_phaeocystis
   ALLOCATE( effetlumiere_day_phaeocystis(kmax,limin:limax,ljmin:ljmax) )
   effetlumiere_day_phaeocystis(:,:,:)=0.0_rsh
#endif
#endif

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

   ! Storage of diagnostic variables within reading order
   it = 0
   DO dimvar = 1,4
     DO isubs = 1,ndiag_tot
       IF (idimv_r(isubs) == dimvar) THEN
         it=it+1
!         irk_diag(it)=isubs
         irk_diag(isubs)=it
#if defined key_N_tracer
         IF(isubs <= ndiag_tot-ndiag_tracerN .or. isubs==ndiag_tot) CALL biolo_init_id(isubs,standard_name_vardiag(isubs))
#elif key_P_tracer
         IF(isubs <= ndiag_tot-ndiag_tracerP .or. isubs==ndiag_tot) CALL biolo_init_id(isubs,standard_name_vardiag(isubs))
#else
         CALL biolo_init_id(isubs,standard_name_vardiag(isubs))
#endif
       END IF
     END DO
   END DO

   IF_MPI (MASTER) THEN
     DO isubs = 1,ndiag_tot
        WRITE(iscreenlog,*) isubs,TRIM(name_vardiag(isubs)),idimv_r(isubs),irk_diag(isubs)
     END DO
   ENDIF_MPI

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

  END SUBROUTINE sub_read_vardiag


 END MODULE

