! You need to typecast only the integers when calling the Fortran subroutines, 
! as R seems to initialise all numbers to float (well, "numeric"), and those seem
! to be handled well enough.
!
! 2014-05-28    restructuring of code, CProx is now global in the main module
!               and also started specifying INTENTs for the variables.
! 2012-06-26	newly created file, moved proxy stuff from T_Updater.f95

!SUBROUTINE getCProx( CProxPar)
!
!  USE TUPDATER_DATA
!
!  DOUBLE PRECISION ::  CProxPar(NPT, *)
!  
!  WHERE( CProxPar(:,1) == 20)
!    CProx(:,5) = CProxPar(:,5)
!  ENDWHERE
!
!  CProxPar(1:NPT,1:5) = CProx(1:NPT,1:5)
!
!END SUBROUTINE getCProx

SUBROUTINE Tau2P_Updater(tau2Pprior, ProxNum, MInst)
  
  USE random
  USE TUPDATER_DATA

  IMPLICIT NONE
  
  INTEGER, INTENT(IN)           ::  ProxNum
  DOUBLE PRECISION, INTENT(IN)  ::  tau2Pprior(*), MInst(1:NLOCS)

  DOUBLE PRECISION  ::  TMatDiff(1:NLOCS, 1:NYEARS)
  REAL              ::  TauAlpha, TauBeta
  !
  !UPDATES the variance of the instrumental error in the main BARCAST code
  !need to find the sum of the squared residuals between the Inst
  !observations and the corresponding Temp values:
  !extract the Inst part of the T,D matrices:
  !also get rid of the time=0 value from the temp mat:
  !we don't actually need the HH_Select matrix: just take the difference
  !between the Temperature matrix and ther Inst Obs matrix:: wherever there
  !is a NaN in the Data matrix, there will be a NaN in the difference. Then
  !find all non NaN entries, square them, add them up.

  IF (CProx(ProxNum, 1, 1) == 1.0) THEN
    WHERE( DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 2:NYEARS) == -99 .OR. &
           DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 1:NYEARS-1) == -99)
      TMatDiff = 0;
    ELSEWHERE  
      TMatDiff = CProx(ProxNum, 4,1) * T_MATRIX(1:NLOCS, 2:NYEARS) + CProx(ProxNum, 3,1) - &
        DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 2:NYEARS) +  &
        CProx(ProxNum, 5,1) * DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 1:NYEARS-1);
    END WHERE
  ELSE
    WHERE( DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 1:NYEARS) == -99)
      TMatDiff = 0;
    ELSEWHERE  
      TMatDiff = CProx(ProxNum, 4,1) * T_MATRIX(1:NLOCS, 1:NYEARS) + CProx(ProxNum, 3,1) - &
        DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 1:NYEARS)
    END WHERE
  END IF
  
  ! sum of the squares of the residuals and calculate the postertior beta parameter:
  TauBeta = SUM(TMatDiff**2)/2.0D0 + tau2Pprior(2)
  !can then calculate the first parameter, alpha, for the posterior inv-gamma dist.:
  TauAlpha = SUM(MInst) / 2.0D0 + tau2Pprior(1)
  
  !make the draw:
  CProx(ProxNum, 2, :) =  TauBeta / random_gamma( TauAlpha , .TRUE.)
  !write(*,*) CProx
  
END SUBROUTINE Tau2P_Updater


SUBROUTINE Beta2_Updater( Beta2Prior, ProxNum)
  !
  ! UPDATES the scaling constant for the past proxy state in the proxy measurement
  ! equation in the main BARCAST code.
  
  USE random
  USE TUPDATER_DATA

  IMPLICIT NONE
  
  INTEGER           ::  ProxNum
  DOUBLE PRECISION  ::  Beta2Prior(*)
  DOUBLE PRECISION  ::  TMatDiff(1:NLOCS, 1:NYEARS)
  DOUBLE PRECISION  ::  Beta2Var, Beta2Mean

  !for the posterior variance, need the sum of squares of ALL Proxy
  !values which correspond to proxy observation in the next year.
  !So: extract the part of the data matrix for which there are prox obs
  !in the next year:
  !must get rid of the time=0 value from the temp mat:
  !find the elements of the proxy data matrix which are not NaNs:
  
  WHERE( DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 1:NYEARS) == -99 )
    TMatDiff = 0;
  ELSEWHERE  
    TMatDiff = DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 1:NYEARS)
  END WHERE
  !now the posterior variance is easy:
  Beta2Var = 1.0D0 / (1.0D0 / Beta2Prior(2) + SUM(TMatDiff**2)/CProx(ProxNum, 2,1))
  
  !For the posterior mean:
  !Remove the additive constant from each Proxy obs, multiple each by the
  !corresponding instrumental value, and sum. 
  
  TMatDiff = 0;
  WHERE( DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 2:NYEARS) == -99 .OR. &
         DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 1:NYEARS-1) == -99 )
    TMatDiff = 0;
  ELSEWHERE  
    TMatDiff = ( - CProx(ProxNum, 4,1) * T_MATRIX(1:NLOCS, 2:NYEARS) + &
      DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 2:NYEARS) - CProx(ProxNum, 3,1)) * &
        DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 1:NYEARS-1)
  END WHERE
  
  !now we need to sum all of the Non NaN entries:
  
  Beta2Mean = Beta2Var*( Beta2Prior(1)/Beta2Prior(2) + 1.0D0/CProx(ProxNum, 2,1) * SUM(TMatDiff) )
  
  ! make the draw and save it in _all_ NLOCS columns for this proxy -> no problem to 
  ! save it for each step...
  CALL getzig( CProx(ProxNum, 5,1), 1)
  CProx(ProxNum, 5, :) = Beta2Mean + CProx(ProxNum, 5,1) * Beta2Var**.5

