funcpdf Subroutine

	subroutine funcpdf(ista,iend,n,m,sd,alpha,dx,dy, &
					   z,xprm,yprm,xrec,yrec,f,pl_opt)

	!$ use omp_lib
	implicit none
	integer  :: n, m
	real(dp), dimension (n) :: xprm, yprm, z
	real(dp), dimension (m) :: xrec, yrec, f
	real(dp) :: dx, dy, sd, alpha, soma
	real(dp) :: dxby2, dyby2, x, y, zi, y1, y2, dg
	!$ real(dp) :: t1, t2
	character(len=5) :: loop
	character(len=5), optional :: pl_opt
	integer :: i, j, ista, iend
	
	f(ista:iend)=0._dp
	dxby2=dx/2.0
	dyby2=dy/2.0 
	loop='outer'
	if (present(pl_opt)) loop=pl_opt
	
	if (loop .ne. 'outer')then
		!$ t1 = omp_get_wtime()  
		do i=ista,iend
			soma=0._dp
			!$OMP PARALLEL PRIVATE(j,y,y1,y2,x,zi,dg)
			!$OMP DO REDUCTION(+:soma) schedule (runtime)
			do j=1,n
				y=yprm(j)-yrec(i)
				y1=dyby2+y
				y2=dyby2-y
				x=xprm(j)-xrec(i)
				zi=z(j)
				dg=gr3dprm(x,dxby2,y1,zi,sd,alpha)
				dg=0.5*(dg+gr3dprm(x,dxby2,y2,zi,sd,alpha))
				soma=soma+dg
			enddo
			!$OMP END DO NOWAIT
			!$OMP END PARALLEL 
			f(i)=soma
		enddo
		!$ t2 = omp_get_wtime()
	else
		!$ t1 = omp_get_wtime()  
		!$OMP PARALLEL PRIVATE(i,j,y,y1,y2,x,zi,dg)
		!$OMP DO
		do i=ista,iend
			do j=1,n
				y=yprm(j)-yrec(i)
				y1=dyby2+y
				y2=dyby2-y
				x=xprm(j)-xrec(i)
				zi=z(j)
				dg=gr3dprm(x,dxby2,y1,zi,sd,alpha)
				dg=0.5*(dg+gr3dprm(x,dxby2,y2,zi,sd,alpha))
				f(i)=f(i)+dg
			enddo
		enddo
		!$OMP END DO NOWAIT
		!$OMP END PARALLEL
		!$ t2 = omp_get_wtime()
	endif
	
	!$ print*, "Forward modeling took", t2-t1, "seconds"

	end subroutine funcpdf