PROGRAM inf_tied1
!
!      INFINITE SLOPE STABILITY ANALYSIS
!-------------------------------------------------------------------------
 USE main
! USE geom
 USE gaf95
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 INTEGER::neq,nsrf,nels,nn,nodof=2,nod=4,i,j,iel,nosurf,limit,ndim=2, &
   nlen,maxfld,nrv,nrvp,coy,ifsa,fsp,rvdf
 CHARACTER(LEN=80)::element='quadrilateral',trash,argv
 REAL(iwp)::gamma,pi,h1,w1,s1,area,w,gam_w,fos_tol,srfm,depth,k_h,r_u,    &
   zero=0.0_iwp,one=1.0_iwp,d4=4.0_iwp,sigy,d180=180.0_iwp,pt999=0.999_iwp&
   ,beta,w_d,small=1.0e-6_iwp,mean,sd,pf,huge=1.0e20_iwp,cytol=1.0e-4_iwp,&
   h2,pore,phi,coh,opt5=1.5_iwp,tpt5=2.5_iwp,hh
!---------------------------- dynamic arrays------------------------------
 REAL(iwp),ALLOCATABLE::coord(:,:),g_coord(:,:),       &
   pbuff(:),mcoy(:),vcoy(:),cyrate(:)
 INTEGER,ALLOCATABLE::num(:),g_num(:,:),fnum(:),rvd(:),mlcoy(:),cyapt(:), &
   nf(:,:)
 REAL, ALLOCATABLE::a(:,:),R(:,:),gfld(:,:),ac(:,:,:),fos(:)
 INTEGER kseed,nsim,jr,nseed,iterm,nrfy
 REAL thy,yl,half,onn,threept5,eh,var,v,vl,ul,bm,bsd,bdist,km,ksd,kdist
 LOGICAL debug,lxfld,lsquare
 CHARACTER varfnc*6
 CHARACTER*1  dig(6)
 DATA dig/'1','2','3','4','5','6'/
 data half/0.5/, onn/1.0/,threept5/3.5/
 iterm = 33
 lsquare=.FALSE.
 debug=.FALSE.
 lxfld=.FALSE.
!-------------------------input and initialisation------------------------
 CALL getname(argv,nlen)
 OPEN(10,FILE=argv(1:nlen)//'.dat')
 OPEN(11,FILE=argv(1:nlen)//'.res')
 OPEN(88,FILE=argv(1:nlen)//'.his')
 OPEN(89,FILE=argv(1:nlen)//'.fos')
!
 READ(10,*)trash
 WRITE(11,'(A)')trash
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)h1
 WRITE(11,'(A)')"Height of soil column"
 WRITE(11,'(f8.2)')h1
! READ(10,*)trash
! READ(10,*)beta
! WRITE(11,'(/a)')"Slope angle"
! WRITE(11,'(f8.2)')beta
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)nels
 WRITE(11,'(A)')"Number of finite elements in column"
 WRITE(11,'(I5)')nels
 WRITE(11,*)
 DO i=1,24
  nrfy=2**i
  IF(nrfy>=nels)EXIT
 ENDDO
 maxfld=3*nrfy/2
 yl=h1/DBLE(nels)*DBLE(nrfy)
 READ(10,*)trash
 READ(10,*)nrv
 WRITE(11,'(A)')"Number of random variables"
 WRITE(11,'(I5)')nrv
 WRITE(11,*)
 pi=d4*ATAN(one)
 h2=h1/nels
 w1=h1/10.0_iwp
 nn=2*(nels+1)
 neq=nels*4
 area=w1*h2
!-------------------------------------------------------------------------
 ALLOCATE(g_coord(ndim,nn),nf(nodof,nn),num(nod),coord(nod,ndim),         &
   g_num(nod,nels),gfld(maxfld,nrv),a(7,nrv),r(nrv,nrv),        &
   vcoy(nels),ac(2,nrfy,nrv),fos(nels))
! prop(1)=TAN(phi), prop(2)=c, prop(3)=gamma
!
 READ(10,*)trash
 DO i=1,nrv
  READ(10,*)a(:,i)
 ENDDO
 WRITE(11,'(A)')"nrv rows of material properties (mean, SD, dist, L, U, m, s)"
 DO i=1,nrv
  WRITE(11,'(7f8.2)')a(:,i)
 ENDDO
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)r
 IF(SUM(ABS(r))>DBLE(nrv))lxfld=.TRUE.
 WRITE(11,'(A)')"Correlation matrix"
 DO i=1,nrv
   WRITE(11,'(6f8.2)')r(i,:)
 END DO
 WRITE(11,*)
 DO i=1,nrv
  CALL propst(a(:,i))
 ENDDO
 READ(10,*)trash
 READ(10,*)thy
 WRITE(11,'(A)')"Spatial correlation length"
 WRITE(11,'(f8.2)')thy
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)bm,bsd,bdist
 WRITE(11,'(A)')"Statistics of TAN(beta)"
 WRITE(11,'(3f8.4)')bm,bsd,bdist
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)km,ksd,kdist
 WRITE(11,'(A)')"Statistics of k_h"
 WRITE(11,'(3f8.2)')km,ksd,kdist
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)kseed,nsim
 WRITE(11,'(A)')"Random seed, number of simulations"
 WRITE(11,'(2i5)')kseed,nsim
 WRITE(11,*)
 kseed = iseed(kseed)
 nseed=kseed
 READ(10,*)trash
 READ(10,'(A)')varfnc
 WRITE(11,'(A)')"Covariance function name"
 WRITE(11,'(A)')varfnc
 WRITE(11,*)
