! Read parameters for soil, including layer thicknesses.

  SUBROUTINE init_soil(nia,nja,npatch,nzg,soil_text, &
                       slbs,slcpd,slcons,slden,Wwilt,PHIsat,alfa_vG,slmsts,nstyp,slz)

!DSM {
!! Variaveis provenientes do BRAMS, todas com: intent(in) 
!  USE leaf_coms, ONLY : &
!  !   - b exponent (dimensionless)
!     slbs,   &  
!  !   - dry soil volumetric heat capacity (J/m3/K)   
!     slcpd,  & 
!  !   - saturation soil hydraulic conductivity (m/s)    
!     slcons, & 
!  !   - dry soil density (kg/m3) (also total soil porosity)    
!     slden, &
!  !   - retirado do arquivo util_soil_moisture_init_global_gpcp_trmm.f90
!     Wwilt, &
!  !   - saturation moisture potential (m) - retirado do arquivo util_soil_moisture_init_global_gpcp_trmm.f90
!     PHIsat,alfa_vG &
!  !   - saturation volumetric moisture content (m3/m3)
!     slmsts   
           
!DSM }

  USE ancil_info, ONLY :  &
!  imported scalars with intent(in)
     land_pts,sm_levels  &
!  imported scalars with intent(out)
    ,soil_pts  &
!  imported arrays with intent(out)
    ,soil_index  &
!   imported arrays with intent(inout) !DSM: Para selecionar os pontos de solo que estao no continente
    ,row_length,land_index


  USE file_utils, ONLY :  &
!  imported procedures
     closeFile,fileUnit,findTag,openFile  &
!  imported arrays with intent(inout)
    ,irecPrev

  USE inout, ONLY :  &
!  imported scalar parameters
     formatAsc,formatBin,formatLen,formatNc,formatPP,jinUnit,tagAscBin,tagNc  &
!  imported scalars with intent(in)
    ,echo,nxIn,nyIn  &
!  imported arrays with intent(in)
    ,mapInLand

  USE jules_netcdf, ONLY :  &
!  imported scalars with intent(in)
     ncType

  USE misc_utils, ONLY :  &
!  imported procedures
     checkVarPos,checkVars,read_list,repeatVal,varInfo,varList,varValue

  USE p_s_parms, ONLY :  &
!  imported arrays with intent(out)
    albsoil,b,hcap,hcon,satcon,sathh  &
   ,smvccl,smvcst,smvcwt

  USE readwrite_mod, ONLY :  &
!  imported procedures
     readvar3dcomp

  USE soil_param, ONLY :  &
!  imported arrays with intent(out)
     dzsoil

  USE switches, ONLY :  &
!  imported scalars with intent(in)
     l_triffid,routeOnly  &
!  imported scalars with intent(out)
    ,l_q10,l_soil_sat_down,l_vg_soil,soilhc_method

!-------------------------------------------------------------------------------
  IMPLICIT NONE

  INTEGER, PARAMETER ::  &!  local SCALARS
    nvarMax = 10  !  the number of possible soil characteristic variables
!                         This may be larger than the number needed for any
!                         particular configuration.
!                         Possible variables:
!                         b,sathh,satcon,smvccl,smvcst,smvcwt,hcap,hcon,albsoil,soilType

  INTEGER ::  &!  local SCALARS
    i,ierr       &!  loop counters/work
   ,inUnit       &!  unit used to connect to input file
   ,iorder       &!  work
!  Every possible variable (nvarMax) must have an "iposX" variables declared here
   ,iposAlb      &!  position in list of all possible variables of albsoil
   ,iposB        &!  position in list of all possible variables of b exponent
   ,iposHcap     &!  position in list of all possible variables of hcap
   ,iposHcon     &!  position in list of all possible variables of hcon
   ,iposSatcon   &!  position in list of all possible variables of satcon
   ,iposSathh    &!  position in list of all possible variables of sathh
   ,iposSMCrit   &!  position in list of all possible variables of critical soil moisture
   ,iposSMSat    &!  position in list of all possible variables of saturated soil moisture
   ,iposSMWilt   &!  position in list of all possible variables of wilting point soil moisture
   ,iposSoilType &!  position in list of all possible variables of soil type (number)
   ,iivar,jvar   &!  work !ST replace variable integer ivar to iivar:conflicting with the optional variable integer array "ivar=varFlagTmp"
   ,l            &!  work
   ,nfieldFile   &!  number of fields per time in a file
   ,nheaderField &!  number of header records before each field in file
   ,nheaderFile  &!  number of header records at start of file
   ,nheaderT     &!  number of headers at start of each time
   ,nlineField   &!  work
   ,nsoilVar     &!  number of soil variables required for chosen configuration
!XX Do we need nsoilVar?? Can we just use nvarIn?
   ,nvarIn   &!  number of soil variables that are read in
!                      Often this =nsoilVar, but may be < nsoilVar if we can derive a field
!                      from other fields. or set by some other assumption.
   ,nz           &!  number of levels of data to read
   ,nzz          &!  work
   ,readT        &!  time level to be read from file
   ,useIndex      !  index in irecPrev

  REAL :: albSoilConst   !  a constant (in space and time) value of soil albedo, for use
!                               when useSoilType=TRUE.

  INTEGER :: &!  local arrays
    ival(nvarMax)           &!  work
   ,soilType(land_pts)      &!  soil type (class)
   ,varFlag(nvarMax)        &!  flag indicating source of data.
