!   Copyright 2014 College of William and Mary
!
!   Licensed under the Apache License, Version 2.0 (the "License");
!   you may not use this file except in compliance with the License.
!   You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
!   Unless required by applicable law or agreed to in writing, software
!   distributed under the License is distributed on an "AS IS" BASIS,
!   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!   See the License for the specific language governing permissions and
!   limitations under the License.

      SUBROUTINE read_sed_input
!--------------------------------------------------------------------!
! This subroutine reads sediment model inputs                        !
!                                                                    !
! Author: Ligia Pinto                                                !
! Date: xx/09/2007                                                   !
!                                                                    !
! History: 2012/12 - F.Ganthy : form homogenisation of sediments     !
!          routines                                                  !
!          2013/01 - F.Ganthy : Introduction of switchs for roughness!
!          predictor                                                 !
!          2013/01 - F.Ganthy : Implementation of avalanching        !
!          2013/05 - F.Ganthy : Updates related to ripple predictor  !
!          2013/05 - F.Ganthy : Updates to the ripple predictor:     !
!                               - Changes on the total bedform       !
!                                 roughness computation              !
!                               - Add wave-ripple computation from   !
!                                 Nielsen (1992)                     !
!          2013/05 - F.Ganthy : Add different sediment behavior:     !
!                                - MUD-like or SAND-like             !
!          2013/05 - F.Ganthy : Re-introduction of the number of bed !
!                               sediment layers within sediment.in   !
!          2020/02 - B.Mengual                                       !
!                  * Implementation of a method to compute and update!
!                    a porosity representative of non-cohesive       !
!                    matrix (poro_option,poro_cst,Awooster,Bwooster) !
!                  * Implementation of several methods to compute    !
!                    the current-induced bottom shear stress         !
!                    (tau_option,tau_max,zstress)                    !
!                  * Wave asymmetry computations: restructuration    !
!                    and introduction of filtering possibilities     !
!                    (iasym,w_asym_max,elfrink_filter,ech_uorb)      !
!                  * Implementation of bedload flux caused by        !
!                    wave acceleration skewness                      !
!                    (bedload_acc,kacc_hoe,                          !
!                    kacc_dub,thresh_acc_opt,acrit)                  !
!                  * Add the possibility to filter and limit         !
!                    bedload fluxes (bedload_filter,bedload_limiter, !
!                    bedload_acc_filter)
!                  * Add the possibility to use different methods    !
!                    to solve the Exner equation in bedchange_bedload!
!                    (IMETH_BED_EVOL)                                !
!                  * Introduction of a maximum active layer thickness!
!                    (actv_max)                                      !
!                  * Morphodynamics: common morph_fac over all       !
!                    sediment classes (morph_fac)                    !
!                  * Initialisation of the vertical discretisation of!
!                    sediment layers (sedlay_ini_opt,toplay_inithick)!
!                  * Bed slope effects for bedload fluxes            !
!                    (slope_formulation=4 with alpha_bs,alpha_bn)    !
!--------------------------------------------------------------------!

      USE schism_glbl, ONLY: rkind,ntrs,iwsett,wsett,irange_tr,ddensed, &
   &grav,rho0,errmsg,in_dir,out_dir,len_in_dir,len_out_dir
      USE schism_msgp, ONLY: myrank,parallel_abort
      USE sed_mod      

      IMPLICIT NONE
      SAVE  

