# I need some help to modify this code (simple 2D model) for three layers with fault angle = 10m; nx= 1000; nz=500 **URL:** \ **Category:** Uncategorized\ **Created:** [July 13, 2022, 8:44am UTC](https://fortran-lang.discourse.group/t/i-need-some-help-to-modify-this-code-simple-2d-model-for-three-layers-with-fault-angle-10m-nx-1000-nz-500/3987 "2022-07-13T08:44:11Z")\ **Posts on this page:** 6\ **Page:** 1
**Author:** ![1988](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/1988/32/1103_2.png) [@1988](https://fortran-lang.discourse.group/u/1988)\ **Post date:** [July 13, 2022, 8:44am UTC](https://fortran-lang.discourse.group/t/i-need-some-help-to-modify-this-code-simple-2d-model-for-three-layers-with-fault-angle-10m-nx-1000-nz-500/3987/1 "2022-07-13T08:44:11Z")
```fortran implicit none integer:: i,j,nx,nz character*1,allocatable,dimension(:,:) :: icode real,parameter :: pi=3.1415926 integer :: mm,ang real :: kk,h,dx,dz nx=1000 !grid number in x direction nz=500 !grid number in z direction dx=1.0 !the horizontal spatial step length (unit: m) dz=1.0 !the vertical spatial step length (unit: m) ang=10.0 !Dipping angle of free surface (just for dipping-free-surface model) h=10.0 !Real thickness of surface layer allocate(icode(1:nx,1:nz)) mm=(h/(cos(ang*pi/180))/dx)+1 kk=tan(ang*pi/180) do i=1,(nx-1)/2 do j=1,nz if(jkk*(nx-i)+1).and.(j=kk*(nx-i)+mm)then icode(i,j)='4' else icode(i,j)='2' endif enddo enddo do i=(nx-1)/2+1,nx do j=1,nz if(jkk*(nx-i)+1).and.(j=kk*(nx-i)+(mm-1)/2+1)then icode(i,j)='4' else icode(i,j)='2' endif enddo enddo open(2012,file='model.txt',status='unknown') do j=1,nz write(2012,1000)(icode(i,j),i=1,nx) enddo close(2012) 1000 format(a1) ``` ![Snipaste_2022-07-13_16-35-57](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/3/39005dee8b25deb8fb950121a6b97aaba85fb4d7.png) ---
**Author:** ![1988](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/1988/32/1103_2.png) [@1988](https://fortran-lang.discourse.group/u/1988)\ **Post date:** [July 13, 2022, 8:53am UTC](https://fortran-lang.discourse.group/t/i-need-some-help-to-modify-this-code-simple-2d-model-for-three-layers-with-fault-angle-10m-nx-1000-nz-500/3987/2 "2022-07-13T08:53:11Z")
but when I run the code, the outcome is like this one ![Snipaste_2022-07-13_16-52-36](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/a/a9e0b0af83dd2af35bcb056214104c6ef932a99d.png) ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [July 13, 2022, 5:11pm UTC](https://fortran-lang.discourse.group/t/i-need-some-help-to-modify-this-code-simple-2d-model-for-three-layers-with-fault-angle-10m-nx-1000-nz-500/3987/3 "2022-07-13T17:11:10Z")
Hi Mohammed, Given your previous questions, - [How to make use of cosine Function to make this shape](https://fortran-lang.discourse.group/t/how-to-make-use-of-cosine-function-to-make-this-shape/3764/8) - [How I can draw a surface wave model like the figure below?](https://fortran-lang.discourse.group/t/how-i-can-draw-a-surface-wave-model-like-the-figure-below/3783/15) - [How can I apply this function Y = d+asin(bx-c) in the code shown below? Where a represents the amplitude](https://fortran-lang.discourse.group/t/how-can-i-apply-this-function-y-d-asin-bx-c-in-the-code-shown-below-where-a-represents-the-amplitude/3787/4) - [Create a topography, two-layer model I need help with what is the wrong with this code](https://fortran-lang.discourse.group/t/create-a-topography-two-layer-model-i-need-help-with-what-is-the-wrong-with-this-code/3813/4) - [How to type this part of the code correctly?](https://fortran-lang.discourse.group/t/how-to-type-this-part-of-the-code-correctly/3860/9) - [Create Rectangular Cavity](https://fortran-lang.discourse.group/t/create-rectangular-cavity/3795/11) - [I need to model an incline rectangle shape](https://fortran-lang.discourse.group/t/i-need-to-model-an-incline-rectangle-shape/3984/14) I’d really encourage you to pursue a more modular and flexible program structure, and also to search for better ways (as in geometrical representations) to describe your models. In most of your threads the convoluted if-elseif blocks and conditions are very difficult to reason about. I have the feeling like you don’t really care about this, since these are probably one-off calculations. I tried to demonstrate to you a different approach yesterday with the [point-in-polygon](https://fortran-lang.discourse.group/t/i-need-to-model-an-incline-rectangle-shape/3984/8) test. Today I’m sharing with you an implicit line representation I once found in the [lectures notes on Analysis of Binary Images](https://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/OWENS/LECT2/node3.html#SECTION00031100000000000000) by [Robyn Owens](https://research-repository.uwa.edu.au/en/persons/robyn-owens) from the University of Western Australia. ```fortran program fault implicit none integer, parameter :: nx = 1001, nz = 501 real, parameter :: pi = 4.*atan(1.) real :: depth, theta_deg, theta, x0, z0, rho real :: dx, dz character(1), allocatable :: icode(:,:) ! grid spacing dx = 1. dz = 1. allocate(icode(1:nx,1:nz)) ! fault depth depth = 100. ! fault angle theta_deg = 10. ! degrees theta = theta_deg * (pi/180.) ! radian ! fault cuts through this point x0 = 500. z0 = -50. ! ! the fault is described by the implicit line ! ! x*sin(phi) - z*cos(phi) + rho = 0 ! ! where phi = pi - theta. The meanings of the parameters are as follows ! ! rho ... distance from the origin (0,0) to the closest point ! on the line ! phi ... angle between the line and the x-axis ! theta ... adjacent angle to phi ! ! a full description of the math describing the line is available at ! https://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/OWENS/LECT2/node3.html#SECTION00031100000000000000 ! ! calculate line constant from the known angle, and a point ! lying on the fault rho = z0*cos(pi - theta) - x0*sin(pi - theta) ! fill icode mask call two_layers_with_fault( & ndim = [nx,nz], & spacing = [dx,dz], & icode = icode, & rho = rho, & theta = theta, & depth = depth) call save_txt('model.txt',icode) call save_pgm('model.pgm',icode,ngray=3) contains !> Fill mask array for model with two horizontal layers and a fault !> !> The mask values are !> - 0 or 1 (top layer), depending on the side of the fault line !> - 2 (bottom layer) !> subroutine two_layers_with_fault(ndim,spacing,icode,rho,theta,depth) integer, intent(in) :: ndim(2) ! (nx, nz) real, intent(in) :: spacing(2) ! (dx, dz) character(1), intent(out) :: icode(:,:) real, intent(in) :: rho, theta, depth ! geometry parameters real :: xx, zz integer :: j, i do j = 1, ndim(2) do i = 1, ndim(1) ! we assume implicitly that point (i=1,j=1) is the origin (0,0), ! increasing j corresponds to negative z-values in ! physical coordinate, xx = (i-1) * spacing(1) zz = 0.0 - (j-1) * spacing(2) if (abs(zz) < depth) then ! fault layer icode(i,j) = merge('0','1', & is_below_fault(xx,zz,rho,theta)) else ! bottom layer icode(i,j) = '2' end if end do end do end subroutine !> Check if a point is below the fault line !> !> The line must be prescribed by the implicit form !> !> x * sin(pi - theta) - y * cos(pi - theta) + rho = 0 !> pure function is_below_fault(xx,zz,rho,theta) result(res) real, intent(in) :: xx, zz ! point of interest real, intent(in) :: rho ! line constant real, intent(in) :: theta ! angle in radians real, parameter :: pi = 4.*atan(1.) logical :: res if (xx*sin(pi - theta) - zz*cos(pi - theta) + rho < 0) then res = .true. else res = .false. end if end function !> Save as dumb ASCII txt file (use .txt extension) subroutine save_txt(filename,icode) character(len=*), intent(in) :: filename character(1), intent(in) :: icode(:,:) integer :: i,j,unit open(newunit=unit,file=filename,status='unknown') do j = 1, nz write(unit,'(*(A1))') (icode(i,j), i = 1, nx) end do close(unit) end subroutine !> Save as plain Portable GrayMap (use .pgm extension) subroutine save_pgm(filename,icode,ngray) character(len=*), intent(in) :: filename character(1), intent(in) :: icode(:,:) integer, intent(in) :: ngray ! number of gray values integer :: i, j, unit open(newunit=unit,file=filename,status='unknown') ! a description of the format is available at ! https://en.wikipedia.org/wiki/Netpbm ! header part write(unit,'(A)') 'P2' write(unit,'(I0,1X,I0)') size(icode,1), size(icode,2) write(unit,'(I0)') ngray ! pixel data, gray values do j = 1, size(icode,2) write(unit,'(*(A1,:,1X))') (icode(i,j), i = 1, size(icode,1)) end do close(unit) end subroutine end program ``` Picking the right mathematical representation can make the job much easier for you. Here’s the output I got in the second try, after adjusting a `<` symbol in the `is_below_fault` function. ![model](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/1/174c4ccddc4a0fb253aba11cb50b568b32a566b8.jpeg) ---
**Author:** ![1988](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/1988/32/1103_2.png) [@1988](https://fortran-lang.discourse.group/u/1988)\ **Post date:** [July 14, 2022, 12:28am UTC](https://fortran-lang.discourse.group/t/i-need-some-help-to-modify-this-code-simple-2d-model-for-three-layers-with-fault-angle-10m-nx-1000-nz-500/3987/4 "2022-07-14T00:28:23Z")
Thank you very much for helping me ---
**Author:** ![1988](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/1988/32/1103_2.png) [@1988](https://fortran-lang.discourse.group/u/1988)\ **Post date:** [July 14, 2022, 12:33am UTC](https://fortran-lang.discourse.group/t/i-need-some-help-to-modify-this-code-simple-2d-model-for-three-layers-with-fault-angle-10m-nx-1000-nz-500/3987/5 "2022-07-14T00:33:41Z")
I study engineering geophysics; I use the surface seismic wave method to analyze the structures near the surface, and this is the first time I have dealt with Fortran language. ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [July 19, 2022, 10:32pm UTC](https://fortran-lang.discourse.group/t/i-need-some-help-to-modify-this-code-simple-2d-model-for-three-layers-with-fault-angle-10m-nx-1000-nz-500/3987/6 "2022-07-19T22:32:44Z")
I hope the snippets shown have been useful. Besides from the book [Computational Seismology](http://www.computational-seismology.org/), 2016, by Heiner Igel, I’ve also seen a few seismologically-suited methods in [A Primer on Radial Basis Functions with Applications to the Geosciences](https://doi.org/10.1137/1.9781611974041), 2015, by Fornberg & Flyer. A slightly [younger paper](https://doi.org/10.1016/j.jcp.2017.01.065), by Martin & Fornberg, 2017, discusses how to derive high-accuracy finite difference approximations for material interfaces (discontinuities). For a state of the art seismology code available today, you can have a look at [SeisSol](http://www.seissol.org/). This is a C++ code which has received a lot of attention, and is able to run efficiently on large supercomputer. There’s always something to learn by peaking underneath the surface (pun intended).