!                       At present this must either be zero (variable not used - little point!)
!                       or >0 (location of first level of each variable in input file, in terms of
!                       field number, or number of xy planes. Set to 1 for SDF files).
   ,varInSort(nvarMax)    &!  a list of which variables in master list are to be read in,
!                                sorted into ascending order of location in file.
!                                Values 1:varIn give a list of variables that are to be
!                                read from file (in terms of location in masterlist of all possible
!                                variables), sorted into ascending order of location
!                                of data in file.
   ,varFlagTmp(nvarMax)        &!  work: flags as read in
   ,varStashCode(nvarMax)      &!  STASH code for each variable
   ,varUse(nvarMax)             !  flag indicating if a variable
!                        is required for current setup, and how it is used
!                        0: not required
!                        1: needed. Can be set indirectly (from other variables).
!                        2: only needed indirectly to set other variables.

  REAL ::  &!  local ARRAYS
    tmpval(land_pts,sm_levels)   !  workspace, large enough for all levels of a variable

  LOGICAL ::  &!  local SCALARS
    checkNames,checkPos  &!  work
   ,constZ       &!  flag indicating that soil characteristics are constant with depth
   ,errFound     &!  flag indicating an error
   ,levels       &!  work
   ,readFile     &!  flag indicating if another file is to be read
   ,sdfFile      &!  TRUE if dealing with an SDF (self-describing file)
   ,summary      &!  work
   ,useSoilType  &!  flag indicating if a map of soil types is to be read
!                      T read a map of soil types, and characteristics for each type
!                      F read maps of soil characteristics
   ,zrev          !  T means the order of the soil layers in the input data will be reversed.
                  !    Soil levels start from that at the ground surface and increase with depth (top to bottom).
                  !    F means the input soil data are arranged top to bottom
                  !    T means the input soil data are arranged bottom to top.

  LOGICAL ::  &!  arrays
    done(nvarMax)      &!  work
   ,foundVar(nvarMax)   !  flag indicating if a variable is listed in run control file

  CHARACTER(len=formatLen) ::  &! local SCALARS
    fileFormat   !  format of file

  CHARACTER(len=150) ::  &! local SCALARS
    fileName     &!  the name of a file
   ,soilLUTfile   !  file containing look-uo table of soil characteristics

  CHARACTER(len=100) ::  &! local ARRAYS
    varDesc(nvarMax)      &!  description of each variable
   ,varExtra(nvarMax)     &!  "extra" information for of each variable
   ,varName(nvarMax)      &!  names of all possible variables - used to identify
!                                    variables in this routine
   ,varNameTmp(nvarMax)   &!  work: variable names as read in
   ,varNameSDF(nvarMax)   &!  names of variables in input file (SDF=self-describing file)
   ,varNameSDFTmp(nvarMax) !  work: names as read in

!DSM {
   INTEGER, INTENT(IN)     :: nia,nja,npatch,nzg
   REAL,    INTENT(IN)     :: soil_text(nzg,nia,nja),slz(nzg)
   INTEGER                 :: j,k,nstyp,isoil
   REAL, dimension(nstyp)  :: slbs,slcpd,slcons,slden,Wwilt,PHIsat,alfa_vG,slmsts

   REAL, dimension(3) ::  beta,p ! Ahuja et al., 1984 Ksat calculation  !DSM

   REAL, dimension(12) :: sand,clay,silt,alfa2vG,n2vG,ssmsts,soilcp,soilef,sscons
      
   REAL                   :: Bp,SATHHp,V_SAT

   REAL,parameter::CAP_CLAY=2.373E6, CAP_SAND=2.133E6, &
                   CAP_SILT=2.133E6, CAP_SNOW=0.63E6, &
                   CON_CLAY=1.16, CON_SAND=1.57,      & 
                   CON_SILT=1.57, CON_SNOW=0.265, CON_AIR=0.025

   DATA  sand  /.97,.92,.80,.57,.60,.65,.35,.48,.50,.30,.25,.20/
   DATA  clay  /.03,.07,.18,.40,.35,.31,.59,.45,.42,.61,.65,.20/
   DATA  silt  /.00,.01,.02,.03,.05,.04,.06,.07,.08,.09,.10,.60/

   DATA  alfa2vG  /0.380,0.837,0.396,0.191,0.246,0.644,0.298,0.392,0.509,0.258,0.463,0.258/
   DATA  n2vG     /2.474,1.672,1.553,1.644,1.461,1.535,1.513,1.437,1.396,1.466,1.514,1.466/
   DATA  ssmsts   /0.410,0.438,0.461,0.601,0.521,0.413,0.586,0.519,0.460,0.570,0.546,0.570/
   DATA  soilcp   /0.037,0.062,0.111,0.223,0.155,0.149,0.267,0.226,0.199,0.278,0.267,0.278/

   DATA beta/29518.,10584.,56540./
   DATA   p/2.4734,3.3545,4.5399/

    soilef(:) = ssmsts(:) - (soilcp(:) + (ssmsts(:)-soilcp(:) )  /&
                            ( 1. + (alfa2vG(:)*33.)**n2vG(:) )**(1.-1./n2vG(:)))
    sscons(:) = 0.05*beta(3)*(soilef(:)**(p(3))) ! 3= escolha da coluna 3
    sscons(:) = sscons(:)*1.e-3*(1./3600.) * 1000.0 ! converte de mm/h para m/s
    sscons(:) = sscons(:)*1000.0  ! DSM  - multiplicando pela densidade da agua para ter kg m-2 s-1
     
!DSM }

