!$id: $
! Copyright (c) 2011-2014 Pierre De Mey, pdm789@gmail.com
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU Lesser General Public License
! as published by the Free Software Foundation, either version
! 3 of the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public
! License along with this software.  If not, see <http://www.gnu.org/licenses/>.

! ----------------------------------------------------------------------
  PROGRAM example_arm
! ----------------------------------------------------------------------

! AUTHOR
! Pierre De Mey

! REVISION HISTORY
! 2014-03: P De Mey V1 SANGOMA release

! PURPOSE
! Example on how to CALL sangoma_arm().  An example Ensemble of 3 members
! is provided, along with an array of 2 observations and a diagonal
! observational ECM.  The results show an array mode spectrum with only
! one array mode (eigenvalue=5.5064) above the observational noise floor.

! NOTES
! 1. Input function calls (get*) should be replaced by code to read your 
!    files

  USE sangoma_base, ONLY: REALPREC, INTPREC

  IMPLICIT NONE

! input to program

  INTEGER, PARAMETER :: n = 2 ! state size
  INTEGER, PARAMETER :: m = 3 ! Ensemble size
  INTEGER, PARAMETER :: p = 2 ! number of observations

  CHARACTER(LEN=255) :: fileroot = "example_arm_output"

  REAL(KIND=REALPREC), DIMENSION(n,m) :: ssamples ! Ensemble departures from mean
  REAL(KIND=REALPREC), DIMENSION(p,m) :: dsamples ! same as ssamples in data space
  REAL(KIND=REALPREC), DIMENSION(p,p) :: R ! observation error covariance matrix

! output by program

  INTEGER, PARAMETER :: ndof = MIN(m,p) ! number of d.o.f.s of problem

  INTEGER :: status
  REAL(KIND=REALPREC), DIMENSION(ndof) :: spect ! array mode spectrum
  REAL(KIND=REALPREC), DIMENSION(p,ndof) :: U ! array modes
  REAL(KIND=REALPREC), DIMENSION(n,ndof) :: rho_mu ! modal representers

! local

  INTEGER :: k, l
  REAL(KIND=REALPREC), DIMENSION(p,n) :: H ! observation operator
  CHARACTER(LEN=255) :: filename

! ----------------------------------------------------------------------

  WRITE (*,'(/a/)') "EXAMPLE_ARM *STARTING UP*"

! get state samples ssamples

  WRITE (*,'(a)') "Getting state samples..."
  CALL getStateSamples

! get observation operator H

  WRITE (*,'(a)') "Getting observation operator..."
  CALL getObsOp

! generate data-space samples dsamples = H*ssamples
! (alternately these can be can be directly generated by the model)

  WRITE (*,'(a)') "Generating data-space samples..."
  CALL dgemm( "N","N",p,m,n,1.0,H,p,ssamples,n,0.0,dsamples,p )

! get observation error covariance matrix

  WRITE (*,'(a)') "Getting observation error covariance matrix..."
  CALL getObsECM

! perform stochastic array mode analysis

  WRITE (*,'(a)') "Performing stochastic array mode analysis..."
  CALL sangoma_arm ( n,m,p,ndof,ssamples,dsamples,R,spect,U,rho_mu,status )

  IF ( status /= 0 ) THEN
    WRITE (*,'(/a,i8)') "...sangoma_arm() returned with error code ",status
    WRITE (*,'(a)') "...aborting!"
    STOP
  ELSE
    WRITE (*,'(/a)') "...sangoma_arm() returned with no error"
  END IF

! calculate number of detectable array modes above observational noise floor

  k = 0; DO l = 1, ndof; IF ( spect(l) >= 1 ) k = k + 1; END DO
  WRITE (*,'(a,i8)') "...number of detectable array modes =",k

! save results

  WRITE (*,'(/a)') "Saving results..."
  CALL saveResults

  WRITE (*,'(/a)') "EXAMPLE_ARM *ALL DONE*"

CONTAINS

! NOTE: for the following set of inputs, you should get 1 detectable array
! mode, with corresponding RM spectrum value of 5.5064 .

! ----------------------------------------------------------------------
  SUBROUTINE getStateSamples
! ----------------------------------------------------------------------

  ! state samples ssamples(n,m) should be centered
  ! (relative to Ensemble mean)

  ssamples(1,1) = 0.24
  ssamples(2,1) = 0.51

  ssamples(1,2) = -0.38
  ssamples(2,2) = -0.75

  ssamples(1,3) = 0.14
  ssamples(2,3) = 0.24

  END SUBROUTINE getStateSamples

! ----------------------------------------------------------------------
  SUBROUTINE getObsOp
! ----------------------------------------------------------------------

  ! observation operator H(p,n)

  H(1,1) = 0.0
  H(1,2) = 1.0

  H(2,1) = 1.0
  H(2,2) = 0.0

  END SUBROUTINE getObsOp

! ----------------------------------------------------------------------
  SUBROUTINE getObsECM
! ----------------------------------------------------------------------

  ! observational error covariance matrix (ECM)

  R = 0.0
  DO k = 1, p; R(k,k) = 0.1; END DO

  END SUBROUTINE getObsECM

! ----------------------------------------------------------------------
  SUBROUTINE saveResults
! ----------------------------------------------------------------------

  ! save RM spectrum spect(ndof)

  OPEN (1,file=TRIM(fileroot)//"--RM_spectrum",status="unknown")
  DO l = 1, ndof
    WRITE (1,*) l,spect(l)
  END DO
  CLOSE (1)
  WRITE (*,'(2a)') "...RM spectrum saved to ",TRIM(fileroot)//"--RM_spectrum"

  ! save array modes U(p,ndof)

  DO l = 1, ndof
    WRITE (filename,'(2a,i4.4)') TRIM(fileroot),"--mode",l
    OPEN (1,file=filename,status="unknown")
    DO k = 1, p
      WRITE (1,*) k,U(k,l)
    END DO
    CLOSE (1)
    WRITE (*,'(a,i5,2a)') "...array mode",l," saved to ",TRIM(filename)
  END DO

  ! save modal representers rho_mu(n,ndof)

  do l = 1, ndof
    WRITE (filename,'(2a,i4.4)') TRIM(fileroot),"--modrep",l
    OPEN (1,file=filename,status="unknown")
    DO k = 1, p
      WRITE (1,*) k,rho_mu(k,l)
    END DO
    CLOSE (1)
    WRITE (*,'(a,i5,2a)') "...modal representer",l," saved to ",TRIM(filename)
  END DO

  END SUBROUTINE saveResults

END PROGRAM example_arm