DispersalKernel.f90

Go to the documentation of this file.

! ==============================================================================
! 
!   Module   :  DispersalKernel.f90
! 
!   Remarks:  contains the SUBROUTINEs + functions
!             CalculateAndStoreDispersalKernelOfSpecies (<GetSpc)
!             DispersalKernel (<CalculateAndStoreDispersalKernelOfSpecies)
!             DoubleKernel (<DispersalKernel )
!             SingleKernel (<DispersalKernel AND/OR <DoubleKernel)
!             (F) CenterToCenterValue (<SingleKernel)
!             (F) KernelValue (<CenterToCenterValue )
! ==============================================================================
! ********************************************
 
! ==============================================================================
! ==============================================================================
SUBROUTINE calculateandstoredispersalkernelofspecies  (isp, doFFT)  
    use filelistmodule, only: file, dispersalkernel_file,dispersalkernelnew_file 
    
    use all_par, only:  spec,  & 
            kerneltype,  readdispersalkernel, cellsidelength   
    
    IMPLICIT NONE 
    
    ! Passed variables>--------------------------------------------------------------------------------
    LOGICAL, INTENT(IN)  :: doFFT     ! dispersal with FFT or direct? 
    INTEGER, INTENT(IN)  :: isp       ! species index    
     
    ! Local variables>---------------------------------------------------------------------------------
    INTEGER ::  ii, jj       
  
    LOGICAL     ::  uniqueDispersalKernel  
    ! Here the SUBROUTINE starts>**********************************************************************
    
    uniquedispersalkernel = .true. 
    !-----  If same parameters already in other species, take that kernel and store it      
    do ii = 1, isp - 1 
        if ((spec(isp)%alfa1 == spec(ii)%alfa1) .AND. &
                (spec(isp)%alfa2 == spec(ii)%alfa2) .AND. &
                (spec(isp)%dispFac == spec(ii)%dispFac)) then       ! kernel already exists
            spec(isp)%rad = spec(ii)%rad
            call allocatespeckernel(spec(isp)%rad, isp)    ! allocate the array for the kernel
!             
            !--- store the kernel values of already calculated kernel to the current kernel;
            spec(isp)%kernel(:, :) = 0
            spec(isp)%kernel = spec(ii)%kernel
            uniquedispersalkernel = .false.
            exit   
       end if !  same parameters already exist
    end do !ii
 
    if (readdispersalkernel) then ! if wished, the dispersal kernel is read in, +- never used
        read (dispersalkernel_file%unit, *) spec(isp)%rad
        call allocatespeckernel(spec(isp)%rad, isp)
        do ii = -spec(isp)%rad, spec(isp)%rad
            read (dispersalkernel_file%unit, *) (spec(isp)%kernel(ii, jj), jj=-spec(isp)%rad, spec(isp)%rad)
        end do
    else
        !-----  If these parameters appear first time, kernel has to be calculated
        if (uniquedispersalkernel) then
            kerneltype = 1
            if (spec(isp)%dispFac  <  1.0)  kerneltype = 2
            call dispersalkernel(isp)
        end if ! if uniqueDispersalKernel
       
        ! Write the kernel to "DispersalKernelNew_File", only for the first species
        if (isp == 1) then
            do ii = -spec(isp)%rad, spec(isp)%rad
                do jj = -spec(isp)%rad, spec(isp)%rad
                    write (dispersalkernelnew_file%unit, '(F14.6,4I8,3F14.6)') &
                            cellsidelength, isp, spec(isp)%rad, ii, jj, cellsidelength*ii, &
                            cellsidelength*jj, spec(isp)%kernel(ii, jj)
                end do
            end do
        end if ! isp ==  1 
    end if !  readDispersalKernel 
    if (dofft) then ! if FFT is switched on
        CALL setdimsforfft(isp)  ! the area dimensions for the FFT convolution are determined and stored
        CALL allocatespectransformedkernel( isp ) ! the area for the transformed kernel is allocated
        !  the Fourier Transform of the kernel is calculated and stored to spec(i)%kernel_FT
        CALL calculatekerneltransformation(isp)    
    end if !doFFT
