GFT_Interact.f90

Go to the documentation of this file.

!===================================================================== 
!===============================================================
SUBROUTINE seedrainbyfft  (nspc, year)
    IMPLICIT NONE
    
    INTEGER,  INTENT(IN) :: nspc, year    ! number of species, year     
    
    INTEGER  ::isp
    !     
    DO isp=1, nspc
        CALL seedrainbyfftforthisspec( isp,year)
    END DO   
    RETURN    
END Subroutine seedrainbyfft
!===================================================================== 
!===============================================================
SUBROUTINE seedrainbyfftforthisspec(isp, year)
    USE gft, only: gft_prec
    USE all_par, only: stategrid, spec, seedrain,  maxlat, maxlon, &
            stockability, seednumberthresh, epskernel    
 
             
    IMPLICIT NONE
    INTEGER,         INTENT(IN) :: isp, year  ! number of species 
    
    REAL(kind=gft_prec) ::  seedprodarea(0:spec(isp)%maxlatEvenext-1, 1: spec(isp)%maxlonEvenext), &
                            seedrainarea(0:spec(isp)%maxlatEvenext-1, 1: spec(isp)%maxlonEvenext) 
    real*8  ::  help(1:maxlat,1:maxlon),  verynewseedsarray(1:maxlat, 1:maxlon)   
    
    
    COMPLEX(kind=GFT_Prec):: kernel_FT(0:spec(isp)%maxlatExt/2,spec(isp)%maxlonExt)
    
    
    INTEGER :: ilat, ilon, latStart_Area, latEnd_Area,lonStart_Area, lonEnd_Area
    REAL::producedSeedssum,producedSeedsav, presThresh  
    
    seedrainarea(:,:) =0 
    
    producedseedssum=sum(stategrid(1:maxlat,1:maxlon)%sp(isp)%producedSeeds)
    if (producedseedssum==0) THEN ! if no produced seeds, the seed rain is zero
        seedrain(:,:)%sp(isp)%verynewSeeds = 0
    ELSE   ! otherwise the convolution is done      
        latstart_area=   0        
        latend_area=     maxlat-1  
        lonstart_area=   1         
        lonend_area=     maxlon    
        help=   stategrid(1:maxlat,1:maxlon)%sp(isp)%producedSeeds    
               !   call PrintMatrixASXY2File( help,    1,maxlat , 1,maxlon, "seedprod" )         
        
        ! store seed production to the small extended area, just for giving it to the convolution subprogram
        seedprodarea(:,:) = 0. 
        seedprodarea(latstart_area:latend_area,lonstart_area:lonend_area)=stategrid(1:maxlat,1:maxlon)%sp(isp)%producedSeeds
       !   call PrintMatrixASXY2File( seedProdArea,    0,spec(isp)%maxlatEvenext-1 , 1,spec(isp)%maxlonEvenext, "seedprodlarge" )         
    
           
        ! store the kernel, just for giving it to the convolution subprogram
        kernel_ft=spec(isp)%kernel_FT 
        
         !   call PrintMatrixASXY2File(ABS(kernel_FT),0,spec(isp)%maxlatExt/2,1,spec(isp)%maxlonExt,"ABSkernel_FT_SR" )         
        
        ! do the convolution 
        CALL gft2dconvolution(isp, &
                spec(isp)%maxlatEvenext,spec(isp)%maxlonEvenext,&
                spec(isp)%maxlatExt,    spec(isp)%maxlonExt, &
                seedprodarea, kernel_ft ,seedrainarea)
        ! This is a trick to "clean" the strange noise coming from the GFT convolution; helps against the noise, but slows down the migration!  
        presthresh = 500! 200. ! 600 is +- okay, a bit slow 
        producedseedsav = producedseedssum/(maxlat*maxlon)
        do ilat =latstart_area, latend_area
            do ilon= lonstart_area, lonend_area
                if (seedrainarea(ilat, ilon)/producedseedsav <  epskernel * presthresh)  &  !200 threshold found out by trying out
                        seedrainarea(ilat, ilon)=0
            END DO
        END DO  
             !   call PrintMatrixASXY2File(seedRainArea,0,spec(isp)%maxlatEvenext-1,1,spec(isp)%maxlonEvenext, "seedRainArea")       
        
        ! store back the calculated seedrain
        verynewseedsarray = seedrainarea((latstart_area ) :(latend_area ), lonstart_area :lonend_area)
                !   call PrintMatrixASXY2File(verynewseedsArray,    1,maxlat, 1,maxlon, "verynewseedsArray")  
        
        seedrain(:,:)%sp(isp)%verynewSeeds = 0.
        seedrain(:,:)%sp(isp)%verynewSeeds = verynewseedsarray
        
        seedrain(:,:)%sp(isp)%verynewSeeds = seedrain(:,:)%sp(isp)%verynewSeeds * stockability(:,:)    
    ENDIF
    RETURN
