 !! ---------------------------------------------------------------------------
 !! ---------------------------------------------------------------------------
 !! ---------------------------------------------------------------------------
 !!
 !!    Copyright (c) 2018-2023, Universita' di Padova, Manuele Faccenda
 !!    All rights reserved.
 !!
 !!    This software package was developed at:
 !!
 !!         Dipartimento di Geoscienze
 !!         Universita' di Padova, Padova         
 !!         via Gradenigo 6,            
 !!         35131 Padova, Italy 
 !!
 !!    project:    ECOMAN
 !!    funded by:  ERC StG 758199 - NEWTON
 !!
 !!    ECOMAN is free software package: you can redistribute it and/or modify
 !!    it under the terms of the GNU General Public License as published
 !!    by the Free Software Foundation, version 3 of the License.
 !!
 !!    ECOMAN is distributed WITHOUT ANY WARRANTY; without even the implied
 !!    warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 !!    See the GNU General Public License for more details.
 !!
 !!    You should have received a copy of the GNU General Public License
 !!    along with ECOMAN. If not, see <http://www.gnu.org/licenses/>.
 !!
 !!
 !!    Contact:
 !!        Manuele Faccenda    [manuele.faccenda@unipd.it]
 !!        Brandon VanderBeek  [brandon.vanderbeek@unipd.it]
 !!
 !!
 !!    Main development team:
 !!        Manuele Faccenda    [manuele.faccenda@unipd.it]
 !!        Brandon VanderBeek  [brandon.vanderbeek@unipd.it]
 !!        Albert de Montserrat Navarro
 !!        Jianfeng Yang   
 !!
 !! ---------------------------------------------------------------------------
 !! ---------------------------------------------------------------------------
 !! ---------------------------------------------------------------------------

   PROGRAM DREX_S  

   USE comvar
   USE hdf5    

   IMPLICIT NONE

   INTEGER :: i,j,m,n,nx(3),ti(1),anisnum,nsave

   DOUBLE PRECISION, DIMENSION(3) :: evals,c2
   DOUBLE PRECISION, DIMENSION (3,3) ::LSij,evects,fseacs,acs1,acs2,Rotm,ee,Rotm1,lrot
   DOUBLE PRECISION :: fractvoigt0,fractdisl,phi1,theta,phi2,a0,w3
   DOUBLE PRECISION :: pphi1,ttheta,pphi2
   DOUBLE PRECISION, DIMENSION(6,6) :: Voigt,Reuss,Mixed,Cstilwe
   DOUBLE PRECISION, DIMENSION(1000,6,6) :: Cdem
   DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: phi_a
   DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: f1_ol,f1_opx
   ! average orientation of a-axis
   DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: perc_a,perc_anis,perc_hexa,perc_tetra,perc_ortho,perc_mono,perc_tri
   ! percentage of S wave anisotropy

   CHARACTER (3) :: dt_str3
   CHARACTER (4) :: dt_str4
   CHARACTER (100) :: fname,str

   INTEGER(HID_T)  :: file_id, group_id, dataspace_id, dataset_id, attr_id, dcpl,memtype ! Handles
   INTEGER     ::   error  ! Error flag