!------------------------------------------------------------------------------------------

! If data are not needed, nothing to do.
  IF ( routeOnly ) RETURN

  if (echo) WRITE(*,"(50('-'),/,a)") 'init_soil'

! Locate the start of this section in input file.
  CALL findTag( jinUnit,'init_soil','>INIT_SOIL' )

! Start to read from run control file.
  READ(jinUnit,*) l_vg_soil
  READ(jinUnit,*) l_soil_sat_down
  READ(jinUnit,*) l_q10
  READ(jinUnit,*) soilhc_method
  READ(jinUnit,*) useSoilType
  READ(jinUnit,*) constZ,zrev

! Establish where data are to be read from.
  READ(jinUnit,*) readFile
  READ(jinUnit,*) fileFormat
  READ(jinUnit,*) fileName
  IF ( useSoilType ) READ(jinUnit,*) soilLUTfile

! Set file format if reading from run control file.
  IF ( .NOT. readFile ) fileFormat = formatAsc

! Insist on zrev=FALSE if reading from run control file.
  IF ( .NOT.readFile .AND. zrev ) THEN
    WRITE(*,*)'ERROR: init_soil: NOT readFile requires zrev=FALSE.'
    STOP
  ENDIF

! If reading a map of soil types, insist on zrev=FALSE (for clarity, hopefully).
  IF ( useSoilType .AND. zrev ) THEN
    WRITE(*,*)'ERROR: init_soil: if using a map of soil types,'
    WRITE(*,*)'zrev must be FALSE.'
    STOP
  ENDIF

! Check soilHc_method recognised.
  IF ( soilHc_method<1 .OR. soilHc_method>2 ) THEN
    WRITE(*,*)'ERROR: init_soil: soilHc_method must be 1 or 2.'
    STOP
  ENDIF

!-------------------------------------------------------------------------------
! Set up list of all possible variables.
! Delayed until now so we know value of l_vg_soil.
! Note that some alternative STASH codes could be used - e.g. STASH 43 is
! saturation moisture after timestep.
!------------------------------------------------------------------------------
  i = 0

  i=i+1; iposB=i; varName(i)='b'
  varStashcode(i) = 8288      !  PP code 1380
  varDesc(i)='Clapp-Hornberger b exponent'
  IF ( l_vg_soil ) THEN
    varStashcode(i) = 8288   !  288 is really Clapp-Hornberger b. PP code 1381
    varDesc(i)='exponent, b=1/(n-1) (m)'
  ENDIF

  i=i+1; iposSathh=i; varName(i)='sathh'
  varStashcode(i) = 8224      !  PP code 342
  varDesc(i)='Saturated soil water pressure (m)'
  IF ( l_vg_soil ) varDesc(i)='1/alpha'

  i=i+1; iposSatcon=i; varName(i)='satcon'
  varStashcode(i) = 8214    !  PP code 333
  varDesc(i)='Saturated hydraulic conductivity (kg m-2 s-1)'

  i=i+1; iposSMSat=i; varName(i)='sm_sat'
  varStashcode(i) = 8213  !  PP code 332
  varDesc(i)='Volumetric soil moisture content at saturation'

  i=i+1; iposSMCrit=i; varName(i)='sm_crit'
  varStashcode(i) = 8211  !  PP code 330
  varDesc(i)='Volumetric soil moisture content at critical point'

  i=i+1; iposSMWilt=i; varName(i)='sm_wilt'
  varStashcode(i) = 8210  !  PP code 329
  varDesc(i)='Volumetric soil moisture content at wilting point'

  i=i+1; iposHcap=i; varName(i)='hcap'
  varStashcode(i) = 8216   !  PP code 335
  varDesc(i)='Soil heat capacity (J K-1 m-3)'

  i=i+1; iposHcon=i; varName(i)='hcon'
  varStashcode(i) = 8217   !  PP code 336
  varDesc(i)='Soil thermal conductivity (W m-1 K-1)'

  i=i+1; iposAlb=i; varName(i)='albsoil'
  varStashcode(i) = 220  !  PP code 1395
  varDesc(i)='Soil albedo'

  i=i+1; iposSoilType=i; varName(i)='soilType'
  varStashcode(i) = -1   !  not known
  varDesc(i)='Soil type number'

!-------------------------------------------------------------------------------
! Indicate the variables that are required for the selected configuration.
!------------------------------------------------------------------------------

! Initialise to show that no variables are required.
  varUse(:) = 0

! Variables that are always required.
  varUse(iposAlb) = 1
  varUse(iposB) = 1
  varUse(iposHcap) = 1
  varUse(iposHcon) = 1
  varUse(iposSatcon) = 1
  varUse(iposSathh) = 1
  varUse(iposSMcrit) = 1
  varUse(iposSMsat) = 1
  varUse(iposSMwilt) = 1

! Variables that are only required under some configurations.
  IF ( useSoilType ) varUse(iposSoilType) = 1

! Count number of variables required for chosen configuration.
! Note that this curently includes "temporary" variables such as soil type.
  nsoilVar = COUNT( varUse(:) > 0 )

!-------------------------------------------------------------------------------
! Only read parameters for the file format indicated.

  SELECT CASE ( fileFormat )

    CASE ( formatAsc,formatBin,formatPP )
