average2_cti.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 average2_cti
!     #######################
!
!!**** *AVERAGE2_CTI* computes the topo index stats
!!
!!    PURPOSE
!!    -------
!!
!!    METHOD
!!    ------
!!   
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!
!!    B. Decharme         Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    06/2009
!!
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
USE modd_surf_par, ONLY : xundef
USE modd_pgdwork,       ONLY : nsize, xsumval, xprec, &
                               xmean_work, xstd_work, xskew_work, &
                               xmin_work, xmax_work 
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
!*    0.2    Declaration of other local variables
!            ------------------------------------
!
REAL, DIMENSION(SIZE(NSIZE,1)) :: ZSIZE
!
REAL :: ZINT
INTEGER :: JI
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!----------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('AVERAGE2_CTI',0,zhook_handle)
zsize(:)=REAL(nsize(:,1))
!
WHERE (nsize(:,1)>=36)
!
!----------------------------------------------------------------------------
!
!*    1.     Mean CTI
!            --------------
!
  xmean_work(:) = xsumval(:,1)/zsize(:)
!
!-------------------------------------------------------------------------------
!
!*    2.     Standard deviation
!            ------------------
!
  WHERE (xmax_work(:)-xmin_work(:)>=1.0) 
    xstd_work(:) = sqrt( max(0.,xsumval(:,2)/nsize(:,1) - xmean_work(:)*xmean_work(:)) )
  ELSEWHERE
    xstd_work(:) = 0.0
  END WHERE
!
!-------------------------------------------------------------------------------
!
!*    3.     Skewness
!            --------
!
  WHERE(xstd_work(:)>0.0)
!          
        xskew_work(:) = xsumval(:,3)-zsize(:)*xmean_work(:)*xmean_work(:)*xmean_work(:) &
                       -3.0*zsize(:)*xmean_work(:)*xstd_work(:)*xstd_work(:)  
!
        xskew_work(:) = xskew_work(:)/(zsize(:)*xstd_work(:)*xstd_work(:)*xstd_work(:))
!
  END WHERE
!
END WHERE
!----------------------------------------------------------------------------
!
DO ji = 1,SIZE(xmean_work,1)

  IF (xmean_work(ji)/=xundef) THEN

    zint = aint(xmean_work(ji))
    IF (xmean_work(ji)/=zint) &
      xmean_work(ji) = zint + anint((xmean_work(ji)-zint)*xprec)/xprec

    zint = aint(xmin_work(ji))
    IF (xmin_work(ji)/=zint) &
      xmin_work(ji) = zint + anint((xmin_work(ji)-zint)*xprec)/xprec

    zint = aint(xmax_work(ji))
    IF (xmax_work(ji)/=zint) &      
      xmax_work(ji) = zint + anint((xmax_work(ji)-zint)*xprec)/xprec

    zint = aint(xstd_work(ji))
    IF (xstd_work(ji)/=zint) &      
    xstd_work(ji) = zint + anint((xstd_work(ji)-zint)*xprec)/xprec

    zint = aint(xskew_work(ji))
    IF (xskew_work(ji)/=zint) &      
    xskew_work(ji) = zint + anint((xskew_work(ji)-zint)*xprec)/xprec
!
  ENDIF
ENDDO
!
IF (lhook) CALL dr_hook('AVERAGE2_CTI',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE average2_cti