GFTConvolution2D.f90

Go to the documentation of this file.

!==============================================================================
!==============================================================================
SUBROUTINE gft2dconvolution(isp, ny,mx,niy,mjx, spatDataIn,kernelTransf,spatDataOut)
  USE gft, only: gft_prec, gft_rcr, gft_set_fft, gft_do_fft, gft_end_fft
  IMPLICIT NONE
  INTEGER, INTENT(IN)::  isp      ! > species index
  INTEGER, INTENT(IN)::  ny,mx,niy,mjx ! > actual and maximal dimensions
  COMPLEX(kind=GFT_Prec), DIMENSION(0:niy/2,mjx), INTENT(IN) ::   kernelTransf 
  REAL(kind=gft_prec),    DIMENSION(0:ny-1,mx), INTENT(IN) ::   spatdatain 
  REAL(kind=gft_prec),    DIMENSION(0:ny-1,mx), INTENT(OUT)::   spatdataout 
 ! > REAL(kind=GFT_Prec),    DIMENSION(0:ny-1,mx)             ::   spatDataOutTest
 
  REAL(kind=gft_prec),    DIMENSION(0:niy-1,mjx) ::   spatdatainhelp
  REAL(kind=gft_prec),    DIMENSION(0:niy-1,mjx) ::   spatdataouthelpwithkernel  
  REAL(kind=gft_prec)                          :: sc , ratio  
  COMPLEX(kind=GFT_Prec), DIMENSION(0:niy/2,mjx) :: spatDataTransf,spatDataTransfWithKernel 
  
  TYPE(gft_rcr)                                :: rcr
   
  INTEGER :: ilat, ilon, nstart_y, nend_y, mstart_x, mend_x, dimx,dimy   
 
 ! > Here, an additional Filter to get rid of the noise can be created. Currently not used
  ! > CALL LowPassFilter( 0,  niy/2,  1, mjx ,  1.,  filter)
                 
 ! >settings     
  dimx=mx    ! > For GFT_set; the FFT will be done on an intermediate area (ny x mx) that is larger than  
  dimy=ny    ! > the simulation area (maxlon x maxlat) and embedded in a very large area (niy x mix) for zero padding
  
 ! > Spatial data file is padded, i.e. extended to a "good" size by zeros 
 
  ! > bounds where intermediate area is stored in the large area; in corner! >   
  nstart_y = 0; nend_y = ny-1   
  mstart_x = 1; mend_x = mx
 
         ! >   call PrintMatrixASXY2File(spatDataIn,    0,ny-1, 1,mx, "spatDataIn" ) 
 
 ! > Here the intermediate area is stored in the corner of the large area.
  spatdatainhelp(0:,:) =0
  spatdatainhelp(nstart_y : nend_y, mstart_x : mend_x  )=spatdatain     
    
            ! >   call PrintMatrixASXY2File(spatDataInHelp,0,niy-1,1,mjx,"spatDataInHelp" )  
   
 ! > Preparation of FFTs
    
  CALL gft_set_fft(nx=dimy,  ny=dimx,  fft=rcr) ! > x and y exchanged on purpose, because Nx refers to the first dimension
    
 ! >...  forward FFT; transformation from REAL input **  kernel kernelTransf is already transformed and passed *************************    
  
  ! > This is the call to transform the spatial seed production data
 
    sc = 1.0_gft_prec/real(dimx*dimy , kind=gft_prec)     
  CALL gft_do_fft(fft=rcr, isign= -1, scale=sc, r_in=spatdatainhelp, c_out=spatdatatransf)
 
! >   print *, "nach forward transformation"
    
              ! >   call PrintMatrixASXY2File(ABS(spatDataTransf), 0, niy/2, 1,mjx, "ABSspatDataTransf" )   
    
      
! > Actual convolution by multiplication of the two transforms  
   spatdatatransfwithkernel = spatdatatransf  * kerneltransf  
! > 
  
              ! >   call PrintMatrixASXY2File(ABS(kernelTransf ),    0, niy/2, 1,mjx,  "ABSkernelTransf" )
              ! >   call PrintMatrixASXY2File(ABS(spatDataTransfWithKernel),  0, niy /2, 1,mjx,  "ABSspatDataTransfWithKernel" )
   
 
