hor_extrapol_surf_cheap.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 hor_extrapol_surf_cheap(KLUOUT,HCOORTYPE,PLAT_IN,PLON_IN,PFIELD_IN, &
                                         plat,plon,pfield,ointerp)  
!     ###################################################################
!
!!**** *HOR_EXTRAPOL_SURF_CHEAP* 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
!!     M. Jidane    11/2013 add OpenMP directives
!!     L. Raynaudl  06/2014 cheap version : nearest point looked for in a
!!                          neighbourhood
!
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
USE modd_surf_par,   ONLY : xundef
USE modd_csts,       ONLY : xpi, xradius
USE modn_prep_surf_atm, ONLY : nhalo_prep
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
INTEGER,               INTENT(IN)     :: kluout   ! output listing logical unit
 CHARACTER(LEN=4), INTENT(IN) :: hcoortype
REAL,   DIMENSION(:),  INTENT(IN)     :: plat_in  ! input lat. of each grid mesh.
REAL,   DIMENSION(:),  INTENT(IN)     :: plon_in  ! input lon. of each grid mesh.
REAL,   DIMENSION(:,:),  INTENT(IN)     :: pfield_in! input field on grid mesh
REAL,   DIMENSION(:),  INTENT(IN)     :: plat     ! latitude of each grid mesh.
REAL,   DIMENSION(:),  INTENT(IN)     :: plon     ! longitude of each grid mesh.
REAL,   DIMENSION(:),  INTENT(INOUT)  :: pfield   ! field on grid mesh
LOGICAL,DIMENSION(:),  INTENT(IN)     :: ointerp  ! .true. where physical value is needed
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
INTEGER  :: ino   ! output  array size
!
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
!
INTEGER  :: ji,jii,jj,ici,icj,ici_min,ici_max,icj_min,icj_max,icompt    ! loop index on points
REAL     :: zlonsc! longitude of valid point
LOGICAL  :: glalo ! flag true is second coordinate is a longitude or pseudo-lon.
                  !      false if metric coordinates
!
REAL(KIND=JPRB) :: zrad ! conversion degrees to radians
!
REAL(KIND=JPRB) :: zhook_handle, zhook_handle_omp
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('HOR_EXTRAPOL_SURF_CHEAP',0,zhook_handle)
!
ino = SIZE(pfield,1)
!
WHERE (.NOT. ointerp(:)) pfield(:) = xundef
!
!-------------------------------------------------------------------------------
!
glalo = hcoortype=='LALO'
!
!-------------------------------------------------------------------------------
!
!*    3.     No data point
!            -------------
!
IF (count(pfield_in(:,:)/=xundef)==0 .AND. lhook) CALL dr_hook('HOR_EXTRAPOL_SURF_CHEAP',1,zhook_handle)
IF (count(pfield_in(:,:)/=xundef)==0) RETURN
!
!-------------------------------------------------------------------------------
!
!*      4.   Loop on points to define
!            ------------------------
!
zrad=xpi/180.0_jprb
!
!$OMP PARALLEL DO SCHEDULE(DYNAMIC,1) PRIVATE(JI,JII,JJ,ZLAT,ZLON,ZFIELD,ZCOSLA,ZLONSC,ZDIST,ZNDIST,ZHOOK_HANDLE_OMP)
DO ji=1,ino

  IF (pfield(ji)/=xundef) cycle
  IF (.NOT. ointerp(ji))  cycle
!
!*      4.1  initialisation
!            --------------
!
  IF (lhook) CALL dr_hook('HOR_EXTRAPOL_SURF_CHEAP OMP',0,zhook_handle_omp)
  zndist=1.e20
  zlat=plat(ji)
  zlon=plon(ji)
  zfield=pfield(ji)
  zcosla=cos(zlat*zrad)

  ! Locate surrounding input points around the interpolated output point
  ! Size of the neighbourhood is given by NHALO_PREP
  ici = 1
  DO jj=SIZE(plon_in),1,-1
    IF (plon_in(jj)<=zlon) THEN
      ici = jj
      EXIT
    ENDIF
  ENDDO
  icj = 1
  DO jj=SIZE(plat_in),1,-1
    IF (plat_in(jj)<=zlat) THEN
      icj = jj
      EXIT
    ENDIF
  ENDDO

  ici_min=max(ici-nhalo_prep,1)
  ici_max=min(ici+nhalo_prep,SIZE(plon_in))
  icj_min=max(icj-nhalo_prep,1)
  icj_max=min(icj+nhalo_prep,SIZE(plat_in))
!
!*      4.2  extrapolation with nearest valid point
!            --------------------------------------

  ! The nearest point is looked for in a neighbourhood
  ! This helps reducing the cost
  icompt=0
  DO jj=icj_min,icj_max
    DO jii=ici_min,ici_max
      IF (pfield_in(jii,jj)/=xundef) THEN
        icompt=icompt+1
        zlonsc = plon_in(jii)
        IF (glalo) THEN
          IF (zlonsc-zlon> 180.) zlonsc = zlonsc - 360.
          IF (zlonsc-zlon<-180.) zlonsc = zlonsc + 360. 
          zdist= (plat_in(jj)-zlat) ** 2 + ((zlonsc-zlon)*zcosla) ** 2
        ELSE
          zdist= (plat_in(jj)-zlat) ** 2 + (zlonsc-zlon) ** 2
        END IF
        IF (zdist<=zndist) THEN
          zfield=pfield_in(jii,jj)
          zndist=zdist
        END IF
      END IF
    ENDDO
  ENDDO
  pfield(ji) = zfield
  IF (icompt==0) THEN
    WRITE(*,*) 'NO EXTRAPOLATION : INCREASE YOUR HALO_PREP IN NAM_PREP_SURF_ATM'
    CALL abor1_sfx("HOR_EXTRAPOL_SURF_CHEAP: INCREASE YOUR HALO_PREP IN NAM_PREP_SURF_ATM")
  ENDIF

  IF (lhook) CALL dr_hook('HOR_EXTRAPOL_SURF_CHEAP OMP',1,zhook_handle_omp)
END DO
!$OMP END PARALLEL DO
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('HOR_EXTRAPOL_SURF_CHEAP',1,zhook_handle)
!
END SUBROUTINE hor_extrapol_surf_cheap