!- Local variables --------------------------------------------------!

      CHARACTER(len=100) :: var1, var2
      INTEGER :: i, j, ierror

      INTEGER :: line, len_str, loc, loc2, lstr_tmp
      CHARACTER(len=90) :: line_str,str_tmp,str_tmp2

      namelist /SED_CORE/Sd50,Erate
      namelist /SED_OPT/Wsed,tau_ce,Srho,iSedtype,newlayer_thick,bedload_coeff, &
     &sed_debug,Cdb_min,Cdb_max,actv_max,poro_option,porosity,Awooster,Bwooster, &
     &sedlay_ini_opt,toplay_inithick,bdldiffu,Nbed,bedload,bedload_filter, &
     &bedload_limiter,suspended_load,slope_formulation,alpha_bs,alpha_bn, &
     &ised_bc_bot,sed_morph,sed_morph_time,morph_fac,drag_formulation, &
     &ddensed,comp_ws,comp_tauce,bedforms_rough,iwave_ripple,irough_bdld, &
     &slope_avalanching,dry_slope_cr,wet_slope_cr,bedmass_filter,bedmass_threshold, &
     &ised_dump,ierosion,alphd,refht,Tbp,im_pick_up,imeth_bed_evol,iasym, &
     &w_asym_max,elfrink_filter,ech_uorb,bedload_acc,bedload_acc_filter, &
     &kacc_hoe,kacc_dub,thresh_acc_opt,acrit,tau_option,tau_max,zstress

!- Start Statement --------------------------------------------------!
      !defined # of classes (in sed_mod) for all routines
      ntr_l=ntrs(5) 

      IF(myrank==0) write(16,*)'Entering read_sed_input routine'
      IF(myrank==0) WRITE(16,*)'reading sediment.nml'
      IF(myrank==0) WRITE(16,'(A,A,I2,A)')'##### Number of Tracers',  &
      &             ' / Sediment Classes required in sediment.in: ', &
      &             ntr_l,' #####'


      ALLOCATE(Sd50(ntr_l),Srho(ntr_l),Wsed(ntr_l),Erate(ntr_l),tau_ce(ntr_l), &
     &iSedtype(ntr_l),stat=i)
      IF(i/=0) CALL parallel_abort('main: Sd50 allocation failure')
       