!!! initialization

   CALL init0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!  LPO  and FSE calculation  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!! Constant timestep   
   dt = 1d-2/epsnot(1)

   write(*,"(a,f8.2)"),' Timestep = ',dt
   write(*,*)

   fractdisl = fractdislrock(rocktype(1))

   CALL strain(1,1,fractdisl)

   STOP

   END PROGRAM DREX_S

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! subroutine INIT0, initialization of variables, arrays and parameters   !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   SUBROUTINE init0

   USE comvar
   USE hdf5  

   IMPLICIT NONE

   INTEGER :: gi,j,i1,i3,i,j1,j2,j3,nrot ! loop counters
   INTEGER :: iph,ith,ips,nbox

   DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: ran0
   DOUBLE PRECISION :: dph,dcosth,dps,ph,costh,ps,th
   DOUBLE PRECISION :: p21,p31,p41
   ! matrix of random numbers used to generate initial random LPO

   DOUBLE PRECISION :: phi1,theta,phi2
   ! eulerian angles

   DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xe1,xe2,xe3
   ! matrixes of initial random eulerian angles

   DOUBLE PRECISION, DIMENSION(3) :: evals
   DOUBLE PRECISION, DIMENSION(3,3) :: evects,ee
   ! eigen values and vectors in jacobi

   DOUBLE PRECISION phi
   DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: phix0,phiy0,phiz0,phiz1
   DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: Rhodum,Vpdum,Vsdum

   INTEGER(HID_T)  :: file_id, group_id, dataspace_id, dataset_id, attr_id, dcpl,memtype ! Handles
   INTEGER     ::   error  ! Error flag

   ! to fix random distribution
   INTEGER, ALLOCATABLE :: seed(:)

   l = 0d0 ; e = 0d0

   pi = 3.141592653589793238462643383279
   deg2rad = pi/180.0

!!! Name of input files from file input_fabric0.dat

   CALL read_input_file

!!! initial size
   size3 = 5
   size = 5!size3**3

   Xol = Xol /100.0

!!! strain rate tensor
   e(1,1,1) = l(1,1,1) ; e(1,3,3) = l(1,3,3) ; e(1,2,2) = l(1,2,2)
   e(1,3,1) = (l(1,1,3)+l(1,3,1))/2d0 ; e(1,1,3) = e(1,3,1)
   e(1,1,2) = (l(1,2,1)+l(1,1,2))/2d0 ; e(1,2,1) = e(1,1,2)
   e(1,2,3) = (l(1,3,2)+l(1,2,3))/2d0 ; e(1,3,2) = e(1,2,3)

!! reference strain rate
   ee = e(1,:,:)
   CALL DSYEVQ3(ee,evects,evals)
   epsnot(1) = MAXVAL(ABS(evals))

   write(*,"(a)"),' VELOCITY GRADIENT TENSOR Dij'
   write(*,*)
   write(*,"(3f8.2)"),(l(1,1,:))
   write(*,"(3f8.2)"),(l(1,2,:))
   write(*,"(3f8.2)"),(l(1,3,:))
   write(*,*)
   write(*,"(a)"),' STRAIN RATE tensor'
   write(*,*)
   write(*,"(3f8.2)"),(e(1,1,:))
   write(*,"(3f8.2)"),(e(1,2,:))
   write(*,"(3f8.2)"),(e(1,3,:))
   write(*,*)
   write(*,"(a)"),'--------------------------------------------------------'
   write(*,*)

   max_strain = 0d0
        
!!! Initial deformation gradient tensor
   Fij = 0d0 ; Fij(1,1,:) = 1d0 ; Fij(2,2,:) = 1d0 ; Fij(3,3,:) = 1d0

!!! tensor \epsilon_{ijk}

   alt=0d0
   alt(1,2,3) = 1d0 ; alt(2,3,1) = 1d0 ; alt(3,1,2) = 1d0
   alt(1,3,2) = -1d0 ; alt(2,1,3) = -1d0 ; alt(3,2,1) = -1d0

!!! tensor \delta_{ij}

   del=0d0
   del(1,1) = 1d0 ; del(2,2) = 1d0 ; del(3,3) = 1d0

!!! tensors of indices

   ijkl(1,1) = 1 ; ijkl(1,2) = 6 ; ijkl(1,3) = 5
   ijkl(2,1) = 6 ; ijkl(2,2) = 2 ; ijkl(2,3) = 4
   ijkl(3,1) = 5 ; ijkl(3,2) = 4 ; ijkl(3,3) = 3

   l1(1) = 1 ; l1(2) = 2 ; l1(3) = 3
   l1(4) = 2 ; l1(5) = 3 ; l1(6) = 1
   l2(1) = 1 ; l2(2) = 2 ; l2(3) = 3
   l2(4) = 3 ; l2(5) = 1 ; l2(6) = 2

   mandel_scale = 1d0

   DO j1=1,6 ; DO j2=1,6
      IF(j1 .GT. 3 .AND. j2 .GT. 3) THEN
         mandel_scale(j1,j2) = 2
         !mandel_scale(j1,j2) = 4
      ELSE IF(j1 .GT. 3 .OR. j2 .GT. 3) THEN
         mandel_scale(j1,j2) = 2**0.5d0
         !mandel_scale(j1,j2) = 2
      END IF
   END DO ; END DO

