Create a topography, two-layer model I need help with what is the wrong with this code

1 

program modelcreater
 implicit none  
  integer :: nx,nz
  integer :: i,j,k,n_rec
  integer :: isx,isz
  real :: sx,sz,dx,dz,offset0,drec
  real,allocatable,dimension(:) :: recx,recz
  character, allocatable,dimension(:,:) :: icode
  real,parameter :: pi=3.1415926
  real :: kk,h
  integer :: mm,ang

   nx=1201                       !grid number in x direction
   nz=601                         !grid number in z direction
   sx=50.0                        !horizontal position of the source (unit: m)
   offset0=1.0                   !the nearest offset in horizontal direction
   drec=1.0                       !the trace interval in horizontal direction
   n_rec=100                    !the total trace numbers
   dx=1.0                             !the horizontal spatial step length (unit: m)
   dz=1.0                             !the vertical spatial step length (unit: m)
   isx=nint(sx/dx)                !the horizontal grid position of the source
   ang=10.0                          !Dipping angle of free surface (just for dipping-free-surface model)
   h=10.0                              !Real thickness of surface layer
   
  allocate(recx(1:n_rec),recz(1:n_rec))
  allocate(icode(1:nx,1:nz))
  
!***********************************************************************************************
!---call the subrotine to configure the matrix of icode  
!call Right_dipping_twolayer(icode,nx,nz,ang,h,pi,dx,dz)     !for twolayer model with right dipping free surface and layer interface
!call Left_dipping_twolayer(icode,nx,nz,ang,h,pi,dx,dz)      !for twolayer model with left dipping free surface and layer interface
!call Left_dipping_twolayer_fault(icode,nx,nz,ang,h,pi,dx,dz) !for twolayer fault model with left dipping free surface and layer interface
call sintopo_twolayer(icode,nx,nz,ang,h,pi,dx,dz)   ! for twolayer model with sin function case topographical free surface
!***********************************************************************************************


!***********************************************************************************************
call threelayer2(icode,nx,nz,ang,h,pi,dx,dz)
! output model the model
  open(1, file='threelayer2.txt',status='unknown')
   do j=1,nz
      write(1,100)(icode(i,j),i=1,nx)
   enddo

   100 format(<nx>a1)
  close(1)
  
  open(2,file='threelayer2_srf.txt',status='unknown')
    do i=1,nx
	   do j=1,nz
	      write(2,200)(i-1)*dx,-1*(j-1)*dz,icode(i,j)
	   enddo
	enddo
  200 format(2e,1x,1a)
  close(2)
  
end program  modelcreater


!!*****************************************************************************************
!!*************************************topography twolayer model***************************

subroutine sintopo_twolayer(icode,nx,nz,ang,h,pi,dx,dz) 
  implicit none
  integer :: nx,nz
  integer :: i,j
  character :: icode(1:nx,1:nz)
  real ::  pi,kk,h,dx,dz
  integer :: ang,mm
           do i=1,nx
               mm=nint(80*sin(pi*(i+200)/400)+0.5)+100
 		     do j=1,nz
 			    if(j<mm)then
 				  icode(i,j)='1'
 				elseif(j>mm.and.j<mm+nint(h/dz))then  !twolayer halfspace with topographic free surface
        		  icode(i,j)='3'
 				elseif(j>=mm+nint(h/dz))then    !twolayer halfspacer with topographic layer interface
 				  icode(i,j)='4'
 				else
 				  icode(i,j)='2'
 				endif
 			 enddo
           enddo
           return
    end subroutine sintopo_twolayer