da_correlation_coeff2d.inc

1 subroutine da_correlation_coeff2d(field1, field2, corr_coeff, &
2                                    rel_acc)
3    
4    !--------------------------------------------------------------------------
5    ! Purpose: Calculate correlation coefficient between two fields.
6    !--------------------------------------------------------------------------
7 
8    implicit none
9    
10    real, intent(in)           :: field1(:,:)     ! input field 1.
11    real, intent(in)           :: field2(:,:)     ! input field 2.
12    real, intent(out)          :: corr_coeff      ! correlation coefficient
13    real, intent(out)          :: rel_acc      ! relative error.
14    
15    integer                    :: iy              ! size of field dim1.
16    integer                    :: jx              ! size of field dim2.
17    real                       :: ij_inv          ! 1/(ij)
18    real                       :: coeff0          ! coefficient.
19    real                       :: coeff1          ! coefficient.
20    real                       :: coeff2          ! coefficient.
21    real                       :: coeff3          ! coefficient.
22    real                       :: field_mean      ! mean of field.
23    real, allocatable          :: data1(:,:)
24    real, allocatable          :: data2(:,:)
25    
26    ! [1.0] Set up scalars:
27    
28    iy = size(field1(:,:),dim=1)
29    jx = size(field1(:,:),dim=2)
30    
31    ij_inv = 1.0 / real(iy*jx)
32    
33    ! [2.0] Calculate mean and remove from field:
34    
35    field_mean = sum(field1(:,:)) * ij_inv
36    allocate(data1(1:iy,1:jx))
37    data1(:,:) = field1(:,:) - field_mean
38    
39    field_mean = sum(field2(:,:)) * ij_inv
40    allocate(data2(1:iy,1:jx))
41    data2(:,:) = field2(:,:) - field_mean
42    
43    ! [3.0] Calculate correlation coefficient:
44    
45    coeff0 = sum(data1(:,:) * data2(:,:))
46    coeff1 = sum(data1(:,:) * data1(:,:))
47    coeff2 = sum(data2(:,:) * data2(:,:))
48    
49    if (coeff1 /= 0.0 .and. coeff2 /= 0.0) then
50       corr_coeff = coeff0 /  sqrt(coeff1 * coeff2)
51    else
52       corr_coeff = 0.0
53    end if
54    
55    ! [4.0] Calculate relative error:
56    
57    coeff3 = sum((data2(:,:) - data1(:,:))**2)
58    if (coeff2 /= 0.0) then
59       rel_acc = min(coeff3/coeff2,1.0)
60    else
61       rel_acc = 0.0
62    end if
63    
64    deallocate(data1)
65    deallocate(data2)
66    
67 end subroutine da_correlation_coeff2d
68 
69