!!! Loading stiffness tensors (GPa)

   !CALL elastic_database(S0,dS0dp,dS0dp2,dS0dt,dS0dpt,dS0dp5,dS0dt5)

!!! allocation of the dimensions of the arrays

   ALLOCATE(odf(1,size))
   ALLOCATE(odf_ens(1,size))

   ALLOCATE(ran0(3*size))

   ALLOCATE(acs(size,3,3,1),acs0(size,3,3))
   ALLOCATE(acs_ens(size,3,3,1))

   
   IF(sbfmod == 0) THEN

!!! initialization of orientations - uniformally random distribution
!!! Rmq cos(theta) used to sample the metric Eulerian space

   ALLOCATE(xe1(size),xe2(size),xe3(size))

   CALL RANDOM_SEED(size=i)
   ALLOCATE(seed(i))
   DO j1 = 1,i
      seed(j1) = 0
   END DO
   CALL RANDOM_SEED(put=seed)
   CALL RANDOM_NUMBER(ran0)
   DEALLOCATE(seed)

   i = 1

   DO j1 =1, size3 ; DO j2 =1, size3 ; DO j3 =1, size3
      xe1(i) = (REAL(j1)-ran0(i))/REAL(size3)*ACOS(-1d0)
      xe2(i) = ACOS(-1d0 + (REAL(j2)-ran0(size+i))/REAL(size3)*2d0)
      xe3(i) = (REAL(j3)-ran0(i+2*size))/REAL(size3)*ACOS(-1d0)
      i = i + 1
   END DO ; END DO ; END DO

   DO i = 1 , size

      phi1 = xe1(i) ; theta = xe2(i) ; phi2 = xe3(i)

!!! Direction cosine matrix
!!! acs(k,j) = cosine of the angle between the kth crystallographic axes and the jth external axes 

      acs0(i,1,1)=COS(phi2)*COS(phi1)-COS(theta)*SIN(phi1)*SIN(phi2)
      acs0(i,1,2)=COS(phi2)*SIN(phi1)+COS(theta)*COS(phi1)*SIN(phi2)
      acs0(i,1,3)=SIN(phi2)*SIN(theta)

      acs0(i,2,1)=-SIN(phi2)*COS(phi1)-COS(theta)*SIN(phi1)*COS(phi2)
      acs0(i,2,2)=-SIN(phi2)*SIN(phi1)+COS(theta)*COS(phi1)*COS(phi2)
      acs0(i,2,3)=COS(phi2)*SIN(theta)

      acs0(i,3,1)=SIN(theta)*SIN(phi1)
      acs0(i,3,2)=-SIN(theta)*COS(phi1)
      acs0(i,3,3)=COS(theta)

   END DO
   
!!! Set random initial LPO and same grain size

   DO i = 1 , 1        

      odf(i,:) = 1d0/5d0!REAL(size3**3)
      odf_ens(i,:) = odf(i,:)
      !acs(:,:,:,i) = acs0
      !acs_ens(:,:,:,i) = acs0

   END DO

   DEALLOCATE(xe1,xe2,xe3,ran0)

   acs(1,:,:,1) = 0.5d0

   acs(2,1,1,1) = 0.1d0
   acs(2,1,2,1) = 0.2d0
   acs(2,1,3,1) = 0.3d0
   acs(2,2,1,1) = 0.4d0
   acs(2,2,2,1) = 0.5d0
   acs(2,2,3,1) = 0.6d0
   acs(2,3,1,1) = 0.7d0
   acs(2,3,2,1) = 0.8d0
   acs(2,3,3,1) = 0.9d0

   DO i = 3 , 5        
      acs(i,:,:,1) = acs(2,:,:,1)
   END DO

   RETURN

   ELSE

   nbox3 = 20
   nbox = nbox3
   nboxnum = nbox3**3
   rocktype(1) = 1
   ALLOCATE(acs0sbf(nboxnum,3,3))

   dph = pi/REAL(nbox)
   dcosth = 2.0/REAL(nbox)
   dps = pi/REAL(nbox)

   DO iph = 1, nbox
      ph = (iph - 0.5)*dph
   DO ith = 1, nbox
      costh = -1.0 + (ith - 0.5)*dcosth
      th = dacos(costh)
   DO ips = 1, nbox
      ps = (ips - 0.5)*dps

      i = ips + (ith-1)*nbox + (iph-1)*nbox*nbox

