da_correlation_coeff1d.inc

1 subroutine da_correlation_coeff1d(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                    :: jx              ! 2nd dimension of field.
16    real                       :: j_inv           ! 1/(jx)
17    real                       :: coeff0          ! coefficient.
18    real                       :: coeff1          ! coefficient.
19    real                       :: coeff2          ! coefficient.
20    real                       :: coeff3          ! coefficient.
21    real                       :: field_mean      ! mean of field.
22    real, allocatable          :: data1(:)
23    real, allocatable          :: data2(:)
24    
25    ! [1.0] Set up scalars:
26    
27    jx = size(field1(:))
28    
29    j_inv = 1.0 / real(jx)
30    
31    ! [2.0] Calculate mean and remove from field:
32    
33    field_mean = sum(field1(:)) * j_inv
34    allocate(data1(1:jx))
35    data1(:) = field1(:) - field_mean
36    
37    field_mean = sum(field2(:)) * j_inv
38    allocate(data2(1:jx))
39    data2(:) = field2(:) - field_mean
40    
41    ! [3.0] Calculate correlation coefficient:
42    
43    coeff0 = sum(data1(:) * data2(:))
44    coeff1 = sum(data1(:) * data1(:))
45    coeff2 = sum(data2(:) * data2(:))
46    
47    if (coeff1 /= 0.0 .and. coeff2 /= 0.0) then
48       corr_coeff = coeff0 /  sqrt(coeff1 * coeff2)
49    else
50       corr_coeff = 0.0
51    end if
52    
53    ! [4.0] Calculate accuracy:
54    
55    coeff3 = sum((data2(:) - data1(:))**2)
56    if (coeff2 /= 0.0) then
57       rel_acc = 1.0 - min(coeff3/coeff2,1.0)
58    else
59       rel_acc = 0.0
60    end if
61    
62    deallocate(data1)
63    deallocate(data2)
64    
65 end subroutine da_correlation_coeff1d
66 
67