!   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 sed2d_init
!--------------------------------------------------------------------
! This subroutine allocates arrays, initalizes variables, computes
! boundary flags and matrix coefficients for the jcg solver. 
! In pre-processing mode, user-defined filters are applied on 
! bathymetry and outputs of various routines are written in sed.out         
!
! Adapted from sediment_v8.F90 (MORSELFE, L. Pinto)
!                                                                   
! Author: guillaume dodet (gdodet01@univ-lr.fr, gdodet@lnec.pt)     
! Date:   06/12/2012                                               
!
! History: 
! 01/2013 - G.Dodet:  Added log file open statement  
! 02/2013 - G.Dodet:  Included sed2d_read
! 03/2013 - G.Dodet:  - Added call to sed2d_check; 
!                     - Added d50.gr3 reading;
! 04/2013 - G.Dodet:  - Added element and side variables 
!                       initialization and allocation;
!                     - Added node-centered  volume control area used
!                       in extrema filter;
! 05/2013 - G.Dodet:  - Corrected a bug in vc_area computation 
!                       (division by nne!);
!                     - Added slope_cr.gr3 reading;
! 06/2013 - G.Dodet:  - Initialisation of bedforms (z0);
!                     - Initialisation of qint for morpho. time step;
! 07/2013 - G.Dodet:  - Initialisation of cflsed
! 03/2014 - T.Guerin: - Removed allocation for d50 and d90: now done
!                       in sed2d_read because of multi-class multi-
!                       layer mode
!                     - Allocation and initialisation of d50moy,
!                       d90moy, h2
! 10/2016 - T.Guerin: Modifications related to the merge of single-
!                     class and multi-class routines
! 04/2017 - T.Guerin: Added imorphogrid.gr3 reading
!--------------------------------------------------------------------

  USE schism_glbl, ONLY: area,dp,iegl,indel,iond_global,ipgl,isbnd,      &
                      &ne_global,nea,elnode,nne,nond_global,nope_global, &
                      &np_global,mnei_p,np,npa,nsa,rho0,rhosed,rkind,  &
                      &xnd,ynd,in_dir,out_dir,len_in_dir,len_out_dir
  USE schism_msgp, ONLY: exchange_p2d,myrank,parallel_abort
  USE sed2d_mod, ONLY: bc_flag,bc_val,bed_del_e,bed_delta,cflsed,    &
                      &dpdxy,dpdxy_e,dpdxy_s,Cdsed,d_class,          &
                      &idsed2d,ifilt,ipre_filt,ipre_flag,mcoefd,qav, &
                      &qb,qb_e,qb_s,qdt_e,qs,qs_e,qs_s,qtot,qtot_e,  &
                      &qtot_s,s,slope_cr,vc_area,z0_e,               &
                      &z0cr_e,z0mr_e,z0sw_e,z0wr_e,u2_tmp,v2_tmp,    &
                      &d50,d90,F_class,h2,h_inf,nb_class,nb_layer,   &
                      &imnp,imorpho
  USE sed2d_filter

  IMPLICIT NONE

!- Local variables --------------------------------------------------
  INTEGER :: i,ie,ip,itmp1,itmp2,j,istat,k,nd
  REAL(rkind) :: aux1,aux2,tmp,xtmp,ytmp
  REAL(rkind), DIMENSION(npa) :: dp_filt
  LOGICAL :: lexist
  CHARACTER(LEN=80) :: filename
!- User-defined parameters ------------------------------------------
  INTEGER, PARAMETER :: fileid = 31  !Temporary file handle
  CHARACTER(LEN=*),PARAMETER :: FMT1='(I0,2X,F14.6,2X,F14.6,2X,F11.3)'
!--------------------------------------------------------------------