!!! Direction cosine matrix
!!! acs(k,j) = cosine of the angle between the kth crystallographic axes and the
!jth external axes 

      !CALL EULER_TO_DIRCOS(ph,th,ps,acs0sbf(i,:,:))  ! direction cosines g_(ij) of crystal axes

!      if(i==1) then
!        print *,i,iph,ith,ips,ph,th,ps
!        write(*,'(3f12.4)'), acs0(i,1,:)
!        write(*,'(3f12.4)'), acs0(i,2,:)
!        write(*,'(3f12.4)'), acs0(i,3,:)
!stop
!      end if
   ENDDO; ENDDO; ENDDO

   !Compute calc_sj for Olivine
   p21 = dlog(tau(1,2)/tau(1,1))
   p31 = dlog(tau(1,3)/tau(1,1))
   p41 = dlog(tau(1,5)/tau(1,1))

   ! Calculate spin amplitudes for the four slip systems
   !CALL SPIN_AMPLITUDES(p21,p31,p41,stressexp(1),biga_ol)

!  write(6,*) 'A_1-A_4 = ', (biga(j), j = 1,4)

! Calculate partial expansion coefficients \mathcal C_{sj} for three slip
! systems of olivine
   !CALL CALCSJ(biga_ol,calc_ol)

   p21 = 1d0
   p31 = 1d0
   p41 = dlog(1d0/40d0)

! Calculate spin amplitudes for the four slip systems
   !CALL SPIN_AMPLITUDES(p21,p31,p41,stressexp(1),biga_opx)

!  write(6,*) 'A_1-A_4 = ', (biga(j), j = 1,4)

! Calculate partial expansion coefficients \mathcal C_{sj} for three slip
! systems of olivine
   !CALL CALCSJ(biga_opx,calc_opx)

   END IF

