outer_product.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
SUBROUTINE outer_product(KENS,KX,KY,PX,PY,PA,PC,OPB_CORRELATIONS,HVAR,HOBS)
!---------------------------------------------------------
!
! Computes the outer product of two vectors X and Y
! to produce a matrix A = XY**T
!
!
!                          Jean-Francois MAHFOUF (11/06)
!--------------------------------------------------------
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: kens, kx, ky
LOGICAL, INTENT(IN) :: opb_correlations
 CHARACTER(LEN=3), DIMENSION(KX), INTENT(IN) :: hvar
 CHARACTER(LEN=3), DIMENSION(KY), INTENT(IN) :: hobs
REAL, DIMENSION(KX,KENS), INTENT(IN) :: px
REAL, DIMENSION(KY,KENS), INTENT(IN) :: py
REAL, DIMENSION(KX,KY),  INTENT(OUT) :: pa
REAL, DIMENSION(KY,KY),  INTENT(OUT) :: pc
!
REAL, DIMENSION(KX,KENS) :: zxpert
REAL, DIMENSION(KY,KENS) :: zypert
REAL, DIMENSION(KENS,KY) :: zyt
REAL, DIMENSION(KX) :: zxm
REAL, DIMENSION(KY) :: zym
!
INTEGER :: i, k, l
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('OUTER_PRODUCT',0,zhook_handle)
!
DO l = 1,kx
  !
  zxm(l) = sum(px(l,:))/REAL(kens)!MEAN
  zxpert(l,:) = px(l,:)-zxm(l)     !ANOMALY
  !
ENDDO
!
DO k = 1,ky
  zym(k) = sum(py(k,:))/REAL(kens)
  zypert(k,:) = py(k,:)-zym(k) 
ENDDO
!
IF (opb_correlations) THEN
  !For 2D EnKF
  pa(:,:)=matmul(zxpert(:,:),transpose(zypert(:,:)))
  pc(:,:)=matmul(zypert(:,:),transpose(zypert(:,:)))
  !
ELSE
  !For 1D EnKF
  DO l = 1,kx
    DO k = 1,ky
      IF (hvar(l).EQ.hobs(k)) THEN
        pa(l,k) = dot_product(zxpert(l,:),zypert(k,:))
      ELSE
        pa(l,k) = 0.0        
      ENDIF
    ENDDO
  ENDDO
ENDIF
!
pa = pa / REAL(kens - 1)
pc = pc / REAL(kens - 1)
!
IF (lhook) CALL dr_hook('OUTER_PRODUCT',1,zhook_handle)
!
END SUBROUTINE outer_product