END SUBROUTINE Beta2_Updater

SUBROUTINE Beta1_Updater( Beta1Prior, ProxNum)
  !
  ! UPDATES the scaling constant in the proxy measurement equation in the main
  ! BARCAST code
  !
  ! NORMAL
  ! 
  
  USE random
  USE TUPDATER_DATA

  IMPLICIT NONE
  
  INTEGER           ::  ProxNum
  DOUBLE PRECISION  ::  Beta1Prior(*)
  DOUBLE PRECISION  ::  TMatDiff(1:NLOCS, 1:NYEARS)
  DOUBLE PRECISION  ::  Beta1Var, Beta1Mean

  !for the posterior variance, need the sum of squares of ALL temperature
  !values which correspond to proxy observation.
  !So: extract the part of the T matrix for which there are prox obs:
  !must get rid of the time=0 value from the temp mat:
  !find the elements of the proxy data matrix which are not NaNs:
  
  WHERE( DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 1:NYEARS) == -99 )
    TMatDiff = 0;
  ELSEWHERE  
    TMatDiff = T_MATRIX(1:NLOCS, 1:NYEARS)
  END WHERE
  !now the posterior variance is easy:
  Beta1Var = 1.0D0 / (1.0D0 / Beta1Prior(2) + SUM(TMatDiff**2)/CProx(ProxNum, 2,1))
  
  !For the posterior mean:
  !Remove the additive constant from each Proxy obs, multiple each by the
  !corresponding instrumental value, and sum. 
  
  TMatDiff = 0;
  IF (CProx(ProxNum, 1,1) == 1.0) THEN
    WHERE( DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 2:NYEARS) == -99 .OR. &
           DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 1:NYEARS-1) == -99)
      TMatDiff = 0;
    ELSEWHERE  
      TMatDiff = T_MATRIX(1:NLOCS, 2:NYEARS) * &
        (DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 2:NYEARS) - CProx(ProxNum, 3,1) - &
          CProx(ProxNum, 5,1) * DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 1:NYEARS-1) )
    END WHERE
  ELSE
    WHERE( DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 1:NYEARS) == -99)
      TMatDiff = 0;
    ELSEWHERE  
      TMatDiff = T_MATRIX(1:NLOCS, 1:NYEARS) * &
        (DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 1:NYEARS) - CProx(ProxNum, 3,1))
    END WHERE
  END IF
  
  !now we need to sum all of the Non NaN entries:
  
  Beta1Mean = Beta1Var*( Beta1Prior(1)/Beta1Prior(2) + 1.0D0/CProx(ProxNum, 2,1) * SUM(TMatDiff) )
  
  !make the draw:
  CALL getzig( CProx(ProxNum, 4,1), 1)
  CProx(ProxNum, 4, :) = Beta1Mean + CProx(ProxNum, 4,1) * Beta1Var**.5
!  CProx(ProxNum, 4) = .5
!  CProx(ProxNum, 4) = .5

END SUBROUTINE Beta1_Updater

SUBROUTINE Beta0_Updater( Beta0prior, ProxNum, MProx)
  !
  !UPDATES the Beta_0 parameter, constant adjustment to proxies, in main
  !BARCAST code
  !
  ! NORMAL
  !
  
  USE random
  USE TUPDATER_DATA

  IMPLICIT NONE
  
  INTEGER           ::  ProxNum
  DOUBLE PRECISION  ::  Beta0Prior(*), MProx(1:NLOCS)
  DOUBLE PRECISION  ::  TMatDiff(1:NLOCS, 1:NYEARS)
  DOUBLE PRECISION  ::  Beta0Var, Beta0Mean

  !the posterior variance is easy:
  Beta0Var = 1.0D0/( 1.0D0/Beta0prior(2) + SUM(MProx) / CProx(ProxNum, 2,1) )
  
  !extract the part of the T matrix for which there are prox obs:
  !must get rid of the time=0 value from the temp mat:
  !find the difference between the Temp Mat, and the adjusted data mat.
  !As before, the NaNs in the Data Matrix will carry over:
  
  TMatDiff = 0;
  IF (CProx(ProxNum, 1,1) == 1.0) THEN
    WHERE( DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 2:NYEARS) == -99 .OR. &
           DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 1:NYEARS-1) == -99 )
      TMatDiff = 0;
    ELSEWHERE  
      TMatDiff = DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 2:NYEARS) - &
        CProx(ProxNum, 4,1)*T_MATRIX(1:NLOCS, 2:NYEARS) - &
        CProx(ProxNum, 5,1) * DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 1:NYEARS-1)
    END WHERE
  ELSE
    WHERE( DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 2:NYEARS) == -99 .OR. &
         DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 1:NYEARS-1) == -99 )
      TMatDiff = 0;
    ELSEWHERE  
      TMatDiff = DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 2:NYEARS) - &
        CProx(ProxNum, 4,1)*T_MATRIX(1:NLOCS, 2:NYEARS)
    END WHERE
  END IF
  
  !and we can find the posterior mean:
  Beta0Mean = Beta0Var * ( Beta0prior(1)/Beta0prior(2) + 1.0D0 / CProx(ProxNum, 2,1) * SUM(TmatDiff))
  
  !make the draw:
  CALL getzig( CProx(ProxNum, 3,1), 1)
  CProx(ProxNum, 3, :) = Beta0Mean + CProx(ProxNum, 3,1) * Beta0Var**.5
!  CProx(ProxNum, 3) = 0
!  CProx(ProxNum, 3) = 0
  
END SUBROUTINE Beta0_Updater


