var/da/da_statistics/da_correlation_coeff2d.inc

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