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    if (trace_use) call da_trace_entry("da_correlation_coeff1d")
26    
27    ! [1.0] Set up scalars:
28    
29    jx = size(field1(:))
30    
31    j_inv = 1.0 / real(jx)
32    
33    ! [2.0] Calculate mean and remove from field:
34    
35    field_mean = sum(field1(:)) * j_inv
36    allocate(data1(1:jx))
37    data1(:) = field1(:) - field_mean
38    
39    field_mean = sum(field2(:)) * j_inv
40    allocate(data2(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 accuracy:
56    
57    coeff3 = sum((data2(:) - data1(:))**2)
58    if (coeff2 /= 0.0) then
59       rel_acc = 1.0 - min(coeff3/coeff2,1.0)
60    else
61       rel_acc = 0.0
62    end if
63    
64    deallocate(data1)
65    deallocate(data2)
66 
67    if (trace_use) call da_trace_exit("da_correlation_coeff1d")
68    
69 end subroutine da_correlation_coeff1d
70 
71