!     Locate the information in run control file.
      CALL findTag( jinUnit,'init_soil',tagAscBin,preInit=.TRUE. )
      READ(jinUnit,*) nheaderFile,nheaderField

!     Read details of variables. Read 2 values from a blank-delimited list.
      CALL read_list( jinUnit,2,nvarMax,'>VARS','>ENDVARS',' ','init_soil'  &
                     ,nvarIn,cvar1=varNameTmp,cvar1Pos=1,ivar=varFlagTmp,ivarPos=2 )
      sdfFile = .FALSE.

    CASE ( formatNc )
!     Locate the information in run control file.
      CALL findTag( jinUnit,'init_soil',tagNc,preInit=.TRUE. )

!     Read details of variables. Read 2 values from a blank-delimited list.
      CALL read_list( jinUnit,2,nvarMax,'>VARS','>ENDVARS',' ','init_soil'  &
                     ,nvarIn,cvar1=varNameTmp,cvar1Pos=1,cvar2=varNameSDFTmp,cvar2Pos=2 )

      sdfFile = .TRUE.

    CASE default
      WRITE(*,*)'ERROR: init_soil: no code for fileFormat=',TRIM(fileFormat)
      STOP

  END SELECT

!-------------------------------------------------------------------------------

! If using run control file, set number of headers to zero.
  IF ( .NOT. readFile ) THEN
    nheaderFile = 0
    nheaderField = 0
  ENDIF

! Check that there are no repeated names in input list.
  IF ( repeatVal( varNameTmp(1:nvarIn) ) ) THEN
    WRITE(*,*)'ERROR: init_soil: repeated variable names.'
    STOP
  ENDIF

! Work out what variables are indicated.
  CALL varList( nvarIn,sdfFile,varNameTmp,varName  &
               ,errFound,foundVar  &
               ,varFlagIn=varFlagTmp,varNameSDFin=varNameSDFtmp  &
               ,varFlag=varFlag,varNameSDF=varNameSDF )
  IF ( errFound ) THEN
    WRITE(*,*)'ERROR: init_soil: error returned by varList'
    STOP
  ENDIF

! If using a map of soil type, indicate that all other variables will be "derived",
! meaning set from soil characteristic file.
  IF ( useSoilType ) THEN
    WRITE(*,*)'WARNING: init_soil: as a soil type map is being used,'
    WRITE(*,*)'all other variables will be set based on soil type.'
    DO iivar=1,nvarMax
      IF ( varName(iivar) /= varName(iposSoilType) ) THEN
        varFlag(iivar) = -2
        varExtra(iivar) = 'using characteristics for each soil type'
       ENDIF
    ENDDO
  ENDIF

! Check that all necessary variables are provided.
  CALL checkVars( nvarMax,.FALSE.,.FALSE.,varUse,foundVar,varDesc,varName  &
                 ,nvarIn,varFlag,errFound )
! Stop if an error was detected.
  IF ( errFound ) THEN
    WRITE(*,*)'ERROR: init_soil: an error was found by checkVars.'
    STOP
  ENDIF

!-------------------------------------------------------------------------------
! Get a sorted list, giving input variables according to position in input file(s).
! This is needed for reading from run control file (and makes ASCII reading easier)
! and is always done.
! Note that we now insist that variables from run control file are listed in the
! correct order, but this sorting code has been left - easier for now!
  done(:) = .TRUE.
  varInSort(:) = 0
! Only consider variables that are to be read in (flag>0). SDF variables have flag set >0.
  WHERE ( varFlag(:) > 0 ) done(:)=.FALSE.
  DO iivar=1,nvarIn
!    IF ( readFile ) THEN
      ival(1:1) = MINLOC( varFlag(:), .NOT.done(:) )
!    ELSE
!     Hack so we can use existing code. If reading from run control file, don't
!     sort the list.
!      ival(1) = iivar
!    ENDIF
    i = ival(1)
    done(i) = .TRUE.
    varInSort(iivar) = i
  ENDDO

!-------------------------------------------------------------------------------
! Some output to screen.
  IF ( echo ) THEN
    DO iivar=1,nvarMax
      IF ( varUse(iivar) > 0 ) CALL varInfo( varName(iivar)  &
                      ,varFlag(iivar),varStashcode(iivar),varNameSDF(iivar)  &
                      ,fileFormat,varDesc=varDesc(iivar) )
    ENDDO
  ENDIF
  IF ( echo ) THEN
    DO iivar=1,nvarMax
      IF ( varUse(iivar) /= 0 ) CALL varInfo( varName(iivar)  &
                      ,varFlag(iivar),varStashcode(iivar),varNameSDF(iivar)  &
                      ,fileFormat,varDesc=varDesc(iivar),extraDesc=varExtra(iivar) )
    ENDDO
  ENDIF
!-------------------------------------------------------------------------------

! Check we have an acceptable file format, and set flags indicating what other
! variables to check.
  checkNames = .FALSE.
  checkPos = .FALSE.
  SELECT CASE ( fileFormat )
    CASE ( formatAsc )
!     ASCII. Need positions of variables.
      checkPos = .TRUE.
    CASE ( formatBin )
!     GrADS. For now, need positions of variables.
      checkPos = .TRUE.
    CASE ( formatNc )
!     netCDF. Need names of variables.
      checkNames = .TRUE.
    CASE ( formatPP )