! READ(10,*)trash
! READ(10,*)k_h
! WRITE(11,'(A)')"Pseudo horizontal acceleration"
! WRITE(11,'(f5.2)')k_h
! WRITE(11,*)
 READ(10,*)trash; READ(10,*)trash; READ(10,*)trash; READ(10,*)trash
 READ(10,*)nosurf
 WRITE(11,'(A)')"Pore pressure control parameter"
 WRITE(11,*)nosurf
 !   0 means no water pressures
 !   1 means read r_u
 !   2 means read the vertical depth below ground level of a parallel free surface
 !   4 means use the random values generated
 IF(nosurf==1)THEN
   READ(10,*)r_u
   WRITE(11,'(A)')"r_u parameter"
   WRITE(11,'(f5.2)')r_u
 END IF
 IF(nosurf==2)THEN
   READ(10,*)w_d
   WRITE(11,'(A)')"Free-surface depth"
   WRITE(11,'(F5.2)')w_d
 END IF
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)gam_w
 READ(10,*)trash
 READ(10,*)nrvp
 WRITE(11,'(A)')"Number of random variables to be printed"
 WRITE(11,*)nrvp
 WRITE(11,*)
 ALLOCATE(fnum(nrvp),rvd(nrvp),mcoy(nsim),mlcoy(nsim),       &
   pbuff(nsim),cyapt(nsim),cyrate(nsim))
 READ(10,*)trash
 READ(10,*)fnum
 WRITE(11,'(A)')"nvrp random variable numbers"
 WRITE(11,*)fnum
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)rvd
 WRITE(11,'(A)')"Grayscale of random variable (1=log, 0=actual)"
 WRITE(11,*)rvd
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)fsp
 WRITE(11,'(A)')"Create grayscale of FS (fsp=1, print; fsp=0, no print)"
 WRITE(11,*)fsp
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)rvdf
 WRITE(11,'(A)')"Grayscale of log values of FS (rvdf=1) or actual values (rvdf=0)"
 WRITE(11,*)rvdf
 WRITE(11,*)
 eh=yl/nrfy
!---------- adjust parameters for linear case --------------------
 DO i=1,nrv
  IF(a(3,i)>threept5)THEN
   DO iel=1,nrfy
    ac(1,iel,i)=a(1,i)*iel*eh+a(2,i)
    ac(2,iel,i)=a(4,i)*ac(1,iel,i)
    ac(2,iel,i)=alog(onn+ac(2,iel,i)*ac(2,iel,i)/(ac(1,iel,i)*ac(1,iel,i)))
    ac(1,iel,i)=alog(ac(1,iel,i))-half*ac(2,iel,i)
    ac(2,iel,i)=sqrt(ac(2,iel,i))
   END DO
  ENDIF
 END DO
 IF(bdist>half.AND.bdist<opt5)THEN
   var=bsd*bsd
 ELSE IF(bdist<tpt5)THEN
   v=bsd/bm
   vl=ALOG(onn+v*v )
   ul=ALOG(bm)-half*vl
 END IF
!---------- loop the element to set up global arrays --------------------
!---------------start of Monte-Carlo simulations-------------------------
 DO  jr = 1, nsim
  IF( jr/=1)nseed = iseed(kseed + jr - 1)
  CALL  sim1dm1(iterm,nrv,yl,nrfy,thy,varfnc,kseed,lxfld,                      &
                a,ac,R,gfld,nrv,maxfld,lsquare,debug)
!-gravity loads generation and element stiffness integration and assembly-

