module curvature
implicit none

integer,private     :: n1,n2,n3,nrad,windx,windz
real   ,private     :: d2,d3,drad,orad,o2,o3

contains

!======================================================================
subroutine curve_init(n1_in,n2_in,n3_in,d2_in,d3_in,o2_in,o3_in,nrad_in,drad_in,orad_in,windx_in,windz_in)
integer     :: n1_in,n2_in,n3_in,nrad_in,windx_in,windz_in
real        :: d2_in,d3_in,drad_in,orad_in,o2_in,o3_in

n1=n1_in;n2=n2_in;n3=n3_in;nrad=nrad_in
d2=d2_in;d3=d3_in;drad=drad_in;orad=orad_in
o2=o2_in;o3=o3_in
windx=windx_in;windz=windz_in
write(0,*) n1,n2,n3,nrad
write(0,*) o2,o3,orad
write(0,*) d2,d3,drad
write(0,*) windx,windz

end subroutine
!======================================================================
subroutine curve_exe(wgcig,semb)
real						:: num,den,z0,center,h,zd,rad
integer						:: i1,i2,i3,irad,izd,minx,maxx,minz,maxz
real						:: wgcig(:,:,:),semb(:,:,:),temp(n1,n3,nrad)
integer, external					:: omp_get_thread_num,omp_get_num_threads

write(0,*) " CURVATURE"
   semb=0.
!!   !$OMP parallel do private(i1,i3,z0,rad,center,num,den,i2,h,zd,izd) ordered schedule(dynamic,1) num_threads(nodes)
   do i1=1,n1
      do i3=1,n3
         z0=(i3-1)*d3+o3
         do irad=1,nrad
            rad=(irad-1)*drad+orad
            center=z0-rad
            num=0.;den=0.
            do i2=1,n2
               h=(i2-1)*d2+o2
               if (rad*rad-h*h < 0) then
                  cycle
               else
                  if (rad .ne. 0) then
                     zd=z0-abs(rad)*sqrt(rad*rad-h*h)/rad
                  else
                     zd=z0
                  end if
                  izd=zd/d3+1
               end if
               if (izd < 1 .or. izd > n3) then
                  cycle
               else
                  num = num + wgcig(i1,i2,izd)
                  den = den + wgcig(i1,i2,izd)**2
               end if
            end do
            semb(i1,i3,irad) = semb(i1,i3,irad) + num**2*den/(den**2+0.01)
         end do
      end do
   end do

   do irad=1,nrad
      do i3=1,n3
         minz=max(1,i3-windz)
         maxz=min(n3,i3+windz)
         do i1=1,n1
            minx=max(1,i1-windx)
            maxx=min(n1,i1+windx)
            temp(i1,i3,irad) = temp(i1,i3,irad) + sum(semb(minx:maxx,minz:maxz,irad))/((maxx-minx+1)*(maxz-minz+1))
         end do
      end do
   end do
   semb=temp

!   do i1=1,n1
!      do i3=1,n3
!                     
!      end do
!   end do
!
!                  rzd = (zd - o3) / d3
!                  idz = rzd
!                  izd = idz + 1
!                  fx = rzd - idz
!                  gx = 1.0 - fx
!               end if
!               if (izd < 1 .or. izd > n3-1) then
!                  cycle
!               else
!                  num = num + fx*wgcig(i1,i2,izd+1) + gx*wgcig(i1,i2,izd)
!                  den = den + fx*wgcig(i1,i2,izd+1)**2 + gx*wgcig(i1,i2,izd)**2
!               end if

end subroutine
!======================================================================
end module