! >  This is the backtransform of the convolution  
 ! >... backward FFT ; transformation back to space domain
 
    sc = 1.0_gft_prec ! >* REAL(1.+ sc_bwd * (precfac - 1.), KIND(GFT)  
    CALL gft_do_fft(fft=rcr, isign=1, scale=sc, c_in=spatdatatransfwithkernel, r_out=spatdataouthelpwithkernel)
     
          ! >   call PrintMatrixASXY2File(spatDataOutHelpWithKernel,0,niy-1,1,mjx, "spatDataOutHelpWithKernel" )    
 
! > Some cleaning: negative and very low values are interpreted as noise and set to zero   
   do ilat = 0, niy-1 ! > y, lat
      do ilon = 1, mjx  ! >1,  m! > x, lon
        if (spatdataouthelpwithkernel(ilat, ilon) < 0 ) spatdataouthelpwithkernel(ilat, ilon)=0   ! > mean trick
        if (spatdataouthelpwithkernel(ilat, ilon) < 1.0e-10 ) spatdataouthelpwithkernel(ilat, ilon)=0   ! > mean trick
      END DO
   END DO   
 
 ! > Get the target area out of the large area 
  spatdataout=spatdataouthelpwithkernel(nstart_y : nend_y, mstart_x : mend_x  )
 
  ! > To make sure that all seeds are distributed the seedrain is scaled that its sum is the sum of the produced seeds
  ratio=sum(spatdatain)/sum(spatdataout)
 
  spatdataout=spatdataout*ratio  ! > Hm, for some reason the seed rain is scaled, probably for distributing exactly all seeds in the area
  
                     ! >   call PrintMatrixASXY2File(spatDataOut, 0,ny-1, 1,mx, "spatDataOut" )
    ! >... Free rcr objects
  CALL gft_end_fft(fft=rcr)
 
  END SUBROUTINE gft2dconvolution
 
!==============================================================================
!==============================================================================
SUBROUTINE gft2dtransformkernel(ny, mx, niy, mjx, y_dim_kernel, x_dim_kernel, kernelIn, kernelTransf )
  USE gft, only: gft_rcr, gft_prec, gft_set_fft, gft_do_fft, gft_end_fft  
  IMPLICIT NONE
  
  INTEGER, INTENT(IN):: ny, mx, niy, mjx, y_dim_kernel, x_dim_kernel
  
  REAL   (kind=gft_prec), DIMENSION(0:(y_dim_kernel-1),x_dim_kernel ), INTENT(IN):: kernelin  
  
  COMPLEX(kind=GFT_Prec), DIMENSION(0:(niy/2),mjx), INTENT(OUT)                  ::  kernelTransf 
  
  REAL   (kind=gft_prec), DIMENSION(0:(niy-1),mjx)                               ::  kernel,kernel1 
  REAL   (kind=gft_prec), DIMENSION(0:(ny-1),mx)                                 ::  mediumsizekernel
  REAL   (kind=gft_prec), DIMENSION(0:(ny-1),mx)                                 ::  reorderedkernel  
  REAL   (kind=gft_prec)   :: sc 
  INTEGER                  :: edge_y, edge_x,edge_y_sm, edge_x_sm, code,  dimy,dimx
 
  TYPE(gft_rcr)            :: rcr
 
    dimy= ny; dimx=mx;
! > print *, "in GFT2DTransformKernel spec, dimy, dimx, niy, mjx  ", dimy, dimx, niy, mjx 
   ! > store small kernel to medium area dimensions
    edge_y = -1  
    edge_x =  0
    kernel(0:,:)=0; kernel1(0:,:)=0;
! >    print *, "edges",edge_y , edge_x,"lengths",niy, mjx, "y_dim_kernel, x_dim_kernel",y_dim_kernel, x_dim_kernel
   ! > reorder kernel, center to corners  
   ! > store shifted kernel central to medium size array of size ny x mx  
    mediumsizekernel(:,:)=0.
    edge_y_sm = (dimy - y_dim_kernel)/2 ! > edges of intermediate area from where on kernel is stored   
    edge_x_sm = (dimx - x_dim_kernel)/2 ! > edges of intermediate area from where on kernel is stored
    mediumsizekernel( (edge_y_sm +1 ):(edge_y_sm + y_dim_kernel) , (edge_x_sm + 1):(edge_x_sm + x_dim_kernel)) = kernelin
                 ! >   call PrintMatrixASXY2File(mediumsizekernel, 0,dimy-1, 1,dimx,  "mediumsizekernel" )
      ! > reorder the kernel on the mediumsized array, center to corners and vice versa
    CALL reorderkernelforgft(dimy,dimx, mediumsizekernel,reorderedkernel ) ! >  
                 ! >   call PrintMatrixASXY2File(reorderedkernel, 0,dimy-1, 1,dimx,  "reorderedkernel" )
  ! > reordered kernel is put in very large extended area, at corners
    kernel((edge_y+1):(edge_y+dimy),(edge_x+1):(edge_x+dimx))= reorderedkernel 
             ! >   call PrintMatrixASXY2File(kernel, 0,niy-1, 1,mjx, "expanded_kernel" ) 
   ! > Preparation of FFT 
 
    CALL gft_set_fft(nx=dimy, ny=dimx, fft=rcr, code=code)    ! >Although the large area is greater, the dimensions are set to those of the medium size. On purpose x and y exchanged! >
   
 ! >...  forward FFT; transformations from REAL  to complex 
    sc = 1.0_gft_prec/real(dimx*dimy , kind=gft_prec) 
! >    print *, "sc", sc
! >  This is the forward transformation of the kernel   
    CALL gft_do_fft(fft=rcr, isign=-1 , scale=sc, r_in=kernel, c_out=kerneltransf) 
             ! >   call PrintMatrixASXY2File(ABS(kernelTransf), 0, niy/2, 1,mjx, "ABSkernelTransf1" )
    CALL gft_end_fft(fft=rcr)
 END SUBROUTINE gft2dtransformkernel
 
!==============================================================================
!==============================================================================
SUBROUTINE setdimsforfft(ispec)
 USE all_par, ONLY:  spec, maxlat, maxlon 
 IMPLICIT NONE                 
 INTEGER, INTENT(IN):: ispec
 INTEGER      ::  safety_margin,   &
                  maxlatEven, maxlonEven, maxlatExt, maxlonExt, maxlatEvenExt, maxlonEvenExt 
 ! > a buffer of one radius of the kernel is put to the simulation area. Not around but only to two sides
  safety_margin =   spec(ispec)%rad 
 ! > next even number including the safety_margin 
  maxlateven=  ((maxlat + safety_margin - 1) / 2) * 2 + 2     
  maxloneven=  ((maxlon + safety_margin - 1) / 2) * 2 + 2     
 ! > next power of two of this buffered area, for the intermediate area
  maxlonevenext=2**ceiling(log(1.*maxloneven)/log(2.))  
  maxlatevenext=2**ceiling(log(1.*maxlateven)/log(2.))  
 ! > and another power of two for the large extended area  
   maxlatext=2**(1+ceiling(log(maxlatevenext*1.)/log(2.)))  ! > next power of two
   maxlonext=2**(1+ceiling(log(maxlonevenext*1.)/log(2.)))  ! > next power of two 
  ! >maxlatExt =2 * maxlatEvenext   ! > next power of two
  ! >maxlonExt =2 * maxlonEvenext   ! > next power of two
  
  print *, "ispec, maxlat, maxlatEvenext(ny),maxlatExt(niy), ",  ispec, maxlat,maxlatevenext, maxlatext 
  print *, "ispec, maxlon, maxlonEvenext(mx),maxlonExt(mjx)",    ispec, maxlon,maxlonevenext, maxlonext
  ! > the dimensions are species specific and stored in the spec object
  spec(ispec)%maxlatExt     = maxlatext
  spec(ispec)%maxlonExt     = maxlonext
  spec(ispec)%maxlatEvenExt = maxlatevenext
  spec(ispec)%maxlonEvenExt = maxlonevenext
 END SUBROUTINE setdimsforfft
 
!==============================================================================
!==============================================================================
 SUBROUTINE lowpassfilter(dimy_start, dimy_end, dimx_start, dimx_end , filterexpo, filter)
  USE gft, only: gft_prec
  IMPLICIT NONE
  
  INTEGER, INTENT(IN):: dimy_start, dimy_end, dimx_start, dimx_end
  REAL, INTENT(IN)                         ::    filterexpo
  COMPLEX(kind=GFT_Prec),DIMENSION(dimy_start:dimy_end,dimx_start:dimx_end), INTENT(OUT)::  filter
   
  REAL(kind=gft_prec)                          ::  distance  
  INTEGER :: iy, jx
 
  DO iy = dimy_start, dimy_end
      DO jx = dimx_start, dimx_end
          distance= ( iy**2 + (jx-1)**2 )**0.5  + 1.
          filter(iy,jx) =  COMPLEX((1./distance**filterexpo), 0.)  
      END DO
  END DO
  filter(:,:) = filter(:,:) /sum(filter(:,:))
  END SUBROUTINE lowpassfilter  
 
 !==============================================================================