!  IF(jr==1)WRITE(11,'(6F12.4)')gfld(1:nels,4)

  IF(bdist<half)THEN
    beta=bm
  ELSE IF(bdist<opt5)THEN
    beta=bm+gausv(var)
  ELSE IF(bdist<tpt5)THEN
    beta=EXP(ul+gausv(vl))
  END IF
  beta=ATAN(beta)
!  beta=beta*d180/pi
!  WRITE(*,'(2F8.2)')beta,gfld(1,2)
!  beta=beta*pi/d180
  s1=w1*TAN(beta)
  elements_1: DO iel=1,nels
    CALL inf_geometry_tied(iel,nels,h1,w1,s1,coord,num)
    g_num(:,iel)=num
    g_coord(:,num)=TRANSPOSE(coord)
  END DO elements_1
  k_h=km
!  gravlo=zero
  hh=zero
  sigy=zero
  elements_2: DO iel=1,nels
    phi=ATAN(gfld(iel,1))
    coh=gfld(iel,2)
    gamma=gfld(iel,3)
    hh=hh+h2
    IF(nosurf==0)THEN
      pore=zero
    ELSE IF(nosurf==1)THEN
      sigy=sigy+gamma*h2
      pore=sigy*r_u
    ELSE IF(nosurf==2)THEN
      IF(w_d<zero)w_d=zero
      depth=h2*iel-w_d
      pore=zero
      IF(depth>zero)pore=depth*gam_w*COS(beta)**2
    ELSE IF(nosurf>2)THEN
      pore=gfld(iel,4)
    END IF
    fos(iel)=((hh*gamma*(COS(beta)-k_h*SIN(beta))*COS(beta)-pore)*       &
      TAN(phi)+coh)/(hh*gamma*(SIN(beta)+k_h*COS(beta))*COS(beta))
    IF(jr==nsim)WRITE(89,'(I5,F8.2)')iel,fos(iel)
  END DO elements_2
  pbuff(jr)=MINVAL(fos)
  WRITE(88,'(I5,F8.2,i5)')jr,pbuff(jr),MINLOC(fos)
  IF(jr==1)THEN
   WRITE(11,*)
   WRITE(11,*)
   WRITE(11,*)"---------------For the first simulation---------------------------------------"
   WRITE(11,'(/A)')"Simulation    Random number        FS"
   WRITE(11,'(I10,I17,F12.3)')jr,nseed,pbuff(jr)
  ENDIF
  IF(jr==nsim)THEN
   WRITE(11,*)
   WRITE(11,*)
   WRITE(11,*)"---------------For the last simulation---------------------------------------"
   WRITE(11,'(/A)')"Simulation    Random number        FS"
   WRITE(11,'(I10,I17,F12.3)')jr,nseed,pbuff(jr)
  ENDIF
!
  IF(jr==nsim)THEN
   DO i=1,nrvp
    CALL dismsh1(0.05_iwp,g_coord,g_num,gfld(:,fnum(i)),rvd(fnum(i)),argv, &
      nlen,fnum(i),0,i+20)
    CALL dismsh1(0.05_iwp,g_coord,g_num,gfld(:,fnum(i)),rvd(fnum(i)),argv, &
      nlen,fnum(i),1,i+14)
   ENDDO
   IF(fsp==1)THEN
    IF(rvdf==1)THEN
      CALL dismsh1(0.05_iwp,g_coord,g_num,fos,1,argv, &
       nlen,0,1,90)
     ELSE
      CALL dismsh1(0.05_iwp,g_coord,g_num,fos,0,argv, &
       nlen,0,1,90)
     ENDIF
   ENDIF
  ENDIF
 ENDDO
 WRITE(11,*)
 WRITE(11,*)
 WRITE(11,*)"-------------For the entire suite of Monte-Carlo simulations------------------"
!  CALL mesh(g_coord,g_num,argv,nlen,12)
 mean=SUM(pbuff)/nsim
 sd=SQRT(DOT_PRODUCT(pbuff-mean,pbuff-mean)/DBLE((nsim-1)))
 pf=DBLE(COUNT(pbuff<=one))/DBLE(nsim)
  WRITE(11,'(/A)')"Simulation    Random number      FS_max"
  WRITE(11,'(I10,I17,F12.3)')MAXLOC(pbuff,1),nseed-nsim+MAXLOC(pbuff,1),MAXVAL(pbuff)
  WRITE(11,'(A)')"Simulation    Random number      FS_min"
  WRITE(11,'(I10,I17,F12.3)')MINLOC(pbuff,1),nseed-nsim+MINLOC(pbuff,1),MINVAL(pbuff)
  WRITE(11,'(/A)')"        Pf         mean(FS)       SD(FS)"
  WRITE(11,'(F10.3,F17.3,F12.3)')pf,mean,sd