!     PP format. For now, need positions of variables.
      checkPos = .TRUE.
    CASE default
      WRITE(*,*)'ERROR: init_soil: no code for fileFormat=',TRIM(fileFormat)
      STOP
  END SELECT

  IF ( checkPos ) THEN
!   Check positions of variables within file.
!   Check for repeated locations.
!XX   Note that we should really be checking that there is no overlap if nz>1 and varPos is FIELD num (not var num).

!   First argument to checkVarPos indicates how important the order of the variables is.
!   0 means order is not important. If reading from run control file, insist that variables
!   are 1,2,3,.... - indicated by iorder=2. Note that for run control file, we interpret
!   varFlag as variable number, not field number, since this is required for current code
!   to pass checkVarPos with iorder=2 (want vars to appear adjacent, which wouldn't be if
!   had multi-level data and gave non-consecutive field numbers. Admittedly, this is
!   becoming a bit of a faff, and needs to be rewritten, one day.).
!   In all cases, values <1 are ignored.
    iorder = 0
    IF ( .NOT. readFile ) iorder = 2

!   We need to load file locations in sorted order.
    ival(:) = 0
    DO iivar=1,nvarIn
      ival(iivar) = varFlag( varInSort(iivar) )
    ENDDO

    ierr = checkVarPos( iorder,ival(:),' init_soil: ival' )
    IF ( ierr < 0 ) THEN
      WRITE(*,*)'ERROR: init_soil: error from checkvarPos.'
      WRITE(*,*)'If error was repeated use of same varPos, but you do in fact want to reuse the'
      WRITE(*,*)'same data for more than one variable, comment out this stop!'
      STOP
    ENDIF
  ENDIF

  IF ( checknames ) THEN
!   Check names only for variables that are to be read in (varFlag>0).
    DO iivar=1,nvarMax
      DO jvar=iivar+1,nvarMax
        IF ( ( varFlag(iivar)>0 .AND. varFlag(jvar)>0 ) .AND.  &
             ( varNameSDF(iivar)==varNameSDF(jvar) ) ) THEN
          WRITE(*,*)'ERROR: init_soil: repeated varNameSDF: ',TRIM(varNameSDF(iivar))
          WRITE(*,*)'If you really do want to use the same variable from file to set values'
          WRITE(*,*)'of more than one FORTRAN variable, comment out this check!'
          WRITE(*,*)'Stopping in init_soil'
          STOP
        ENDIF
      ENDDO
    ENDDO
  ENDIF

!------------------------------------------------------------------------------
! Read layer thicknesses.
!------------------------------------------------------------------------------

! Find section in the run control file.
  inUnit = jinUnit
  CALL findTag( inUnit,'init_soil dzsoil','>DATA_DZSOIL',preInit=.TRUE. )

! Read data.
  READ(jinUnit,*) dzsoil(:)
  
!--- Utilizando o definido no RAMSIN {
  dzsoil(1)=-1.0*slz(nzg)
  DO k=2,nzg
      dzsoil(k)=-1.0*slz(nzg-k+1) + slz(nzg-k+2)
  ENDDO
!---}

! Read constant soil albedo, if required.
  IF ( useSoilType ) READ(jinUnit,*) albSoilConst

!------------------------------------------------------------------------------
! Open file.
  IF ( readFile ) THEN
    inUnit = fileUnit( fileFormat )  !  get unit
    !DSM CALL openFile( 1,.FALSE.,inUnit,'read',fileFormat,fileName,'old','init_soil',ncType )
  ELSE
    WRITE(*,*)'Reading soil characteristics from the run control file.'
    WRITE(*,*)'Data must be in fields 1,2,3,.... (i.e. sequential from 1).'
    inUnit = jinUnit
    CALL findTag( inUnit,'init_soil char','>DATA',preInit=.TRUE. )
  ENDIF

! Set the number of levels of data to read.
  nz = sm_levels
  IF ( constZ ) nz = 1

!-------------------------------------------------------------------------------