!- Open log file for sed2d ------------------------------------------
! This file is used to store information generated by sed2d
!--------------------------------------------------------------------
  IF(myrank == 0) THEN
    OPEN(idsed2d,FILE=out_dir(1:len_out_dir)//'sed2d.out',STATUS='REPLACE')
    WRITE(idsed2d,*)'***********************************************'
    WRITE(idsed2d,*)'*              Sed2d log file                 *'
    WRITE(idsed2d,*)'***********************************************'
  ENDIF

!- Define constant --------------------------------------------------
  s = rhosed/rho0
  if(s<=1) call parallel_abort('SED2D_INIT: s<=1')

!- Allocate arrays --------------------------------------------------
  ALLOCATE(bc_flag(npa),bc_val(npa),bed_del_e(nea),bed_delta(npa),    &
           cflsed(npa),Cdsed(npa),dpdxy(npa,2),                       &
           dpdxy_e(nea,2),dpdxy_s(nsa,2),mcoefd(0:mnei_p,np),         &
           qav(npa,2),qb(npa,2),qb_e(nea,2),qb_s(nsa,2),qdt_e(nea,2), &
           qs(npa,2),qs_e(nea,2),qs_s(nsa,2),                         &
           qtot(npa,2),qtot_e(nea,2),qtot_s(nsa,2),slope_cr(npa),     &
           vc_area(npa),u2_tmp(npa),v2_tmp(npa),z0_e(nea),z0cr_e(nea),&
           z0mr_e(nea),z0sw_e(nea),z0wr_e(nea),imnp(npa),stat=istat)
  IF(istat/=0) CALL parallel_abort('sed2d_init: allocation failure')

!- Read sed2d input parameters --------------------------------------
  CALL sed2d_read
  IF(myrank == 0)  WRITE(16,*)'done reading sed2d input parameters'

!- Allocate arrays related to multi-class multi-layer mode ----------
  ALLOCATE(d50(npa,nb_layer),d90(npa,nb_layer),h2(npa))

  IF ((nb_class.GT.1).and.(nb_layer.GT.1)) THEN
    DO i=1,nb_layer
      d50(:,i) = d_class(1)**F_class(:,1,i)
      DO j=2,nb_class
        d50(:,i) = d50(:,i) * d_class(j)**F_class(:,j,i)
      ENDDO
      d90(:,i) = 2.5d0*d50(:,i) !to be improved
    ENDDO
  ELSE
    d50(:,:) = d_class(1)
    d90(:,:) = 2.5d0*d50(:,:) !to be improved
  ENDIF

  h2(:) = h_inf

!- Initialize variables ---------------------------------------------
  Cdsed = 0.0025 !for first step in schism_step
  bed_del_e = 0.d0
  cflsed = 0.d0
  qav = 0.d0
  qb = 0.d0
  qb_e = 0.d0
  qb_s = 0.d0
  qdt_e = 0.d0
  qs = 0.d0
  qs_e = 0.d0
  qs_s = 0.d0
  qtot = 0.d0
  qtot_e = 0.d0
  qtot_s = 0.d0
  slope_cr = 0.d0
  z0_e = 0.d0 !z0_e = d50(1,1)/12.d0
  z0cr_e = 0.d0
  z0mr_e = 0.d0
  z0sw_e = 0.d0
  z0wr_e = 0.d0
  u2_tmp = 0.d0
  v2_tmp = 0.d0
  imnp = 1.d0

!- Read spatially variable grain size input files -------------------
!  IF(d50(1).LT.0.d0)THEN
!    INQUIRE(FILE='d50.gr3',EXIST=lexist)
!    IF(lexist) THEN
!      OPEN(fileid,FILE=in_dir(1:len_in_dir)//'d50.gr3',STATUS='OLD')
!      IF(myrank.EQ.0)WRITE(idsed2d,*)'d50 are read in d50.gr3'
!      READ(fileid,*)
!      READ(fileid,*)itmp1,itmp2
!      IF(itmp2/=np_global)CALL parallel_abort('sed2d: check d50.gr3')
!      DO i = 1,np_global
!         READ(fileid,*)j,xtmp,ytmp,tmp
!         IF(tmp<0) CALL parallel_abort('sed2d: d50 <0 !')
!         IF(ipgl(i)%rank == myrank) d50(ipgl(i)%id) = tmp
!         IF(ipgl(i)%rank == myrank) d90(ipgl(i)%id) = 2.5d0*tmp
!      ENDDO !np_global
!      CLOSE(fileid)
!    ELSE
!      IF(myrank.EQ.0)CALL parallel_abort('sed2d: d50.gr3 not found')
!    ENDIF !lexist
!    IF(myrank == 0)  WRITE(16,*)'done reading d50.gr3'
!  ENDIF !d50<0

!- Read file with critical slopes -----------------------------------
  IF(ifilt >= 2 .OR. ipre_filt == 2) THEN
    INQUIRE(FILE=in_dir(1:len_in_dir)//'slope_cr.gr3',EXIST=lexist)
    IF(lexist) THEN
      OPEN(fileid,FILE=in_dir(1:len_in_dir)//'slope_cr.gr3',STATUS='OLD')
      IF(myrank.EQ.0)WRITE(idsed2d,*)'Max. slopes are read in slope_c&
r.gr3'
      READ(fileid,*)
      READ(fileid,*)itmp1,itmp2
      IF(itmp2/=np_global)CALL parallel_abort('sed2d: Check slope_cr.&
gr3')
      DO i = 1,np_global
         READ(fileid,*)j,xtmp,ytmp,tmp
         IF(tmp<0) CALL parallel_abort('sed2d: max. slope <0 !')
         IF(ipgl(i)%rank == myrank) slope_cr(ipgl(i)%id) = tmp
      ENDDO !np_global
      CLOSE(fileid)
    ELSE
      IF(myrank.EQ.0)WRITE(idsed2d,*)'Default slopes values are used'
    ENDIF !lexist
  ENDIF !ifilt

!- Read morphological ramp file ------------------------------------
  IF (imorpho==1) THEN
    INQUIRE(FILE=in_dir(1:len_in_dir)//'imorphogrid.gr3',EXIST=lexist)
    IF (lexist) THEN
      OPEN(fileid,FILE=in_dir(1:len_in_dir)//'imorphogrid.gr3',STATUS='OLD')
      IF(myrank.EQ.0)WRITE(idsed2d,*)'Morphological ramp values are &
&read in imorphogrid.gr3'
      READ(fileid,*)
      READ(fileid,*)itmp1,itmp2
      IF(itmp2/=np_global)CALL parallel_abort('sed2d: Check np_global &
&in imorphogrid.gr3')
      DO i = 1,np_global
         READ(fileid,*)j,xtmp,ytmp,tmp
         IF(tmp<0.or.tmp>1) CALL parallel_abort('sed2d: local imorpho <0 or >1!')
         IF(ipgl(i)%rank == myrank) imnp(ipgl(i)%id) = tmp
      ENDDO !np_global
      CLOSE(fileid)
    ELSE
      IF(myrank.EQ.0)WRITE(idsed2d,*)'no morphological ramp'
    ENDIF !lexist
  ENDIF !imorpho

!- Compute "2D control volume" area at each node --------------------
  vc_area = 0.d0
  DO i = 1,np
     DO j = 1,nne(i)
        vc_area(i) = vc_area(i)+area(indel(j,i))/3.d0
     ENDDO
  ENDDO
  CALL exchange_p2d(vc_area)

!- Compute matrix coefficients --------------------------------------
! These variables are used in sed2d_morpho
!--------------------------------------------------------------------
!Error: YJZ - not working for quads
  mcoefd = 0.d0
  aux1 = 22.d0/108.d0
  aux2 = 7.d0/108.d0
  DO i = 1,np !residents
     DO j = 1,nne(i)
        ie = indel(j,i)
        mcoefd(0,i) = mcoefd(0,i)+area(ie)*aux1 !diagonal
        mcoefd(j,i) = mcoefd(j,i)+area(ie)*aux2
        IF(isbnd(1,i) == 0.and.j == nne(i)) THEN !internal ball
          mcoefd(1,i) = mcoefd(1,i)+area(ie)*aux2
        ELSE
          mcoefd(j+1,i) = mcoefd(j+1,i)+area(ie)*aux2
        ENDIF
     ENDDO !nne
  ENDDO !np

!- Compute b.c. values and b.c. flags -------------------------------
! These variables are used in sed2d_morpho and in sed2d_filter
!--------------------------------------------------------------------
  bc_val = -9999 ! flags
  DO i = 1,nope_global
     DO j = 1,nond_global(i)
        nd = iond_global(i,j) !global
        IF(ipgl(nd)%rank == myrank) THEN
          ip = ipgl(nd)%id
          bc_val(ip) = 0 !no flux b.c. at the moment
        ENDIF !ipgl(nd)%rank==myrank
     ENDDO !nond_global
  ENDDO !nope_global

!- Compute b.c. flag for all nodes 
  bc_flag = .false.
  DO i = 1,np
     IF(bc_val(i) > -9998) bc_flag(i) = .true.
  ENDDO
  IF(myrank == 0) WRITE(16,*)'done initializing sed2d'
  IF(ipre_flag == 0) RETURN

!- The rest is pre-proc - Apply user-defined filter to bathymetry --------------------------
! If ipre_flag = 1, filters are applied to hgrid.gr3 and the new 
! bathymetry is written in hgrid_filt#.gr3 (where # depends on the 
! chosen filter, ie. ipre_filt)
!--------------------------------------------------------------------
  IF(myrank > 0) CALL parallel_abort('IPRE only works on 1 proc')

  WRITE(idsed2d,*)'----------------- FILTERING ---------------------'
  SELECT CASE(ipre_filt)
    CASE(0) !No filter
      dp_filt = dp
    CASE(1)
      CALL sed2d_filter_extrema(dp,dp_filt)
    CASE(2)
      CALL sed2d_filter_slope(dp,dp_filt)
    CASE(3) 
      CALL sed2d_filter_diffu(dp,dp_filt,npa)
  END SELECT
  dp = dp_filt
  
!- Write filtered bathymetry in hgrid_filt*.gr3 ---------------------
  WRITE(filename,'(A10,I1.1,A4)')'hgrid_filt',ipre_filt,'.gr3'

  OPEN(fileid,FILE=out_dir(1:len_out_dir)//TRIM(filename),STATUS='REPLACE')
  WRITE(fileid,*)TRIM(filename)
  WRITE(fileid,*)ne_global,np_global
  DO i = 1,np_global
     WRITE(fileid,FMT1)i,xnd(i),ynd(i),dp(i)
  ENDDO
  DO i = 1,ne_global
     WRITE(fileid,*)i,3,(elnode(k,i),k=1,3)
  ENDDO
  CLOSE(fileid)

  WRITE(idsed2d,*)'Filtered bathymetry written in: ',TRIM(filename)
  WRITE(idsed2d,*)'-------------------------------------------------'
  WRITE(idsed2d,*)''

  WRITE(16,*)'done filtering bathymetry in sed2d'

!- Check values in sed2d routines -----------------------------------
  CALL sed2d_check

  WRITE(16,*)'done checking values in sed2d'
!--------------------------------------------------------------------
  call mpi_finalize(itmp1)
  STOP

END SUBROUTINE sed2d_init