subroutine sim2sd1ss(istat,iterm,verbos,cfld,phifld,psifld,gamfld, & efld,vfld,c,phi,psi,gam,e,v,R,lxfld,thx,thy, & nrfx,nrfy,dx,dy,dmpfld,nfld,jfld,ifld,job,sub1, & sub2,varfnc,kseed,debug,ltanfi,storan) real cfld(nrfx,*), phifld(nrfx,*), psifld(nrfx,*) real gamfld(nrfx,*), efld(nrfx,*), vfld(nrfx,*) real c(*), phi(*), psi(*), gam(*), e(*), v(*), R(6,*), thx, thy,storan(*) integer nrfx, nrfy, nfld, ifld character*(*) job, sub1, sub2, varfnc logical verbos, dmpfld, debug, liid, shofld, lxfld, ltanfi real*8 dvar, dpb, dthx, dthy, dthz, ddx, ddy real*8 dlavx2, dlsep2, dlspx2, dlafr2, dlsfr2 real*8 dmin1, dble save XL, YL, ienter, liid external dlavx2, dlsep2, dlspx2, dlafr2, dlsfr2 common/dparam/ dvar, dpb, dthx, dthy, dthz data zero/0.0/, half/0.5/, one/1.0/, onept5/1.5/, twopt5/2.5/ data twopi/6.2831853071795864769/ data ienter/0/ 1 format(a,a) !-------------------------------------- initialize --------------------- ! compute required field size (once) ienter = ienter + 1 if( ienter .eq. 1 ) then liid = ((thx .eq. zero) .and. (thy .eq. zero)) XL = float(nrfx)*dx YL = float(nrfy)*dy if( debug ) write(istat,1)'SIM2SD: setting field parameters...' if( .not. liid ) then ! set variance fnc parameters dvar = 1.d0 dthx = thx dthy = thy if( varfnc .eq. 'dlafr2' .or. varfnc .eq. 'dlsfr2' ) then dpb = dx if( dy .lt. dx ) dpb = dy endif ! initialize LAS2G if( varfnc .eq. 'dlavx2' ) then call las2g(cfld,nrfx,nrfy,XL,YL,dlavx2,kseed,-1,istat) elseif( varfnc .eq. 'dlsep2' ) then call las2g(cfld,nrfx,nrfy,XL,YL,dlsep2,kseed,-1,istat) elseif( varfnc .eq. 'dlspx2' ) then call las2g(cfld,nrfx,nrfy,XL,YL,dlspx2,kseed,-1,istat) elseif( varfnc .eq. 'dlafr2' ) then call las2g(cfld,nrfx,nrfy,XL,YL,dlafr2,kseed,-1,istat) elseif( varfnc .eq. 'dlsfr2' ) then call las2g(cfld,nrfx,nrfy,XL,YL,dlsfr2,kseed,-1,istat) else if( verbos ) then write(iterm,1) & '*** Error: unknown variance function name: ', varfnc endif write(istat,1) & '*** Error: unknown variance function name: ', varfnc stop endif endif endif ! are we going to plot anything? shofld = dmpfld .and. (nfld .eq. ienter) ! now produce 5 standard normal fields ! for cohesion... if( c(3) .lt. half ) then ! cohesion is deterministic do 10 j = 1, nrfy do 10 i = 1, nrfx cfld(i,j) = zero ! we'll add mean back on later 10 continue elseif( liid ) then ! white noise field do 20 j = 1, nrfy do 20 i = 1, nrfx cfld(i,j) = gausv(one) 20 continue else ! generate standard normal RF if( varfnc .eq. 'dlavx2' ) then call las2gran(cfld,nrfx,nrfy,XL,YL,dlavx2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlsep2' ) then call las2gran(cfld,nrfx,nrfy,XL,YL,dlsep2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlspx2' ) then call las2gran(cfld,nrfx,nrfy,XL,YL,dlspx2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlafr2' ) then call las2gran(cfld,nrfx,nrfy,XL,YL,dlafr2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlsfr2' ) then call las2gran(cfld,nrfx,nrfy,XL,YL,dlsfr2,kseed,0,istat,storan) endif endif ! for friction angle... if( phi(3) .lt. half ) then ! friction is deterministic do 30 j = 1, nrfy do 30 i = 1, nrfx phifld(i,j) = zero ! we'll add mean back on later 30 continue elseif( liid ) then ! white noise field do 40 j = 1, nrfy do 40 i = 1, nrfx phifld(i,j) = gausv(one) 40 continue else ! generate standard normal RF if( varfnc .eq. 'dlavx2' ) then call las2gran(phifld,nrfx,nrfy,XL,YL,dlavx2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlsep2' ) then call las2gran(phifld,nrfx,nrfy,XL,YL,dlsep2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlspx2' ) then call las2gran(phifld,nrfx,nrfy,XL,YL,dlspx2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlafr2' ) then call las2gran(phifld,nrfx,nrfy,XL,YL,dlafr2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlsfr2' ) then call las2gran(phifld,nrfx,nrfy,XL,YL,dlsfr2,kseed,0,istat,storan) endif endif ! for dilation angle... if( psi(3) .lt. half ) then ! dilation is deterministic do 50 j = 1, nrfy do 50 i = 1, nrfx psifld(i,j) = zero ! we'll add mean back on later 50 continue elseif( liid ) then ! white noise field do 60 j = 1, nrfy do 60 i = 1, nrfx psifld(i,j) = gausv(one) 60 continue else ! generate standard normal RF if( varfnc .eq. 'dlavx2' ) then call las2gran(psifld,nrfx,nrfy,XL,YL,dlavx2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlsep2' ) then call las2gran(psifld,nrfx,nrfy,XL,YL,dlsep2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlspx2' ) then call las2gran(psifld,nrfx,nrfy,XL,YL,dlspx2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlafr2' ) then call las2gran(psifld,nrfx,nrfy,XL,YL,dlafr2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlsfr2' ) then call las2gran(psifld,nrfx,nrfy,XL,YL,dlsfr2,kseed,0,istat,storan) endif endif ! for unit weight... if( gam(3) .lt. half ) then ! unit weight is deterministic do 62 j = 1, nrfy do 62 i = 1, nrfx gamfld(i,j) = zero ! we'll add mean back on later 62 continue elseif( liid ) then ! white noise field do 64 j = 1, nrfy do 64 i = 1, nrfx gamfld(i,j) = gausv(one) 64 continue else ! generate standard normal RF if( varfnc .eq. 'dlavx2' ) then call las2gran(gamfld,nrfx,nrfy,XL,YL,dlavx2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlsep2' ) then call las2gran(gamfld,nrfx,nrfy,XL,YL,dlsep2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlspx2' ) then call las2gran(gamfld,nrfx,nrfy,XL,YL,dlspx2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlafr2' ) then call las2gran(gamfld,nrfx,nrfy,XL,YL,dlafr2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlsfr2' ) then call las2gran(gamfld,nrfx,nrfy,XL,YL,dlsfr2,kseed,0,istat,storan) endif endif ! for elastic modulus... if( e(3) .lt. half ) then ! modulus is deterministic do 70 j = 1, nrfy do 70 i = 1, nrfx efld(i,j) = zero ! we'll add mean back on later 70 continue elseif( liid ) then ! white noise field do 80 j = 1, nrfy do 80 i = 1, nrfx efld(i,j) = gausv(one) 80 continue else ! generate standard normal RF if( varfnc .eq. 'dlavx2' ) then call las2gran(efld,nrfx,nrfy,XL,YL,dlavx2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlsep2' ) then call las2gran(efld,nrfx,nrfy,XL,YL,dlsep2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlspx2' ) then call las2gran(efld,nrfx,nrfy,XL,YL,dlspx2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlafr2' ) then call las2gran(efld,nrfx,nrfy,XL,YL,dlafr2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlsfr2' ) then call las2gran(efld,nrfx,nrfy,XL,YL,dlsfr2,kseed,0,istat,storan) endif endif ! for Poisson's ratio... if( v(3) .lt. half ) then ! ratio is deterministic do 90 j = 1, nrfy do 90 i = 1, nrfx vfld(i,j) = zero ! we'll add mean back on later 90 continue elseif( liid ) then ! white noise field do 100 j = 1, nrfy do 100 i = 1, nrfx vfld(i,j) = gausv(one) 100 continue else ! generate standard normal RF if( varfnc .eq. 'dlavx2' ) then call las2gran(vfld,nrfx,nrfy,XL,YL,dlavx2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlsep2' ) then call las2gran(vfld,nrfx,nrfy,XL,YL,dlsep2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlspx2' ) then call las2gran(vfld,nrfx,nrfy,XL,YL,dlspx2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlafr2' ) then call las2gran(vfld,nrfx,nrfy,XL,YL,dlafr2,kseed,0,istat,storan) elseif( varfnc .eq. 'dlsfr2' ) then call las2gran(vfld,nrfx,nrfy,XL,YL,dlsfr2,kseed,0,istat,storan) endif endif ! combine for correlated fields... if( lxfld ) then do 110 j = 1, nrfy do 110 i = 1, nrfx vfld(i,j) = R(1,6)*cfld(i,j) + R(2,6)*phifld(i,j) & + R(3,6)*psifld(i,j) + R(4,6)*gamfld(i,j) & + R(5,6)*efld(i,j) + R(6,6)*vfld(i,j) efld(i,j) = R(1,5)*cfld(i,j) + R(2,5)*phifld(i,j) & + R(3,5)*psifld(i,j) + R(4,5)*gamfld(i,j) & + R(5,5)*efld(i,j) gamfld(i,j) = R(1,4)*cfld(i,j) + R(2,4)*phifld(i,j) & + R(3,4)*psifld(i,j) + R(4,4)*gamfld(i,j) psifld(i,j) = R(1,3)*cfld(i,j) + R(2,3)*phifld(i,j) & + R(3,3)*psifld(i,j) phifld(i,j) = R(1,2)*cfld(i,j) + R(2,2)*phifld(i,j) 110 continue endif ! plot the random field? if( shofld ) then if( jfld .eq. 1 ) then call pltfld( job, sub1, sub2, cfld, nrfx,nrfx,nrfy,XL,YL, & '(Underlying) Cohesion Field', ifld ) elseif( jfld .eq. 2 ) then if( ltanfi ) then call pltfld( job,sub1,sub2,phifld,nrfx,nrfx,nrfy,XL,YL, & '(Underlying) Friction Angle Field', ifld ) else call pltfld( job,sub1,sub2,phifld,nrfx,nrfx,nrfy,XL,YL, & '(Underlying) tan(Friction Angle) Field',ifld) endif elseif( jfld .eq. 3 ) then call pltfld( job, sub1, sub2, psifld,nrfx,nrfx,nrfy,XL,YL, & '(Underlying) Dilation Angle Field', ifld ) elseif( jfld .eq. 4 ) then call pltfld( job, sub1, sub2, gamfld,nrfx,nrfx,nrfy,XL,YL, & '(Underlying) Unit Weight Field', ifld ) elseif( jfld .eq. 5 ) then call pltfld( job, sub1, sub2, efld,nrfx,nrfx,nrfy,XL,YL, & '(Underlying) Elastic Modulus Field', ifld ) elseif( jfld .eq. 6 ) then call pltfld( job, sub1, sub2, vfld,nrfx,nrfx,nrfy,XL,YL, & '(Underlying) Poisson''s Ratio Field', ifld ) endif endif ! convert to final fields ! for cohesion... if( c(3) .lt. half ) then ! deterministic do 120 j = 1, nrfy do 120 i = 1, nrfx cfld(i,j) = c(1) 120 continue elseif( c(3) .lt. onept5 ) then ! normal do 130 j = 1, nrfy do 130 i = 1, nrfx cfld(i,j) = c(1) + c(2)*cfld(i,j) 130 continue elseif( c(3) .lt. twopt5 ) then ! lognormal do 140 j = 1, nrfy do 140 i = 1, nrfx cfld(i,j) = exp( c(4) + c(5)*cfld(i,j) ) 140 continue else ! bounded do 150 j = 1, nrfy do 150 i = 1, nrfx cfld(i,j) = c(4) + half*(c(5)-c(4)) & *( one + tanh((c(6)+c(7)*cfld(i,j))/twopi) ) 150 continue endif ! for friction angle... if( phi(3) .lt. half ) then ! deterministic do 160 j = 1, nrfy do 160 i = 1, nrfx phifld(i,j) = phi(1) 160 continue elseif( phi(3) .lt. onept5 ) then ! normal do 170 j = 1, nrfy do 170 i = 1, nrfx phifld(i,j) = phi(1) + phi(2)*phifld(i,j) 170 continue elseif( phi(3) .lt. twopt5 ) then ! lognormal do 180 j = 1, nrfy do 180 i = 1, nrfx phifld(i,j) = exp( phi(4) + phi(5)*phifld(i,j) ) 180 continue else ! bounded do 190 j = 1, nrfy do 190 i = 1, nrfx phifld(i,j) = phi(4) + half*(phi(5)-phi(4)) & *( one + tanh((phi(6)+phi(7)*phifld(i,j))/twopi) ) 190 continue endif ! for dilation angle... if( psi(3) .lt. half ) then ! deterministic do 200 j = 1, nrfy do 200 i = 1, nrfx psifld(i,j) = psi(1) 200 continue elseif( psi(3) .lt. onept5 ) then ! normal do 210 j = 1, nrfy do 210 i = 1, nrfx psifld(i,j) = psi(1) + psi(2)*psifld(i,j) 210 continue elseif( psi(3) .lt. twopt5 ) then ! lognormal do 220 j = 1, nrfy do 220 i = 1, nrfx psifld(i,j) = exp( psi(4) + psi(5)*psifld(i,j) ) 220 continue else ! bounded do 230 j = 1, nrfy do 230 i = 1, nrfx psifld(i,j) = psi(4) + half*(psi(5)-psi(4)) & *( one + tanh((psi(6)+psi(7)*psifld(i,j))/twopi) ) 230 continue endif ! for unit weight... if( gam(3) .lt. half ) then ! deterministic do 232 j = 1, nrfy do 232 i = 1, nrfx gamfld(i,j) = gam(1) 232 continue elseif( gam(3) .lt. onept5 ) then ! normal do 234 j = 1, nrfy do 234 i = 1, nrfx gamfld(i,j) = gam(1) + gam(2)*gamfld(i,j) 234 continue elseif( gam(3) .lt. twopt5 ) then ! lognormal do 236 j = 1, nrfy do 236 i = 1, nrfx gamfld(i,j) = exp( gam(4) + gam(5)*gamfld(i,j) ) 236 continue else ! bounded do 238 j = 1, nrfy do 238 i = 1, nrfx gamfld(i,j) = gam(4) + half*(gam(5)-gam(4)) & *( one + tanh((gam(6)+gam(7)*gamfld(i,j))/twopi) ) 238 continue endif ! for elastic modulus... if( e(3) .lt. half ) then ! deterministic do 240 j = 1, nrfy do 240 i = 1, nrfx efld(i,j) = e(1) 240 continue elseif( e(3) .lt. onept5 ) then ! normal do 250 j = 1, nrfy do 250 i = 1, nrfx efld(i,j) = e(1) + e(2)*efld(i,j) 250 continue elseif( e(3) .lt. twopt5 ) then ! lognormal do 260 j = 1, nrfy do 260 i = 1, nrfx efld(i,j) = exp( e(4) + e(5)*efld(i,j) ) 260 continue else ! bounded do 270 j = 1, nrfy do 270 i = 1, nrfx efld(i,j) = e(4) + half*(e(5)-e(4)) & *( one + tanh((e(6)+e(7)*efld(i,j))/twopi) ) 270 continue endif ! for Poisson's ratio... if( v(3) .lt. half ) then ! deterministic do 280 j = 1, nrfy do 280 i = 1, nrfx vfld(i,j) = v(1) 280 continue elseif( v(3) .lt. onept5 ) then ! normal do 290 j = 1, nrfy do 290 i = 1, nrfx vfld(i,j) = v(1) + v(2)*vfld(i,j) 290 continue elseif( v(3) .lt. twopt5 ) then ! lognormal do 300 j = 1, nrfy do 300 i = 1, nrfx vfld(i,j) = exp( v(4) + v(5)*vfld(i,j) ) 300 continue else ! bounded do 310 j = 1, nrfy do 310 i = 1, nrfx vfld(i,j) = v(4) + half*(v(5)-v(4)) & *( one + tanh((v(6)+v(7)*vfld(i,j))/twopi) ) 310 continue endif return end