oi_hor_extrapol_surf.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 oi_hor_extrapol_surf(NDIM,PLAT_IN,PLON_IN,PFIELD_IN, &
                                         plat,plon,pfield,ointerp,pzs,ndim2)  
!     ###################################################################
!
!!**** *OI_HOR_EXTRAPOL_SURF* extrapolate a surface field
!!
!!    PURPOSE
!!    -------
!!
!!
!!    METHOD
!!    ------
!!
!!    For each point to interpolate, the nearest valid point value is set.
!!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!
!!    V. Masson          Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original     01/12/98
!!   V. Masson     01/2004 extrapolation in latitude and longitude
!!   J.-F. Mahfouf 03/2010 adaptation for OI soil analysis 
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
USE modd_surf_par,   ONLY : xundef
USE modd_csts,       ONLY : xpi
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
USE modi_abor1_sfx
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
INTEGER,               INTENT(IN)     :: ndim  ! dimension of arrays
REAL,   DIMENSION(NDIM),  INTENT(IN)     :: plat_in  ! input lat. of each grid mesh.
REAL,   DIMENSION(NDIM),  INTENT(IN)     :: plon_in  ! input lon. of each grid mesh.
REAL,   DIMENSION(NDIM),  INTENT(IN)     :: pfield_in! input field on grid mesh
REAL,   DIMENSION(NDIM),  INTENT(IN)     :: plat     ! latitude of each grid mesh.
REAL,   DIMENSION(NDIM),  INTENT(IN)     :: plon     ! longitude of each grid mesh.
REAL,   DIMENSION(NDIM),  INTENT(INOUT)  :: pfield   ! field on grid mesh
LOGICAL,DIMENSION(NDIM),  INTENT(IN)     :: ointerp  ! .true. where physical value is needed
REAL,   DIMENSION(NDIM), OPTIONAL, INTENT(IN) :: pzs      ! surface height
INTEGER,                 OPTIONAL, INTENT(IN) :: ndim2 ! Optional subdomain to search in
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
REAL     :: zlat    ! latitude of point to define
REAL     :: zlon    ! longitude of point to define
REAL     :: zdist   ! current distance to valid point (in lat/lon grid)
REAL     :: zfield  ! current found field value
REAL     :: zndist  ! smallest distance to valid point
REAL     :: zcosla  ! cosine of latitude
REAL     :: zzs_out ! altitude of nearest grid point
REAL,PARAMETER :: zlimmax = 1.  ! Maximum distance allowed (in degrees)
!
INTEGER  :: ji    ! loop index on points
INTEGER  :: jisc  ! loop index on valid points
INTEGER  :: jisc1,jisc2,jzlimcnt
!
REAL     :: zlonsc, zdlat, zdlon, zconv, zr_earth
REAL(KIND=JPRB) :: zhook_handle
!
LOGICAL :: lndim2
!
! Earth radius
!
IF (lhook) CALL dr_hook('OI_HOR_EXTRAPOL_SURF',0,zhook_handle)
zr_earth = 6371598.0 
!
! Angle conversion factor
!  
zconv = xpi/180.0
!-------------------------------------------------------------------------------
!
!*    3.     No data point
!            -------------
!
IF (count(pfield_in(:)/=xundef)==0 .AND. lhook) CALL dr_hook('OI_HOR_EXTRAPOL_SURF',1,zhook_handle)
IF (count(pfield_in(:)/=xundef)==0) RETURN
!
!-------------------------------------------------------------------------------
!
!*      4.   Loop on points to define
!            ------------------------
!
jzlimcnt = 0
!
DO ji=1,ndim
  IF (pfield(ji)/=xundef) cycle
  IF (.NOT. ointerp(ji))  cycle
!
!*      4.1  initialisation
!            --------------
!
  zndist=xundef
  zlat=plat(ji)
  zlon=plon(ji)
  zfield=pfield(ji)
  zcosla=cos(zlat*zconv)
  IF (present(pzs)) zzs_out=pzs(ji)
  IF ( present(ndim2)) THEN
    jisc1=max((ji-ndim2),1)
    jisc2=min((ji+ndim2),ndim)
    lndim2=.true.
  ELSE
    jisc1=1
    jisc2=ndim
    lndim2=.false.
  ENDIF
!
!*      4.2  extrapolation with nearest valid point
!            --------------------------------------
!
  DO jisc=jisc1,jisc2
    IF (pfield_in(jisc)/=xundef) THEN
      zlonsc = plon_in(jisc)
      IF (zlonsc-zlon> 180.) zlonsc = zlonsc - 360.0
      IF (zlonsc-zlon<-180.) zlonsc = zlonsc + 360.0
      zdlat = (plat_in(jisc)-zlat)*zconv
      zdlon = (zlonsc-zlon)*zconv
      zdist = zdlat*zdlat + zdlon*zdlon*zcosla*zcosla
      IF (zdist<=zndist) THEN
        zfield=pfield_in(jisc)
        IF (present(pzs)) zzs_out=pzs(jisc)
        zndist=zdist
      END IF
    END IF
  END DO

  ! Check if we got values
  IF ( zndist  == xundef ) THEN
    CALL abor1_sfx("Extrapolated point is undefined! No nearby point found.")      
  ELSEIF ( zndist > (zlimmax*zconv) ) THEN
    IF ( lndim2 ) &
    & CALL abor1_sfx("Distance to extrapolated point is to large. Increase ZLIMMAX or NDIM2")     
    jzlimcnt = jzlimcnt + 1
  ENDIF
  IF (present(pzs)) THEN
    pfield(ji) = zfield + (zzs_out - pzs(ji))*0.0065
  ELSE
    pfield(ji) = zfield
  ENDIF  

END DO

IF ( jzlimcnt > 0 ) THEN
  print *,'Points with extrapolation distance > ',zlimmax,' degrees are ',jzlimcnt
ENDIF

IF (lhook) CALL dr_hook('OI_HOR_EXTRAPOL_SURF',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE oi_hor_extrapol_surf