var/da/da_statistics/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