!--------------------------------------------------------------------
! - Scan sediment.in.
!--------------------------------------------------------------------

      !Take some old parameters out of sediment.in
      g=grav
      rhom=rho0
      vonKar=0.4

      !Init default parameters. Not init'ed are (which will be computed below): Wsed,tau_ce
      Sd50=-huge(1.d0); Erate=-huge(1.d0)

      Srho(:)=2650.d0; iSedtype(:)=1; 
      newlayer_thick=100.d0; bedload_coeff=1
      sed_debug=0; Cdb_min=1.d-6; Cdb_max=1.d-2; actv_max=1.d0 
      poro_option=1; porosity=0.4d0; Awooster=0.42d0; Bwooster=-0.458d0
      sedlay_ini_opt=0; toplay_inithick=10.0d-2; bdldiffu=0.5d0; Nbed=1
      bedload=1; bedload_filter=0; bedload_limiter=0; suspended_load=1
      slope_formulation=4; alpha_bs=1.0d0; alpha_bn=1.5d0; ised_bc_bot=1
      sed_morph=0; sed_morph_time=1.d0; morph_fac=1.d0; drag_formulation=1
      ddensed=0; comp_ws=1; comp_tauce=1; bedforms_rough=0; iwave_ripple=1;
      irough_bdld=1; slope_avalanching=1; dry_slope_cr=1.d0; wet_slope_cr=0.3d0
      bedmass_filter=0; bedmass_threshold=0.025d0; ised_dump=0; ierosion=0
      alphd=1.d0; refht=0.75d0; Tbp=1.d2; im_pick_up=4 
      imeth_bed_evol=2; iasym=0; w_asym_max=0.4d0; elfrink_filter=0;
      ech_uorb=200; bedload_acc=0; bedload_acc_filter=0; kacc_hoe=1.4d-4
      kacc_dub=0.631d-4; thresh_acc_opt=2; acrit=0.2d0; tau_option=1
      tau_max=10.d0; zstress=0.2d0

      OPEN(32,file=in_dir(1:len_in_dir)//'sediment.nml',delim='apostrophe',status='old')
      REWIND(32) !go to beginning of file
      read(32,nml=SED_CORE)
      read(32,nml=SED_OPT)
      CLOSE(32)

!     Check core parameters
      if(minval(Sd50)<0.d0.or.minval(Erate)<0.d0) then
        write(errmsg,*)'READ_SED, core pars. not read in:',Sd50,Erate
        call parallel_abort(errmsg)
      endif

!     Check optional parameters (incomplete)
      if(ised_bc_bot<0.or.ised_dump<0.or.ierosion<0) then
        write(errmsg,*)'READ_SED, some pars. not rite:',ised_bc_bot,ised_dump,ierosion
        call parallel_abort(errmsg)
      endif 

      IF ((tau_option .EQ. 0) .OR. (tau_option .GT. 3)) THEN
        WRITE(errmsg,*)'READ_SED: Invalid tau_option (from 1 to 3)'
        CALL parallel_abort(errmsg)
      END IF

!      line=0
!      ! - PS: start reading
!      DO
!        ! - PS: read line an break on '!'
!        line=line+1
!        READ(5,'(a)',END=99)line_str
!        line_str=ADJUSTL(line_str)
!        len_str=LEN_TRIM(line_str)
!        IF(len_str==0.OR.line_str(1:1)=='!') CYCLE

!        ! - PS: locate '==' and '!'
!        loc=INDEX(line_str,'==')
!        loc2=INDEX(line_str,'!')
!        IF(loc2/=0.AND.loc2-1<loc+1)THEN
!          CALL parallel_abort('READ_PARAM: ! before =')
!        ENDIF

!        ! - PS: get name of the variable
!        str_tmp=''
!        str_tmp(1:loc-1)=line_str(1:loc-1)
!        str_tmp=TRIM(str_tmp)
!        lstr_tmp=LEN_TRIM(str_tmp)

!        ! PS: get the value
!        !if(loc2/=0) then
!        !  str_tmp2=line_str(loc+2:loc2-1)
!        !else
!        !  str_tmp2=line_str(loc+2:len_str)
!        !endif

!        !str_tmp2=adjustl(str_tmp2)
!        !str_tmp2=trim(str_tmp2)

!        ! - PS: switch between variables and set values
!        SELECT CASE(str_tmp)

!          CASE('NEWLAYER_THICK')
!            READ(line_str,*) var1, var2, newlayer_thick

!          CASE('BEDLOAD_COEFF')
!            READ(line_str,*) var1, var2, bedload_coeff

!          CASE('SAND_SD50')
!            READ(line_str,*) var1, var2, (Sd50(i), i=1,ntr_l)

!          CASE('SAND_SRHO')
!            READ(line_str,*) var1, var2, (Srho(i), i=1,ntr_l)

!          CASE('SAND_WSED')
!            READ(line_str,*) var1, var2, (Wsed(i), i=1,ntr_l)

!          CASE('SAND_ERATE')
!            READ(line_str,*) var1, var2, (Erate(i), i=1,ntr_l)

!          CASE('SAND_TAU_CE')
!            READ(line_str,*) var1, var2, (tau_ce(i), i=1,ntr_l)

!          CASE('SED_TYPE')
!            READ(line_str,*) var1, var2, (iSedtype(i), i=1,ntr_l)

!          CASE('sed_debug')
!            READ(line_str,*) var1, var2, sed_debug

!          CASE('g')
!            READ(line_str,*) var1, var2, g

!          CASE('vonKar')
!            READ(line_str,*) var1, var2, vonKar

!          CASE('Cdb_min')
!            READ(line_str,*) var1, var2, Cdb_min

!          CASE('Cdb_max')
!            READ(line_str,*) var1, var2, Cdb_max

!          CASE('actv_max')
!            READ(line_str,*) var1, var2, actv_max

!          CASE('rhom')
!            READ(line_str,*) var1, var2, rhom

!          CASE('poro_option')
!            READ(line_str,*) var1, var2, poro_option

!          CASE('poro_cst')
!            READ(line_str,*) var1, var2, porosity

!          CASE('Awooster')
!            READ(line_str,*) var1, var2, Awooster

!          CASE('Bwooster')
!            READ(line_str,*) var1, var2, Bwooster

!          CASE('sedlay_ini_opt')
!            READ(line_str,*) var1, var2, sedlay_ini_opt

!          CASE('toplay_inithick')
!            READ(line_str,*) var1, var2, toplay_inithick

!          CASE('bdldiffu')
!            READ(line_str,*) var1, var2, bdldiffu

!          CASE('Nbed')
!            READ(line_str,*) var1, var2, Nbed

!          CASE('bedload')
!            READ(line_str,*) var1, var2, bedload

!          CASE('bedload_filter')
!            READ(line_str,*) var1, var2, bedload_filter

!          CASE('bedload_limiter')
!            READ(line_str,*) var1, var2, bedload_limiter

!          CASE('suspended_load')
!            READ(line_str,*) var1, var2, suspended_load

!          CASE('slope_formulation')
!            READ(line_str,*) var1, var2, slope_formulation

!          CASE('alpha_bs')
!            READ(line_str,*) var1, var2, alpha_bs

!          CASE('alpha_bn')
!            READ(line_str,*) var1, var2, alpha_bn

!          CASE('ised_bc_bot')
!            READ(line_str,*) var1, var2, ised_bc_bot

!          CASE('sed_morph')
!            READ(line_str,*) var1, var2, sed_morph

!          !CASE('depo_scale')
!          !  READ(line_str,*) var1, var2, depo_scale

!          CASE('sed_morph_time')
!            READ(line_str,*) var1, var2, sed_morph_time !in days

!          CASE('sed_morph_fac')
!            READ(line_str,*) var1, var2, morph_fac ! morphological factor [-]

!          CASE('drag_formulation')
!            READ(line_str,*) var1, var2, drag_formulation

!          CASE('ddensed')
!            READ(line_str,*) var1, var2, ddensed

!          CASE('comp_ws')
!            READ(line_str,*) var1, var2, comp_ws

!          CASE('comp_tauce')
!            READ(line_str,*) var1, var2, comp_tauce

!          CASE('bedforms_rough')
!            READ(line_str,*) var1, var2, bedforms_rough

!          CASE('iwave_ripple')
!            READ(line_str,*) var1, var2, iwave_ripple

!          CASE('irough_bdld')
!            READ(line_str,*) var1, var2, irough_bdld

!          CASE('slope_avalanching')
!            READ(line_str,*) var1, var2, slope_avalanching

!          CASE('dry_slope_cr')
!            READ(line_str,*) var1, var2, dry_slope_cr
! 
!          CASE('wet_slope_cr')
!            READ(line_str,*) var1, var2, wet_slope_cr

!          CASE('bedmass_filter')
!            READ(line_str,*) var1, var2, bedmass_filter

!          CASE('bedmass_threshold')
!            READ(line_str,*) var1, var2, bedmass_threshold  

!          CASE('relath')
!            READ(line_str,*) var1, var2, relath

!          CASE('ised_dump')
!            READ(line_str,*) var1, var2, ised_dump

!          CASE('ierosion')
!            READ(line_str,*) var1, var2, ierosion

!          CASE('alphd') !1120:+alphd,refht,Tbp,im_pick_up
!            READ(line_str,*) var1, var2, alphd

!          CASE('refht')
!            READ(line_str,*) var1, var2, refht

!          CASE('Tbp')
!            READ(line_str,*) var1, var2, Tbp

!          CASE('im_pick_up')
!            READ(line_str,*) var1, var2, im_pick_up

!          CASE('IMETH_BED_EVOL')
!            READ(line_str,*) var1, var2, imeth_bed_evol

!          CASE('iasym')
!            READ(line_str,*) var1, var2, iasym

!          CASE('w_asym_max')
!            READ(line_str,*) var1, var2, w_asym_max

!          CASE('elfrink_filter')
!            READ(line_str,*) var1, var2, elfrink_filter

!          CASE('ech_uorb')
!            READ(line_str,*) var1, var2, ech_uorb

!          CASE('bedload_acc')
!            READ(line_str,*) var1, var2, bedload_acc

!          CASE('bedload_acc_filter')
!            READ(line_str,*) var1, var2, bedload_acc_filter

!          CASE('kacc_hoe')
!            READ(line_str,*) var1, var2, kacc_hoe

!          CASE('kacc_dub')
!            READ(line_str,*) var1, var2, kacc_dub

!          CASE('thresh_acc_opt')
!            READ(line_str,*) var1, var2, thresh_acc_opt

!          CASE('acrit')
!            READ(line_str,*) var1, var2, acrit

!          CASE('tau_option')
!            READ(line_str,*) var1, var2, tau_option

!          CASE('tau_max')
!            READ(line_str,*) var1, var2, tau_max

!          CASE('zstress')
!            READ(line_str,*) var1, var2, zstress

!        END SELECT

!      ENDDO !scan sediment.in
!99    CLOSE(5)


!      g = g*1.0d0
!      vonKar = vonKar*1.0d0
!      Cdb_min = Cdb_min*1.0d0
!      Cdb_max = Cdb_max*1.0d0
!      rhom = rhom*1.0d0


!---------------------------------------------------------------------
! - Scale input parameters
!
! Particel settling velocity (Wsed) input is in [mm/s]
! Median sediment grain diameter (Sd50) input is in [mm]
! Critical shear for erosion and deposition (Tau_ce) input is in [N/m2]
!---------------------------------------------------------------------

      ! Conversion of threshold from mm to m
      bedmass_threshold = bedmass_threshold/1000.0d0

      DO i=1,ntr_l !ntracers
        Sd50(i) = Sd50(i)*0.001d0 !to [m]
        IF(comp_ws.EQ.1.AND.iSedtype(i).EQ.1) THEN
          ! Sed type is SAND
          CALL sed_settleveloc(i)
        ELSE
          ! Sed type is other (MUD or GRAVEL)
          Wsed(i) = Wsed(i)*0.001d0 !to m/s
        ENDIF
        IF(comp_tauce.EQ.1.AND.iSedtype(i).EQ.1)THEN
          ! Sed type is SAND
          CALL sed_taucrit(i)
        ENDIF

        ! tau_ce already read in for other Sed type (MUD or GRAVEL); scale to
        ! get m^2/s/s
        tau_ce(i) = tau_ce(i)/rhom
      ENDDO !i

      !Update global array for settling vel.
      do i=1,ntr_l
        wsett(i-1+irange_tr(1,5),:,:)=Wsed(i)
        iwsett(i-1+irange_tr(1,5))=1
      enddo !i

!---------------------------------------------------------------------
! - Screen output of sediment-model parameters
!---------------------------------------------------------------------
      if(myrank==0) then
        open(32,file=out_dir(1:len_out_dir)//'sediment.out.nml',status='replace')
        write(32,nml=SED_CORE)
        write(32,nml=SED_OPT)
        close(32)

        WRITE(16,*)
        WRITE(16,*) 'Sediment Parameters:'
        WRITE(16,*) '--------------------'
        WRITE(16,*) 'sed_debug: ', sed_debug
        WRITE(16,*) 'Nbed: ',Nbed
        IF (poro_option .EQ. 1) THEN
          WRITE(16,*) 'Constant porosity: ',porosity
        ELSE
          WRITE(16,*) 'Variable porosity poro_option > 1'
        END IF
        IF (sed_morph .EQ. 1) THEN
          WRITE(16,*) 'Accounting for morphodynamics'
          WRITE(16,*) ' > sed_morph_time (ramp in days)=',sed_morph_time
          WRITE(16,*) ' > morphological factor: ',morph_fac
        END IF
!        WRITE(16,*) 'bedthick_overall(1): ',bedthick_overall(1)
        WRITE(16,*) 'New layer thickness',newlayer_thick
        WRITE(16,*) '----------------------------------------------'
        WRITE(16,*) 'Sed type  ;  Sd50 [m]  ;  Srho [kg/m3]  ;',     &
        &          '  Wsed [m/s]  ;  Erate ;  tau_ce', &
        &          ' [Pa]  ;  ierosion '
        DO i = 1,ntr_l !ntracers
          WRITE(16,*) iSedtype(i),Sd50(i),Srho(i),Wsed(i),Erate(i),   &
          &           tau_ce(i)*rhom,ierosion
        ENDDO
        WRITE(16,*)
      ENDIF

      IF(myrank==0) WRITE(16,*)'Leaving read_sed_input routine'
!--------------------------------------------------------------------!
       END SUBROUTINE read_sed_input