!!! Cartesian coordinates of a spherical object    

   degstp=15 !!! step in degrees
   nxy=360/degstp
   nz=90/degstp + 1
    
   ALLOCATE(phix0(nxy),phiy0(nxy),phiz0(nz),phiz1(nz))
   ALLOCATE(phix(nz,nxy),phiy(nz,nxy),phiz(nz,nxy))

   DO i = 1 , nz 
      phi = degstp*(i-1)*acos(-1.0)/180
      phiz0(i) = cos(phi); 
      phiz1(i) = sin(phi); 
   END DO 

   DO i = 1 , nxy
      phi = degstp*(i-1)*acos(-1.0)/180
      phix0(i) = cos(phi); 
      phiy0(i) = sin(phi); 
   END DO 
  
   DO i = 1 , nz
      DO j = 1 , nxy
      phix(i,j) = phix0(j) * phiz1(i)
      phiy(i,j) = phiy0(j) * phiz1(i)
      phiz(i,j) = phiz0(i)           
      END DO
   END DO

   DEALLOCATE(phix0,phiy0,phiz0,phiz1)

   !!! Load mantle density database made with PERPLE_X
   rho(1) = 3353d0
   IF(ptmod > 0) THEN

      !P-T domain grid
      tknum = 85 ; tkmin = 300; tkstp = 50
      pbnum = 1401 ; pbmin = 0d0; pbstp = 1d3
      ptnum = tknum*pbnum

      ALLOCATE(Rhodum(ptnum),Vpdum(ptnum),Vsdum(ptnum))

      CALL H5open_f (error)

      IF(eosmod == 1) CALL H5Fopen_f("../DATABASES/MMA-EoS/dunite.h5", H5F_ACC_RDONLY_F, file_id, error)
      IF(eosmod == 2) CALL H5Fopen_f("../DATABASES/MMA-EoS/hartzburgite.h5", H5F_ACC_RDONLY_F, file_id, error)
      IF(eosmod == 3) CALL H5Fopen_f("../DATABASES/MMA-EoS/pyrolite.h5", H5F_ACC_RDONLY_F, file_id, error)
      IF(eosmod == 4) CALL H5Fopen_f("../DATABASES/MMA-EoS/morb.h5", H5F_ACC_RDONLY_F, file_id, error)
      IF(eosmod == 5) CALL H5Fopen_f("../DATABASES/MMA-EoS/pyroxenite.h5", H5F_ACC_RDONLY_F, file_id, error)

      CALL H5Gopen_f(file_id, "/rock", group_id, error)

      CALL loadsave_double(0,1,group_id,ptnum,H5T_NATIVE_DOUBLE,Rhodum,'Rho',0)
      CALL loadsave_double(0,1,group_id,ptnum,H5T_NATIVE_DOUBLE,Vpdum,'Vp',0)
      CALL loadsave_double(0,1,group_id,ptnum,H5T_NATIVE_DOUBLE,Vsdum,'Vs',0)

      CALL H5Gclose_f(group_id, error)
      CALL H5Fclose_f(file_id, error)
      CALL H5close_f(error)

      ALLOCATE(td_rho(tknum,pbnum),td_vp(tknum,pbnum),td_vs(tknum,pbnum))

      DO j = 1 , tknum ; DO i = 1 , pbnum
 
         gi = j + (i-1)*tknum
         td_rho(j,i) = Rhodum(gi) !kg/m^3
         td_vp(j,i)  = Vpdum(gi)  !km/s
         td_vs(j,i)  = Vsdum(gi)  !km/s

      END DO ; END DO

      DEALLOCATE(Rhodum,Vpdum,Vsdum)

   END IF

   RETURN

   END SUBROUTINE init0
  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!  Rotate tensor parallel to shortes axis of FSE  !!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


   SUBROUTINE tensorrot_aggr(m,acs1)

   USE comvar

   IMPLICIT NONE

   INTEGER :: i,j,k,ll,p,q,r,ss,m
   DOUBLE PRECISION, DIMENSION (3,3,3,3) :: C0,Cav
   DOUBLE PRECISION, DIMENSION(3,3) :: acs1

   !Angles are commonly defined according to the right-hand rule. Namely, they
   !have positive values when they represent a rotation that appears clockwise
   !when looking in the positive direction of the rotating axis, and negative
   !values when
   !the rotation appears counter-clockwise. 

   Cav = 0d0 ; C0 = 0d0

   !Convert to Cijkl 4th order tensor
   DO i = 1 , 3 ; DO j = 1 , 3 ; DO k = 1 , 3 ; DO ll = 1 , 3
      C0(i,j,k,ll) = Sav(ijkl(i,j),ijkl(k,ll),m)
   END DO ; END DO ; END DO ; END DO

   !Rotate 4th order tensor
   DO i = 1 , 3 ; DO j = 1 , 3 ; DO k = 1 , 3 ; DO ll = 1 , 3
      DO p = 1 , 3 ; DO q = 1 , 3 ; DO r = 1 , 3 ; DO ss = 1 , 3
          Cav(i,j,k,ll) = Cav(i,j,k,ll) + acs1(i,p)*acs1(j,q)*acs1(k,r)*acs1(ll,ss)*C0(p,q,r,ss)
      END DO ; END DO ; END DO ; END DO
   END DO ; END DO ; END DO ; END DO

   !Convert to Voigt notation
   DO i = 1 , 6 ; DO j = 1 , 6
      Sav(i,j,m) = Cav(l1(i),l2(i),l1(j),l2(j))
   END DO ; END DO

   RETURN

   END SUBROUTINE tensorrot_aggr

