! 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

  REAL( kind = dp ) ::  CProxPar(NPT, 1:5, NLOCS)
  
  ! copy the chose ADM to the structure in the TUPDATER_DATA module
  ! no problem doing this here, as this is not used in the FORTRAN
  ! code at any time.
  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_mv(tau2Pprior, ProxNum, MInst)
  
  USE random
  USE TUPDATER_DATA

  IMPLICIT NONE
  
  INTEGER, INTENT(IN)           ::  ProxNum
  REAL( kind = dp ), INTENT(IN)  ::  tau2Pprior(*), MInst(1:NLOCS)

  REAL( kind = dp )  ::  TMatDiff(1:NLOCS, 1:NYEARS)
  REAL( kind = dp )  ::  TauBeta(1:NLOCS)
  REAL               ::  TauAlpha(1:NLOCS)
  INTEGER            ::  I
  !
  !UPDATES the variance of the proxy 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 = SPREAD(CProx(ProxNum, 4, :),2,NYEARS-1) * T_MATRIX(1:NLOCS, 2:NYEARS) + &
        SPREAD(CProx(ProxNum, 3, :),2,NYEARS-1) - &
        DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 2:NYEARS) +  &
        SPREAD(CProx(ProxNum, 5, :),2,NYEARS-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 = SPREAD(CProx(ProxNum, 4, :),2,NYEARS) * T_MATRIX(1:NLOCS, 1:NYEARS) + &
        SPREAD(CProx(ProxNum, 3, :),2,NYEARS) - &
        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)/2.0_dp + tau2Pprior(2)
  !can then calculate the first parameter, alpha, for the posterior inv-gamma dist.:
  TauAlpha = MInst(1:NLOCS) / 2.0_dp + tau2Pprior(1)
  
  !make the draw:
  DO I=1,NLOCS 
    CProx(ProxNum, 2, I) =  TauBeta(I) / random_gamma( TauAlpha( I) , .TRUE.)
  ENDDO
  WHERE ( MInst == 0 )
    CProx(ProxNum, 2, :) = -99
  END WHERE
  WHERE ( CProx(ProxNum, 2, :) > 5)
    CProx(ProxNum, 2, :) = 5
  END WHERE
  
END SUBROUTINE Tau2P_Updater_mv


SUBROUTINE Beta2_Updater_mv( 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
  REAL( kind = dp )  ::  Beta2Prior(*)
  REAL( kind = dp )  ::  TMatDiff(1:NLOCS, 1:NYEARS)
  REAL( kind = dp )  ::  Beta2Var(1:NLOCS), Beta2Mean(1:NLOCS)

  !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.0_dp / (1.0_dp / Beta2Prior(2) + SUM(TMatDiff**2,2)/CProx(ProxNum, 2, :))
  WHERE( SUM(TMatDiff,2) == 0) 
    Beta2Var = 0 
  END WHERE
  !WRITE( *, '(10F4.1)') Beta2Var
  
  !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 = ( - SPREAD(CProx(ProxNum, 4, :),2,NYEARS-1) * T_MATRIX(1:NLOCS, 2:NYEARS) + &
      DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 2:NYEARS) - &
      SPREAD( CProx(ProxNum, 3, :),2,NYEARS-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.0_dp/CProx(ProxNum, 2, :) * SUM(TMatDiff,2) )
  
  !make the draw:
  CALL getzig( CProx(ProxNum, 5, :), NLOCS)
  CProx(ProxNum, 5, :) = Beta2Mean + CProx(ProxNum, 5, :) * Beta2Var**.5

END SUBROUTINE Beta2_Updater_mv

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

  IMPLICIT NONE
  
  INTEGER           ::  ProxNum
  REAL( kind = dp )  ::  Beta1Prior(*)
  REAL( kind = dp )  ::  TMatDiff(1:NLOCS, 1:NYEARS)
  REAL( kind = dp )  ::  Beta1Var(1:NLOCS), Beta1Mean(1:NLOCS)

  !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.0_dp / (1.0_dp / Beta1Prior(2) + SUM(TMatDiff**2,2)/CProx(ProxNum, 2, :))
  
  !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) 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) - &
        SPREAD( CProx(ProxNum, 3, :), 2, NYEARS-1) - &
        SPREAD(CProx(ProxNum, 5, :), 2, NYEARS-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) -&
        SPREAD(CProx(ProxNum, 3, :), 2,NYEARS) )
    END WHERE
  END IF
  
  !now we need to sum all of the Non NaN entries:
  
  Beta1Mean = Beta1Var*( Beta1Prior(1)/Beta1Prior(2) + 1.0_dp/CProx(ProxNum, 2, :) * SUM(TMatDiff,2) )
  
  !make the draw:
  CALL getzig( CProx(ProxNum, 4, :), NLOCS)
  CProx(ProxNum, 4, :) = Beta1Mean + CProx(ProxNum, 4, :) * Beta1Var**.5
  WHERE( SUM(TMatDiff,2) == 0) 
    CProx(ProxNum, 4, :) = 0
  END WHERE
   
END SUBROUTINE Beta1_Updater_mv

SUBROUTINE Beta0_Updater_mv( 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
  REAL( kind = dp )  ::  Beta0Prior(*), MProx(1:NLOCS)
  REAL( kind = dp )  ::  TMatDiff(1:NLOCS, 1:NYEARS-1)
  REAL( kind = dp )  ::  Beta0Var(1:NLOCS), Beta0Mean(1:NLOCS)

  !the posterior variance is easy:
  Beta0Var = 1.0_dp/( 1.0_dp/Beta0prior(2) + MProx / CProx(ProxNum, 2, :) )
  WHERE (MProx == 0)
    Beta0Var = 0
  END WHERE
  
  !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) 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) - &
        SPREAD(CProx(ProxNum, 4, :), 2, NYEARS-1) *T_MATRIX(1:NLOCS, 2:NYEARS) - &
        SPREAD(CProx(ProxNum, 5, :), 2, NYEARS-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 )
      TMatDiff = 0;
    ELSEWHERE  
      TMatDiff = DATA_ALL(NLOCS*ProxNum + 1:NLOCS*(ProxNum + 1), 2:NYEARS) - &
        SPREAD(CProx(ProxNum, 4, :), 2, NYEARS-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.0_dp / CProx(ProxNum, 2, :) * &
    SUM(TmatDiff,2))
  
  !make the draw:
  CALL getzig( CProx(ProxNum, 3, :), NLOCS)
  CProx(ProxNum, 3, :) = Beta0Mean + CProx(ProxNum, 3, :) * Beta0Var**.5
  WHERE( SUM(TMatDiff,2) == 0) 
    CProx(ProxNum, 3, :) = 0
  END WHERE
  
END SUBROUTINE Beta0_Updater_mv