END SUBROUTINE  calculateandstoredispersalkernelofspecies  
 
 
! ==============================================================================
!  ==============================================================================
SUBROUTINE dispersalkernel(isp)
    ! CALL modules for TYPE definitions and +- fixed variables>---------------------------------------------------------------
    use all_par, only: dispersaldifferent, alpha_all, cellsidelength, epskernel, kerneltype, spec, kernel
    
    IMPLICIT NONE
    
    ! <Passed variables>--------------------------------------------------------------------------------
    INTEGER, INTENT(in) :: isp  ! species index
    
    ! <Local variables>---------------------------------------------------------------------------------
    INTEGER :: radi
    REAL :: alfa
    
    ! <Here the SUBROUTINE starts>**********************************************************************
    if (.NOT. dispersaldifferent) then
        call singlekernel(alpha_all, cellsidelength, epskernel, radi)
    elseif (kerneltype == 1) then
        alfa = spec(isp)%alfa1
        call singlekernel(alfa, cellsidelength, epskernel, radi)
    elseif (kerneltype == 2) then
        call doublekernel(spec(isp)%dispFac, spec(isp)%alfa1, spec(isp)%alfa2, cellsidelength, epskernel, radi)
    end if
    call allocatespeckernel(radi, isp)
    spec(isp)%kernel = kernel
    deallocate (kernel)
    spec(isp)%rad = radi
end SUBROUTINE dispersalkernel
 
! ==============================================================================
! ==============================================================================
SUBROUTINE singlekernel(alpha, cellSideLength, epsKernel, rad_c)
    Use all_par, Only: kernel, kernelfine 
    IMPLICIT NONE
    ! <Passed variables>--------------------------------------------------------------------------------
    
    REAL, INTENT(in)  :: alpha, cellSideLength, epsKernel
    INTEGER, INTENT(out) :: rad_c ! > radius of kernel of actual discretization in grid cell number
    ! <Local variables>---------------------------------------------------------------------------------
    INTEGER :: latwin, latwingr, latwinstart, latwinend, &    ! > index of relative window latitude
            lonwin, lonwingr, lonwinstart, lonwinend, &    ! > index of relative window longitude
            icellSideLengthFak, radMax! >, irad                               ! > for switching between min. and actual discretization
    REAL    :: krnlsm, radFac,rad_m, &  
            cellSideLengthMin, &
            CenterToCenterValue             
    
    ! <Here the SUBROUTINE starts>**********************************************************************
    ! ---- determine radius of actual kernel so that kernel reaches only to the value  epskernel
    !      kernel (x) = (1/alpha)* exp(-x/alpha) ==> epskernel = (1/alpha) * exp(- x_epskernel / alpha) 
    !      ==> epskernel * alpha =   exp(- x_epskernel / alpha) ==> ln(epskernel * alpha)= -  x_epskernel / alpha 
    !      ==> x_epskernel = - alpha * ln(epskernel * alpha)
    radfac = min(20., -log(epskernel))                               ! > at which distance (in m) has the kernel function decreased down to the epskernel value?
    rad_m  = alpha * radfac                                            ! > radius of actual discretization (in meters)
    rad_c  = int((rad_m)/cellsidelength + 1.0)                       ! > radius of actual discretization (in grid cell number)
 
    ! ----- Determine settings of fine discretization kernel
    cellsidelengthmin = 1000.0/32.0                                  ! > side length of totally finest discretization, about 0.0001 of kernel cut off, ...
    cellsidelengthmin = min(cellsidelengthmin, cellsidelength)       ! > ... not larger than actual side length
    icellsidelengthfak = cellsidelength/cellsidelengthmin            ! > ratio actual to finest discretization
    radmax = rad_m / cellsidelengthmin + 1                             ! > radius of finest discretization in number of grid cells (The 20 is sufficient at small gridsizes)
    ! ---- allocate and initialize the fine and actual kernel    
    call allocatekernel(radmax, "kernelFine")
    call allocatekernel(rad_c,    "kernel    ")
    
    kernel = 0.0
    kernelfine = 0.0
    ! ----- fine discretization kernel, calculate one quadrant
    kernellat: do latwin = 0, radmax
        kernellon: do lonwin = 0, radmax
            kernelfine(latwin, lonwin) = centertocentervalue(alpha, cellsidelengthmin, latwin, lonwin) ! > this  calculates the values of the kernel
        end do kernellon
    end do kernellat
    ! ----- fine discretization kernel, copy this one quadrant to  other quadrants
    kernellat2: do latwin = -radmax, radmax
        kernellon2: do lonwin = -radmax, radmax
            if ((latwin < 0) .or. (lonwin < 0)) &
                    kernelfine(latwin, lonwin) = kernelfine(abs(latwin), abs(lonwin))
        end do kernellon2
    end do kernellat2
    ! ----- The kernel on the coarse (i.e. actual) grid is calculated by
    ! ----- summing over the kernel values in all fine-grid cells in each coarse-grid cell
    kernellat1: do latwingr = -rad_c, rad_c
        kernellon1: do lonwingr = -rad_c, rad_c
            latwinstart = max(min(radmax, latwingr*icellsidelengthfak - icellsidelengthfak/2), -radmax)
            latwinend   = max(min(radmax, latwingr*icellsidelengthfak + icellsidelengthfak/2), -radmax)
            lonwinstart = max(min(radmax, lonwingr*icellsidelengthfak - icellsidelengthfak/2), -radmax)
            lonwinend   = max(min(radmax, lonwingr*icellsidelengthfak + icellsidelengthfak/2), -radmax)
            do latwin   = latwinstart, latwinend
                do lonwin = lonwinstart, lonwinend
                    kernel(latwingr, lonwingr) = kernel(latwingr, lonwingr) + kernelfine(latwin, lonwin)/2
                end do ! >----- lonwin
            end do ! >----- latwin
        end do kernellon1
    end do kernellat1
    deallocate (kernelfine)
    ! ----- Normalize kernel, sum = 1
    krnlsm = 0
    do latwingr = -rad_c, rad_c
        do lonwingr = -rad_c, rad_c
            krnlsm = krnlsm + kernel(latwingr, lonwingr)
        end do
    end do
    kernel = kernel/krnlsm