!
STOP
CONTAINS
!
   SUBROUTINE inf_geometry_tied(iel,nels,h1,w1,s1,coord,num)
!
!      this subroutine forms the coordinates and steering vector
!      for 4-node quadrilaterals in a 'infintite slope' geometry
!      (elements and nodes numbered in the y-direction)
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::h1,w1,s1
 INTEGER,INTENT(IN)::iel,nels
 REAL(iwp),INTENT(OUT)::coord(:,:)
 REAL(iwp)::zero=0.0_iwp
 INTEGER,INTENT(OUT)::num(:)
 INTEGER::iq
 iq=iel
 num(1)=iq*2+1
 num(2)=num(1)-2
 num(3)=num(2)+1
 num(4)=num(1)+1
 coord(1,1)=zero
 coord(2,1)=zero
 coord(3,1)=w1
 coord(4,1)=w1
 coord(1,2)=    -iq*h1/nels
 coord(2,2)=-(iq-1)*h1/nels
 coord(3,2)=coord(2,2)-s1
 coord(4,2)=coord(1,2)-s1
RETURN
END SUBROUTINE inf_geometry_tied
!
SUBROUTINE dismsh1(ratmax,g_coord,g_num,prop,rvd,argv,nlen,fnum,rsgl,ips)
!
! This subroutine produces a PostScript output file "*.dis" displaying
! the deformed finite element mesh with random property field shown as grey-scale
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::g_coord(:,:),ratmax
 REAL prop(:)
 INTEGER,INTENT(IN)::g_num(:,:),ips,nlen,fnum,rvd,rsgl
 CHARACTER(*),INTENT(IN)::argv
 REAL(iwp)::width,height,scale=72,sxy,xo,yo,x,y,dismag,vmax
 REAL(iwp)::xmin,xmax,ymin,ymax,dmax,zero=0.0_iwp,pt5=0.5_iwp,            &
   opt5=1.5_iwp,fpt5=5.5_iwp,d8=8.0_iwp,ept5=8.5_iwp,d11=11.0_iwp,pmin,   &
   pmax,ap,scgr,one=1.0_iwp,gr,small=1.0e-6_iwp,huge=1.0e20_iwp
 INTEGER::i,ii,j,jj,nn,nel,nod
 character*1  dig(6)
 data dig/'1','2','3','4','5','6'/
 IF(ips>20.AND.ips<90)THEN
   OPEN(ips,FILE=argv(1:nlen)//'_'//dig(fnum)//'.mgr')
  ELSEIF(ips>=90)THEN
   OPEN(ips,FILE=argv(1:nlen)//'_fs.mgr')
  ELSE
   OPEN(ips,FILE=argv(1:nlen)//'_'//dig(fnum)//'.dis')
 ENDIF
!
 nn=UBOUND(g_coord,2)
 nel=UBOUND(g_num,2)
 nod=UBOUND(g_num,1)
 xmin=g_coord(1,1)
 xmax=g_coord(1,1)
 ymin=g_coord(2,1)
 ymax=g_coord(2,1)


!                   get range in greyscales
    IF(rvd==1)THEN
      pmin = ALOG(prop(1))
      pmax = ALOG(prop(1))
      DO i = 1, nel
         ap = ALOG(prop(i))
         IF( ap .lt. pmin ) pmin = ap
         IF( ap .gt. pmax ) pmax = ap
      ENDDO
     ELSE
       pmin = prop(1)
       pmax = prop(1)
       DO i = 1, nel
          if( prop(i) .lt. pmin ) pmin = prop(i)
          if( prop(i) .gt. pmax ) pmax = prop(i)
       ENDDO
      ENDIF
      scgr = one/(pmax-pmin)
!
 DO i=2,nn
   IF(g_coord(1,i)<xmin)xmin=g_coord(1,i)
   IF(g_coord(1,i)>xmax)xmax=g_coord(1,i)
   IF(g_coord(2,i)<ymin)ymin=g_coord(2,i)
   IF(g_coord(2,i)>ymax)ymax=g_coord(2,i)
 END DO
 width=xmax-xmin
 height=ymax-ymin
 dmax=ratmax*width
 IF(height>width)dmax=ratmax*height
!
! vmax=zero
! DO i=1,nn
!   DO j=1,2
!     IF(ABS(loads(nf(j,i)))>vmax)vmax=ABS(loads(nf(j,i)))
!   END DO
! END DO
!
! IF(vmax<=small)vmax=huge
!
! dismag=dmax/vmax
!
! xmin=g_coord(1,1)
! xmax=g_coord(1,1)
! ymin=g_coord(2,1)
! ymax=g_coord(2,1)
! DO i=1,nn
!   IF(g_coord(1,i)+dismag*loads(nf(1,i))<xmin)                            &
!     xmin=g_coord(1,i)+dismag*loads(nf(1,i))
!   IF(g_coord(1,i)+dismag*loads(nf(1,i))>xmax)                            &
!     xmax=g_coord(1,i)+dismag*loads(nf(1,i))
!   IF(g_coord(2,i)+dismag*loads(nf(2,i))<ymin)                            &
!     ymin=g_coord(2,i)+dismag*loads(nf(2,i))
!   IF(g_coord(2,i)+dismag*loads(nf(2,i))>ymax)                            &
!     ymax=g_coord(2,i)+dismag*loads(nf(2,i))
!
!   IF(g_coord(1,i)<xmin)xmin=g_coord(1,i)
!   IF(g_coord(1,i)>xmax)xmax=g_coord(1,i)
!   IF(g_coord(2,i)<ymin)ymin=g_coord(2,i)
!   IF(g_coord(2,i)>ymax)ymax=g_coord(2,i)
! END DO
!
 width =xmax-xmin
 height=ymax-ymin
!
!                       allow 1.5" margin minimum on each side of figure
!
!                       portrait mode
!
 IF(height.GE.d11/ept5*width)THEN
!
!                       height governs the scale
!
   sxy=scale*d8/height
   xo=scale*pt5*(ept5-d8*width/height)
   yo=scale*opt5
 ELSE
!
!                       width governs the scale
!
   sxy=scale*fpt5/width
   xo=scale*opt5
   yo=scale*pt5*(d11-fpt5*height/width)
 END IF
!
!                       start PostScript output
!
 WRITE(ips,'(A)')'%!PS-Adobe-1.0'
 WRITE(ips,'(A)')'%%DocumentFonts: none'
 WRITE(ips,'(A)')'%%Pages: 1'
 WRITE(ips,'(A)')'%%EndComments'
 WRITE(ips,'(A)')'/m {moveto} def'
 WRITE(ips,'(A)')'/l {lineto} def'
 WRITE(ips,'(A)')'/s {stroke} def'
 WRITE(ips,'(A)')'/c {closepath} def'
 WRITE(ips,'(A)')'%%EndProlog'
 WRITE(ips,'(A)')'%%Page: 0 1'
 WRITE(ips,'(A)')'gsave'
!
!                       draw the deformed mesh
!
 WRITE(ips,'(2f9.2,a)') xo, yo, ' translate'
 WRITE(ips,'(f9.2,a)') 0.5, ' setlinewidth'
 IF(nod==9)nod=8
 IF(nod==10)nod=9
 IF(nod==15)nod=12
 DO i=1,nel
 !                   overlay greyscale of random property?
    IF( rvd==1 ) THEN       ! map large values to black
         gr = scgr*(pmax - ALOG(prop(i)))
     ELSE
         gr = scgr*(pmax - prop(i))
    ENDIF
    WRITE(ips,'(f9.3,a)')gr,'   setgray'
   ii=g_num(1,i)
   IF(ii==0)CYCLE
   x=sxy*(g_coord(1,ii)-xmin)
   y=sxy*(g_coord(2,ii)-ymin)

   WRITE(ips,'(2f9.2,a)') x, y,' m'
   DO j=2,nod
     jj=g_num(j,i)
     x=sxy*(g_coord(1,jj)-xmin)
     y=sxy*(g_coord(2,jj)-ymin)
     WRITE(ips,'(2f9.2,a)') x, y,' l'
   END DO

    IF(rsgl==1)THEN
!     IF(gr<=small)THEN
!      WRITE(ips, '(A)')'gsave c 0.0 1.0 0.0 setrgbcolor fill grestore'
     IF(gr>=one-small)THEN
      WRITE(ips, '(A)')'gsave c 1.0 0.0 0.0 setrgbcolor fill grestore'
     ELSE
      WRITE(ips, '(A,F5.3,A)')'gsave c ',gr,' setgray fill grestore'
     ENDIF
    ELSE
     WRITE(ips, '(A,F5.3,A)')'gsave c ',gr,' setgray fill grestore'
    ENDIF
  WRITE(ips,'(A)')'c s'
 END DO
 WRITE(ips,'(A)')'grestore'
 WRITE(ips,'(A)')'showpage'
 CLOSE(ips)
!
RETURN
END SUBROUTINE dismsh1
END PROGRAM inf_tied1
