PDAF_diag_crps.F90 Source File


Source Code

! Copyright (c) 2012-2024 Lars Nerger, lars.nerger@awi.de
!
! This routine 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 code 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/>.
!
!$Id$

!> Computation of CRPS
!!
!! This routine computes the continuous ranked probability
!! score (CRPS) and its decomposition into uncertainty and
!! resolution: CRPS = RELI + RESOL. In addition the uncertainty
!! is computed.
!! A perfectly reliable system gives RELI=0.
!! An informative system gives RESOL << UNCERT.
!!
!! The computation follows H. Hersbach, Weather and Forecasting
!! 15(2000) 599-570. Here, RESOL is equivalent to CPRS_pot.
!!
!! __Revision history:__
!! * 2021-05 - Lars Nerger - Initial code based on sangoma_ComputeCRPS
!! * Later revisions - see repository log
!!
SUBROUTINE PDAF_diag_CRPS(dim, dim_ens, element, oens, obs, &
     CRPS, reli, resol, uncert, status)!

  IMPLICIT NONE

! *** Arguments ***
  INTEGER, INTENT(in) :: dim                !< PE-local state dimension
  INTEGER, INTENT(in) :: dim_ens            !< Ensemble size
  INTEGER, INTENT(in) :: element            !< ID of element to be used
       !< If element=0, mean values over all elements are computed
  REAL, INTENT(in)    :: oens(dim, dim_ens) !< State ensemble
  REAL, INTENT(in)    :: obs(dim)           !< State ensemble
  REAL, INTENT(out)   :: CRPS               !< CRPS
  REAL, INTENT(out)   :: reli               !< Reliability
  REAL, INTENT(out)   :: resol              !< resolution
  REAL, INTENT(out)   :: uncert             !< uncertainty
  INTEGER, INTENT(out) :: status            !< Status flag (0=success)

! *** local variables ***
  INTEGER :: i, elem, istart                ! Counters
  INTEGER :: imax
  REAL :: dinv
  REAL :: oneobs
  REAL :: gi, oi, pval
  REAL, ALLOCATABLE :: oneens(:)
  REAL, ALLOCATABLE :: allobs(:)
  REAL, ALLOCATABLE :: c_a(:), c_b(:)


! ********************
! *** Compute CRPS ***
! ********************

  ! Set number of element over which CPRS is computed
  IF (element==0) THEN
     imax = dim
     dinv = 1.0/REAL(dim)
     istart = 1
  ELSEIF (element<=dim) THEN
     imax = element
     istart = element
     dinv = 1.0
  ELSE
     imax = 0
     status = 100
     dinv = 1.0
     istart = 1
  END IF

  ALLOCATE(oneens(dim_ens))
  ALLOCATE(c_a(0:dim_ens))
  ALLOCATE(c_b(0:dim_ens))

  c_a = 0.0
  c_b = 0.0

  ! Loop over elements
  DO elem = istart, imax

     ! Get observation for current element
     oneobs = obs(elem)

     ! Get sorted ensemble for current element
     DO i = 1, dim_ens
        oneens(i) = oens(elem, i)
     END DO
     CALL PDAF_sisort(dim_ens, oneens)

     IF (oneobs < oneens(1)) THEN
        ! Case 1: obs < all ensemble members
        c_a(0) = c_a(0) + dinv*(oneens(1) - oneobs)
        c_b(0) = c_b(0) + dinv
     ELSEIF (oneobs > oneens(dim_ens)) THEN
        ! Case 2: obs > all ensemble members
        c_a(dim_ens) = c_a(dim_ens) + dinv
        c_b(dim_ens) = c_b(dim_ens) + dinv*(oneobs - oneens(dim_ens))
     END IF
     DO i = 1, dim_ens-1
        c_a(i) = c_a(i) + dinv*MAX(oneens(i+1) - MAX(oneobs, oneens(i)), 0.0)
        c_b(i) = c_b(i) + dinv*MAX( MIN(oneobs, oneens(i+1)) - oneens(i), 0.0)
     END DO
  END DO

  haveobs: IF (imax>0) THEN
     
     ! Reinterpretation of c_a(dim_ens) and c_b(0)
     IF (c_b(0) /= 0.0) c_b(0) = c_a(0) * (1.0/c_b(0) - 1.0)
     IF (c_a(dim_ens) /= 0.0) c_a(dim_ens) = c_b(dim_ens) * (1.0/c_a(dim_ens) - 1.0)

     ! Calculate untertainty
     ALLOCATE(allobs(dim))
     allobs = obs

     CALL PDAF_sisort(dim, allobs)
     uncert = 0
     pval = 0.0

     DO i = 1, dim-1
        pval = pval + dinv
        uncert = uncert + (allobs(i+1) - allobs(i)) * pval*(1.0-pval)
     END DO

     ! Complete computation of CPRS, reliability and resolution
     crps = 0.0
     reli = 0.0
     resol = 0.0

     oi = 0.0
     DO i = 1, dim_ens
        gi = c_a(i) + c_b(i)
        IF (gi /= 0.0) THEN
           oi = c_a(i) / gi
        ELSE
           oi = 0.0
        END IF
        pval = REAL(i) / REAL(dim_ens)

        crps = crps + c_b(i)*pval*pval + c_a(i)*(1.0-pval)*(1.0-pval)
        reli = reli + gi * (oi-pval)*(oi-pval)
        resol = resol + gi * oi * (1.0-oi)

     END DO

  ELSE
     crps = 0.0
     reli = 0.0
     resol = 0.0
     uncert = 0.0
  END IF haveobs


! ****************
! *** Clean up ***
! ****************

  DEALLOCATE(oneens, c_a, c_b)

END SUBROUTINE PDAF_diag_CRPS


!--------------------------------------------------------
!> Sorting routine 
!!
!! Sorts a vector veca(1:n) into ascending numerical order, by
!! straight insertion.
!!
!! For large vectors, this routine will not be efficient.
!!
SUBROUTINE PDAF_sisort(n, veca)

  IMPLICIT NONE

  INTEGER, INTENT(in) :: n 
  REAL, INTENT(inout) :: veca(n)

  INTEGER :: i, j, k 
  REAL :: tmpa
  LOGICAL :: eflag

  DO j = 2, n

     eflag = .FALSE.

     tmpa = veca(j) 

     sortloop: DO i = j-1, 1, -1 
        k = i

        IF(veca(i) <= tmpa) THEN
           eflag = .TRUE.
           EXIT sortloop
        END IF

        veca(i+1) = veca(i) 
     ENDDO sortloop

     IF (.NOT.eflag) k=0

     veca(k+1) = tmpa 

  ENDDO 

END SUBROUTINE PDAF_sisort