end SUBROUTINE singlekernel
 
! ==============================================================================
!    ==============================================================================
SUBROUTINE doublekernel(dispFac, alfa1, alfa2, cellSideLength, epsKernel, rad_c)
    Use all_par, Only: kernel, kernel1, kernel2
    ! <CALL modules for TYPE definitions and +- fixed variables>---------------------------------------------------------------
    IMPLICIT NONE
    
    REAL, INTENT(in)  :: alfa1, alfa2, dispFac, cellSideLength, epsKernel
    ! REAL, allocatable, INTENT(out) :: kernel
    INTEGER, INTENT(out) :: rad_c
    ! <Local variables>---------------------------------------------------------------------------------
    INTEGER :: rad1, rad2
    ! <Here the SUBROUTINE starts >**********************************************************************
    if ((1.0 - dispfac) .gt. 0.0) then
        ! > calculate and store wider kernel
        call singlekernel(alfa2, cellsidelength, epskernel, rad2)
        call allocatekernel(rad2, "kernel2   ")
        kernel2 = kernel
        ! > calculate and store narrower kernel
        call singlekernel(alfa1, cellsidelength, epskernel, rad1)! > narrower kernel
        call allocatekernel(rad1, "kernel1   ") ! > narrower kernel
        kernel1 = kernel
        ! > Put together the double kernel:  dispfac * narrower kernel   + ( 1-dispfac )* wider kernel
        rad_c = rad2
        call allocatekernel(rad_c, "kernel    ")
        kernel = (1.-dispfac)*kernel2  ! > the wider kernel
        kernel(-rad1:rad1, -rad1:rad1) = kernel(-rad1:rad1, -rad1:rad1) + dispfac*kernel1  ! > ...plus the smaller kernel
        deallocate (kernel1)
        deallocate (kernel2)
    else
        call singlekernel(alfa1, cellsidelength, epskernel, rad_c)
    end if
end SUBROUTINE doublekernel
 
! ==============================================================================
! ==============================================================================
 
REAL function centertocentervalue(alpha, cellSideLength, latwin, lonwin)
IMPLICIT NONE
INTEGER, INTENT(in) :: latwin, lonwin
REAL, INTENT(in) :: alpha, cellsidelength
REAL :: dist, &                 !> distance from source to target
        kernelvalue
dist = sqrt((real(latwin)**2) + (real(lonwin)**2))*cellsidelength
centertocentervalue = kernelvalue(dist, alpha)
end function centertocentervalue
 
! ==============================================================================
! ==============================================================================
REAL function kernelvalue(dist, alpha)
IMPLICIT NONE
REAL, INTENT(in) :: dist, alpha
kernelvalue = exp(-dist/alpha)   !  should be set to zero if smaller than epskernel
end function kernelvalue