END Subroutine seedrainbyfftforthisspec
!===================================================================== 
!===============================================================
SUBROUTINE printmatrixasxy2file(array, dimystart, dimyend, dimxstart, dimxend, varname)
    USE gft, only: gft_prec
    IMPLICIT NONE
    REAL(kind=gft_prec), DIMENSION(dimystart:dimyend,dimxstart:dimxend), INTENT (in):: array
    INTEGER, INTENT(in):: dimystart,dimyend,dimxstart,dimxend
    CHARACTER(*), INTENT(in):: varname !, filename
    
    CHARACTER*50  ::   filenamelong
    
    INTEGER:: i,j
    filenamelong = "AreaOutputs/"//varname//".txt" 
    print *,  varname, " printing" 
    open(17,file=  trim(filenamelong) )
    write(17,*) "i,j,", varname 
    do i= dimystart,dimyend  
        do j= dimxstart,dimxend
            write(17,*)  i,",",j,",",array(i,j) 
        end do             
    end do 
    close(17)
    return
END SUBROUTINE printmatrixasxy2file
 
!==============================================================================
!==============================================================================
Subroutine calculatekerneltransformation(isp)
    ! <CALL modules for TYPEdefinitions>---------------------------------------------------------------    
    USE gft, only: gft_prec
    USE all_par, ONLY: spec                                                   
    
    IMPLICIT NONE
    
    ! <Passed variables>--------------------------------------------------------------------------------
    INTEGER, INTENT(IN) :: isp  
    
    ! <Local variables>---------------------------------------------------------------------------------
    INTEGER ::  kerneldim, y_dim_kernel_orig,x_dim_kernel_orig,y_dim_kernel,x_dim_kernel, rad      
    
    REAL   (kind=gft_prec), allocatable, dimension(:,:)::    kernel_shifted, kernel_orig 
    COMPLEX(kind=GFT_Prec), allocatable, dimension(:,:)::    kernel_FT
    
    
    ! calculate  the dimensions for the kernel arrays
    rad = spec(isp)%rad
    kerneldim= 2 * rad+1
    y_dim_kernel = ((kerneldim-1)/2) * 2 +2  ! next even number
    x_dim_kernel = ((kerneldim-1)/2) * 2 +2  ! next even number
 
    y_dim_kernel_orig=kerneldim
    x_dim_kernel_orig=kerneldim
  ! allocate the kernel arrays   
    allocate(        kernel_shifted(0:y_dim_kernel-1,1:x_dim_kernel))
    allocate(        kernel_orig(kerneldim,kerneldim))
    allocate(        kernel_ft(0 : spec(isp)%maxlatExt/2, 1 : spec(isp)%maxlonExt))
  
    kernel_orig(:,:) =0.; kernel_shifted(:,:)=0.; kernel_ft(0:,:)=0.; spec(isp)%kernel_FT (0: ,:)=0.
    
    kernel_orig=spec(isp)%kernel(-rad:rad, -rad:rad)
 
                 !   call PrintMatrixASXY2File(kernel_orig, 1,kerneldim, 1,kerneldim, "kernel_orig1" )
    
    ! The kernel dimensions need to be even, therefore shifted
    
    CALL shiftkerneltoeven(isp,y_dim_kernel_orig,x_dim_kernel_orig,&
            y_dim_kernel,     x_dim_kernel,&
            kernel_orig,kernel_shifted)       
                 !   call PrintMatrixASXY2File(kernel_shifted, 0,y_dim_kernel-1,  1,x_dim_kernel, "kernel_shifted" )
    
  
   ! this is the transformation
    CALL gft2dtransformkernel(spec(isp)%maxlatEvenExt,spec(isp)%maxlonEvenExt, &
                              spec(isp)%maxlatExt,    spec(isp)%maxlonExt  ,&
                              y_dim_kernel,x_dim_kernel, kernel_shifted, kernel_ft)
  
   ! here the result, i.e. the kernel in frequency domain is stored
    spec(isp)%kernel_FT=kernel_ft
    
   ! free the kernel arrays      
    deallocate(kernel_shifted)
    deallocate(kernel_ft)
    deallocate(kernel_orig)
    return
END Subroutine  calculatekerneltransformation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!