! Read data.
! All levels of a variable are read into temporary space, then loaded into final variables.
! Only data for the required points are returned.
! Only need to do anything for variables that are to be read in, since currently all (required)
! variables must be read in (no option to set as constant).
! For non-SDF files, variables are read in the order in which they appear in file - this is
! done so that values can be read from the run control file (which is possibly open on stdIn and
! therefore can't be backspaced).

  DO iivar=1,nvarIn

!   Get location in master list of next variable to be read.
    jvar = varInSort(iivar)

    IF ( varUse(jvar)>0 .AND. varFlag(jvar)>0 ) THEN

!     The netCDF reading code can't cope with a call from here to a soil with more
!     than one dimension (Land). So, for now, stop if hoping to read more than one level!
      IF ( fileFormat==formatNC .AND. .NOT.constZ ) THEN
        WRITE(*,*)'ERROR: init_soil: can''t read a netCDF file with more than one'
        WRITE(*,*)'dimension from init_soil. Needs more code! Sorry.'
        STOP
      ENDIF

      nzz = nz
      IF ( jvar == iposAlb ) nzz=1

      IF ( inUnit==jinUnit .AND. nxIn*nyIn==1 .AND. fileFormat==formatAsc) THEN
!       If the input grid is a single point and reading from run control file, expect no
!       new line between fields (eg all on one line).
!       Calling readVar means we could cope with headers in the run control file.
!       But there's no need since annotation is already simple in this case.
        READ(jinUnit,*) tmpval(:,1:nzz)

      ELSE

!       Simplifying assumptions regarding input file. Otherwise have to read these in.
        readT = 1             !  time level to read from file
        nfieldFile = varFlag(jvar)+nzz-1 !  # of fields in file.
!                                           Set to last field number - OK while readT=1
        nheaderT = 0          !  no headers at top of each time
        nlineField = 0        !  0 means will not attempt to read ASCII line-by-line

!       Set index to use for irecPrev with netCDF files - this irecPrev isn't changed,
!       but need to keep index within bounds.
        useIndex = inUnit
        IF ( fileFormat == formatNC ) useIndex = 1

        !DSM CALL readVar3dComp( readT,varFlag(jvar),varStashcode(jvar),.FALSE.,irecPrev(useIndex)  &
        !DSM                  ,nfieldFile,nheaderFile,nheaderT  &
        !DSM                  ,nheaderField,nlineField,nxIn,nyIn,nzz,zrev  &
        !DSM                  ,inUnit,varNameSDF(jvar)  &
        !DSM                  ,mapInLand(:),(/(i,i=1,land_pts)/)  &
        !DSM                  ,fileFormat,tmpval(:,1:nzz),'init_soil','init_soil',ncType )

        !DSM {
         
         DO k=1,nzz
            DO l=1,land_pts
               j = ( land_index(l)-1 ) / row_length + 1
               i = land_index(l) - ( j-1 ) * row_length
               
               isoil=int(soil_text(k,i,j))
               IF ( jvar == iposAlb ) THEN
                  tmpval(l,k)=0.110000

                 !--Imprimindo para utilizar no off-line---
                 !write(42,*) 'soil_Alb=>', k,l,tmpval(l,k)
                 
               ELSEIF ( jvar == iposB ) THEN
                  IF (l_vg_soil) THEN
                     tmpval(l,k)=1/(n2vG(isoil)-1)
                  ELSE
                     tmpval(l,k)=3.10+15.70*CLAY(isoil)-0.3*SAND(isoil)
                  !tmpval(l,k)=slbs(int(soil_text(k,i,j))) ! utilizando sempre patch=2, pois soil_text nao possui variacao em relacao ao patch 
                  !tmpval=6.63  !=slbs
                  ENDIF

                 !--Imprimindo para utilizar no off-line---
                 !write(42,*) 'soil_B=>', k,l,tmpval(l,k)
                 
               ELSEIF ( jvar == iposHcap ) THEN
                  V_SAT=0.505-0.037*CLAY(isoil)-0.142*SAND(isoil)
                  tmpval(l,k)=(1-V_SAT)*(CLAY(isoil)*CAP_CLAY+SAND(isoil)*CAP_SAND+SILT(isoil)*CAP_SILT)
                  
                  !print*,'AAAA-hcap=',isoil,tmpval(l,k)
                  if (isoil==6) tmpval(l,k)=4200000                  !from eq.=1321725.
                  !tmpval(l,k)=slcpd(int(soil_text(k,i,j))) 
                  !tmpval=1185676.  !=slcpd

                 !--Imprimindo para utilizar no off-line---
                 !write(42,*) 'soil_Hcap=>', k,l,tmpval(l,k)
                 
               ELSEIF ( jvar == iposHcon ) THEN
                  V_SAT=0.505-0.037*CLAY(isoil)-0.142*SAND(isoil)
                  tmpval(l,k)=(CON_AIR**V_SAT) *(CON_CLAY**(CLAY(isoil)*(1-V_SAT))) &
                                               *(CON_SAND**(SAND(isoil)*(1-V_SAT))) &
                                               *(CON_SILT**(SILT(isoil)*(1-V_SAT)))

                  !tmpval(l,k)=0.226873  !=???

                 !--Imprimindo para utilizar no off-line---
                 !write(42,*) 'soil_Hcon=>', k,l,tmpval(l,k)
                 
               ELSEIF ( jvar == iposSatcon ) THEN
                  IF (l_vg_soil) THEN
                     tmpval(l,k)=sscons(isoil)  ! utilizando l_vg_soil=TRUE - alfa_vG em (1/m)
                  ELSE
                     tmpval(l,k)=10.0**(-0.60-0.64*CLAY(isoil)+1.26*SAND(isoil))*25.4/3600.0
                     
                     !tmpval(l,k)=slcons(int(soil_text(k,i,j)))*slden(int(soil_text(k,i,j))) 
                     !tmpval=0.004715 !=slcons*slden
                  ENDIF

                 !--Imprimindo para utilizar no off-line---
                 !write(42,*) 'soil_Satcon=>', k,l,tmpval(l,k)
                 
               ELSEIF ( jvar == iposSathh ) THEN
                  IF (l_vg_soil) THEN
                     ! convert alfa to 1/m  -  vG parameter (1/kPa) to 1/m (1m=9.8kPa)
                     tmpval(l,k)=9.8*alfa2vG(isoil)  ! utilizando l_vg_soil=TRUE - alfa_vG em (1/m)
                  ELSE
                     tmpval(l,k)=0.01*10.0**(2.17-0.63*CLAY(isoil)-1.58*SAND(isoil))
                     !tmpval(l,k)=PHIsat(int(soil_text(k,i,j)))  ! utilizando l_vg_soil=FALSE - PHIsat em (m)
                     !tmpval=0.049460 
                  ENDIF

                 !--Imprimindo para utilizar no off-line---
                 !write(42,*) 'soil_Sathh=>', k,l,tmpval(l,k)
                 
               ELSEIF ( jvar == iposSMCrit ) THEN
                  SATHHp=0.01*10.0**(2.17-0.63*CLAY(isoil)-1.58*SAND(isoil))
                  Bp=3.10+15.70*CLAY(isoil)-0.3*SAND(isoil)
                  V_SAT=0.505-0.037*CLAY(isoil)-0.142*SAND(isoil)
                  tmpval(l,k)=V_SAT*(SATHHp/3.364)**(1.0/Bp)
                  !tmpval(l,k)=0.242433 !=???

                 !--Imprimindo para utilizar no off-line---
                 !write(42,*) 'soil_SMCrit=>', k,l,tmpval(l,k)
                 
               ELSEIF ( jvar == iposSMSat ) THEN
                  tmpval(l,k)=0.505-0.037*CLAY(isoil)-0.142*SAND(isoil)
                  !tmpval(l,k)=slmsts(int(soil_text(k,i,j)))
                  !tmpval=0.458150

                 !--Imprimindo para utilizar no off-line---
                 !write(42,*) 'soil_SMSat=>', k,l,tmpval(l,k)
                 
               ELSEIF ( jvar == iposSMWilt ) THEN
                  SATHHp=0.01*10.0**(2.17-0.63*CLAY(isoil)-1.58*SAND(isoil))
                  Bp=3.10+15.70*CLAY(isoil)-0.3*SAND(isoil)
                  V_SAT=0.505-0.037*CLAY(isoil)-0.142*SAND(isoil)
                  tmpval(l,k)=V_SAT*(SATHHp/152.9)**(1.0/Bp)
                  !tmpval(l,k)=Wwilt(int(soil_text(k,i,j)))  !=Wwilt em util_soil_moisture_init_global_gpcp_trmm.f90
                  !tmpval=0.136328

                 !--Imprimindo para utilizar no off-line---
                 !write(42,*) 'soil_SMWilt=>', k,l,tmpval(l,k)
                 
               ELSE
                  WRITE(*,*)'ERROR: init_soil: no code to load jvar=',jvar
                  STOP
               ENDIF
               

               
             ENDDO  !l=1,land_pts
          ENDDO  !k=1,nzz
        !DSM }
      ENDIF

!     Copy data into final variable, if necessary also filling all levels with
!     values from a single level.
      IF ( jvar == iposAlb ) THEN
        albSoil(:) = tmpval(:,1)
      ELSEIF ( jvar == iposB ) THEN
        b(:,1:nzz) = tmpval(:,1:nzz)
        b(:,nzz+1:sm_levels) = SPREAD( b(:,1), 2, sm_levels-nzz )
      ELSEIF ( jvar == iposHcap ) THEN
        hcap(:,1:nzz) = tmpval(:,1:nzz)
        hcap(:,nzz+1:sm_levels) = SPREAD( hcap(:,1), 2, sm_levels-nzz )
      ELSEIF ( jvar == iposHcon ) THEN
        hcon(:,1:nzz) = tmpval(:,1:nzz)
        hcon(:,nzz+1:sm_levels) = SPREAD( hcon(:,1), 2, sm_levels-nzz )
      ELSEIF ( jvar == iposSatcon ) THEN
        satcon(:,1:nzz) = tmpval(:,1:nzz)
        satcon(:,nzz+1:sm_levels) = SPREAD( satcon(:,1), 2, sm_levels-nzz )
      ELSEIF ( jvar == iposSathh ) THEN
        sathh(:,1:nzz) = tmpval(:,1:nzz)
        sathh(:,nzz+1:sm_levels) = SPREAD( sathh(:,1), 2, sm_levels-nzz )
      ELSEIF ( jvar == iposSMCrit ) THEN
        smvccl(:,1:nzz) = tmpval(:,1:nzz)
        smvccl(:,nzz+1:sm_levels) = SPREAD( smvccl(:,1), 2, sm_levels-nzz )
      ELSEIF ( jvar == iposSMSat ) THEN
        smvcst(:,1:nzz) = tmpval(:,1:nzz)
        smvcst(:,nzz+1:sm_levels) = SPREAD( smvcst(:,1), 2, sm_levels-nzz )
      ELSEIF ( jvar == iposSMWilt ) THEN
        smvcwt(:,1:nzz) = tmpval(:,1:nzz)
        smvcwt(:,nzz+1:sm_levels) = SPREAD( smvcwt(:,1), 2, sm_levels-nzz )
      ELSEIF ( jvar == iposSoilType ) THEN
        soilType(:) = NINT( tmpval(:,1) )
      ELSE
        WRITE(*,*)'ERROR: init_soil: no code to load jvar=',jvar
        STOP
      ENDIF

!   Here we would have endif for varpos<0 (meaning derive or otherwise set field)
    ENDIF   !  read variable

  ENDDO  !  variables

! Close file if it is not the JULES control file
  IF ( inUnit /= jinUnit ) CALL closeFile( inUnit,fileFormat )

!-------------------------------------------------------------------------------
! Read look-up table of characteristics for each soil type and set fields.
!------------------------------------------------------------------------------
  IF ( useSoilType ) CALL init_soil_lut( land_pts,nsoilVar,sm_levels  &
                          ,albSoilConst,soilLUTfile,soilType )

!-------------------------------------------------------------------------------
! Set surface values.
  hcon(:,0) = hcon(:,1)
  satcon(:,0) = satcon(:,1)

! Check that sathh>=0 - a common error!
  IF ( ANY( sathh(:,:)<0.0 ) ) THEN
    WRITE(*,*)'sathh < 0.0 detected. Should be > 0.'
    WRITE(*,*)'For JULES, sathh is abs(saturated head).'
    WRITE(*,*)'Stopping in init_soil'
    STOP
  ENDIF

!-------------------------------------------------------------------------------
! Detect soil points.
! If top layer saturation moisture content>0, this is a soil point.
! Note that land ice points are no longer assigned here.
!-------------------------------------------------------------------------------
  SOIL_PTS=0
  soil_index(:) = 0
  DO L=1,LAND_PTS
    IF( SMVCST(L,1) > EPSILON(SMVCST(L,1)) ) THEN
      SOIL_PTS=SOIL_PTS+1
      SOIL_INDEX(SOIL_PTS)=L
    ENDIF
  ENDDO

  if (echo) WRITE(*,*)'Number of soil points=',soil_pts

  IF ( soil_pts == 0 ) THEN
    WRITE(*,"(10('*'),a,10('*'))")' WARNING: init_soil: there are no soil points. '
    WRITE(*,"(10('*'),a,10('*'))")' Any land points are land ice points. '
  ENDIF

!------------------------------------------------------------------------------
! Set l_q10 to TRUE if not using TRIFFID.
! This is anyway done later in MICROBE, so do here too for ease.
!------------------------------------------------------------------------------
  IF ( .NOT. l_triffid ) l_q10 = .TRUE.

!------------------------------------------------------------------------------
! Optional write to screen.
!------------------------------------------------------------------------------
  IF ( echo ) THEN

    IF ( l_vg_soil ) WRITE(*,*) 'van Genuchten model will be used'
    IF ( l_soil_sat_down ) THEN
      WRITE(*,*) 'l_soil_sat_down=TRUE, meaning excess water is pushed down'
    ELSE
      WRITE(*,*) 'l_soil_sat_down=FALSE, meaning excess water is pushed up'
    ENDIF
    IF ( soilHc_method == 1 ) THEN
      WRITE(*,*)'soilHc_method=1:  following Cox et al (1999)'
    ELSE
      WRITE(*,*)'soilHc_method=2:  following simplified Johansen (1975)'
    ENDIF
    IF ( l_q10 ) THEN
      WRITE(*,*) 'Q10 equation will be used for soil respiration.'
    ELSE
      WRITE(*,*) 'RothC equations will be used for soil respiration.'
    ENDIF

    WRITE(*,*)'Soil layer thicknesses (m):'
    DO i=1,sm_levels
      WRITE(*,"(a,i2,a,(10f5.2))")'iz=',i,' dzsoil=',dzsoil(i)
    ENDDO

!-------------------------------------------------------------------------------
!   Write fields to screen.
!-------------------------------------------------------------------------------
    DO i=1,2
!     i = 1 writes full field, level by level
!     i = 2 writes summary statistics for whole field
      IF ( nxIn*nyIn==1 .AND. i==2 ) EXIT
      IF ( i == 1 ) THEN
        summary = .FALSE.
        levels = .TRUE.
        if (echo) WRITE(*,*)'###############################################'
      ELSE
        summary = .TRUE.
        levels = .FALSE.
        WRITE(*,*)'###############################################'
        WRITE(*,*)'Summary of input.'
        WRITE(*,*)'NB The ranges below include any ice points.'
      ENDIF

      DO iivar=1,nvarMax

        IF ( varUse(iivar) == 1 ) THEN

          IF ( iivar == iposAlb ) THEN
            CALL varValue( summary,albsoil,varFormat='f6.3',varName=varName(iivar) )

          ELSEIF ( iivar == iposB ) THEN
            CALL varValue( levels,summary,b,varFormat='f7.3',varName=varName(iivar) )

          ELSEIF ( iivar == iposHcap ) THEN
            CALL varValue( levels,summary,hcap,varFormat='es10.2',varName=varName(iivar) )

          ELSEIF ( iivar == iposHcon ) THEN
            CALL varValue( levels,summary,hcon,varFormat='f6.3',varName=varName(iivar) )

          ELSEIF ( iivar == iposSatcon ) THEN
            CALL varValue( levels,summary,satcon,varFormat='es10.2',varName=varName(iivar) )
            IF ( soil_pts /= land_pts ) WRITE(*,"(/,a,/)")  &
                   'NB satcon will be reset to zero at land ice points.'

          ELSEIF ( iivar == iposSathh ) THEN
            CALL varValue( levels,summary,sathh,varFormat='f7.3',varName=varName(iivar) )

          ELSEIF ( iivar == iposSMsat ) THEN
            CALL varValue( levels,summary,smvcst,varFormat='f6.3',varName=varName(iivar) )

          ELSEIF ( iivar == iposSMcrit ) THEN
            CALL varValue( levels,summary,smvccl,varFormat='f6.3',varName=varName(iivar) )

          ELSEIF ( iivar == iposSMwilt ) THEN
            CALL varValue( levels,summary,smvcwt,varFormat='f6.3',varName=varName(iivar) )

          ENDIF

        ENDIF  !  varUse

      ENDDO  !  iivar
    ENDDO  !  i

  ENDIF  !  echo

  END SUBROUTINE init_soil
!###############################################################################
!###############################################################################
!###############################################################################
