hor_interpol_rotlatlon.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_interpol_rotlatlon(KLUOUT,PFIELDIN,PFIELDOUT)
!     #################################################################################
!
!!****  *HOR_INTERPOL_ROTLATLON * - Interpolation from a rotated lat/lon grid
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     U. Andrae
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    10/2007
!!------------------------------------------------------------------
!
!
USE modd_surfex_mpi, ONLY : nrank, npio, ncomm, nproc, idx_i
USE modd_prep,       ONLY : xlat_out, xlon_out
USE modd_grid_rotlatlon
USE modd_surf_par,   ONLY : xundef
USE modd_grid_grib,  ONLY : nni

!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
#ifdef SFX_MPI
include "mpif.h"
#endif
!
!*      0.1    declarations of arguments
!
INTEGER,              INTENT(IN)  :: KLUOUT    ! logical unit of output listing
REAL, DIMENSION(:,:), INTENT(IN)  :: PFIELDIN  ! field to interpolate horizontally
REAL, DIMENSION(:,:), INTENT(OUT) :: PFIELDOUT ! interpolated field
!
!*      0.2    declarations of local variables
!
#ifdef SFX_MPI
INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ISTATUS
#endif
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZW, ZFIELDIN,ZFGET
INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: IJD,IJGET
REAL, DIMENSION(:), ALLOCATABLE :: XLAT_IND, XLON_IND, XRAT_OUT, XRON_OUT
!
INTEGER :: INFOMPI
INTEGER :: J, I, JI, INO, JL, INL, ILON, ILAT, ISIZE
REAL    :: ZWX, ZWY, ZWSUM
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('HOR_INTERPOL_ROTLATLON',0,zhook_handle)
!
WRITE(kluout,'(A)')' | Running rotated latlon interpolation'
!
ino = SIZE(xlat_out)
inl = SIZE(pfieldout,2)
!
!
!*      1.    Allocations
!
ALLOCATE(xrat_out(ino),xron_out(ino),xlat_ind(ino),xlon_ind(ino),ijd(ino,inl,4),zw(ino,inl,4))  
!
!*  Transformation of latitudes/longitudes into rotated coordinates
!
 CALL regrot(xlon_out,xlat_out,xron_out,xrat_out,xrlop,xrlap,1)  
!
DO j=1,ino
  xlat_ind(j) = ( xrat_out(j) - xrila1) / xrdy + 1.
  xlon_ind(j) = ( xron_out(j) - xrilo1) / xrdx + 1.
ENDDO
!
pfieldout(:,:) = xundef
!
zw(:,:,:) = 0.
DO jl=1,inl
  !
  DO ji=1,ino
    !
    ilon = int(xlon_ind(ji))
    ilat = int(xlat_ind(ji))
    !
    zwx = xlon_ind(ji) - float(ilon)
    zwy = xlat_ind(ji) - float(ilat)
    !
    zw(ji,jl,1) = (1.-zwx)*(1.-zwy)
    zw(ji,jl,2) = (1.-zwx)*    zwy
    zw(ji,jl,3) =     zwx *(1.-zwy)
    zw(ji,jl,4) =     zwx *    zwy
    !
    ijd(ji,jl,1) = ilon   + nrx*(ilat   -1)
    ijd(ji,jl,2) = ilon   + nrx*(ilat+1 -1)
    ijd(ji,jl,3) = ilon+1 + nrx*(ilat   -1)
    ijd(ji,jl,4) = ilon+1 + nrx*(ilat+1 -1)    
    !
  ENDDO
  !
ENDDO
!
ALLOCATE(zfieldin(ino,inl,4))
!
IF (nrank/=npio) THEN
  !
  idx_i = idx_i + 1
#ifdef SFX_MPI  
  CALL mpi_send(ino,kind(ino)/4,mpi_integer,npio,idx_i,ncomm,infompi)
#endif
  !
  idx_i = idx_i + 1
#ifdef SFX_MPI  
  CALL mpi_send(ijd,SIZE(ijd)*kind(ijd)/4,mpi_integer,npio,idx_i,ncomm,infompi)
#endif  
  !
  idx_i = idx_i + 1
#ifdef SFX_MPI  
  CALL mpi_recv(zfieldin,SIZE(zfieldin)*kind(zfieldin)/4,mpi_real,npio,idx_i,ncomm,istatus,infompi)
#endif  
  !
ELSE
  !
  DO j=0,nproc-1
    !
    IF (j/=npio) THEN
#ifdef SFX_MPI            
      CALL mpi_recv(isize,kind(isize)/4,mpi_integer,j,idx_i+1,ncomm,istatus,infompi)
#endif      
      ALLOCATE(ijget(isize,inl,4))
#ifdef SFX_MPI      
      CALL mpi_recv(ijget,SIZE(ijget)*kind(ijget)/4,mpi_integer,j,idx_i+2,ncomm,istatus,infompi)
#endif      
    ELSE
      isize = ino
      ALLOCATE(ijget(isize,inl,4))
      ijget(:,:,:) = ijd(:,:,:)
    ENDIF
    !
    ALLOCATE(zfget(isize,inl,4))
    !
    zfget(:,:,:) = 0.
    DO jl = 1,inl
      DO ji = 1,isize
        DO i = 1,4
          zfget(ji,jl,i) = pfieldin(ijget(ji,jl,i),jl)
        ENDDO
      ENDDO
    ENDDO
    !
    IF (j/=npio) THEN
#ifdef SFX_MPI            
      CALL mpi_send(zfget,SIZE(zfget)*kind(zfget)/4,mpi_real,j,idx_i+3,ncomm,infompi)
#endif      
    ELSE
      zfieldin(:,:,:) = zfget(:,:,:)
    ENDIF
    !
    DEALLOCATE(ijget,zfget)
    !
  ENDDO
  !
  idx_i = idx_i + 3
  !
ENDIF 
!
DO jl = 1,inl
  DO ji = 1,ino
    !
    DO i = 1,4
      IF (abs(zfieldin(ji,jl,i)-xundef)<1.e-6) zw(ji,jl,i) = 0.
    ENDDO
    !
    zwsum = zw(ji,jl,1) + zw(ji,jl,2) + zw(ji,jl,3) + zw(ji,jl,4)
    !
    IF ( abs(zwsum)<1.e-6 ) THEN
      zw(ji,jl,1) = 1.
    ELSE
      DO i = 1,4
        zw(ji,jl,i) = zw(ji,jl,i)/zwsum
      ENDDO
    ENDIF
    !
    pfieldout(ji,jl) = zw(ji,jl,1)*zfieldin(ji,jl,1) + zw(ji,jl,2)*zfieldin(ji,jl,2) + &
                       zw(ji,jl,3)*zfieldin(ji,jl,3) + zw(ji,jl,4)*zfieldin(ji,jl,4)
    !
  ENDDO
ENDDO
!
!*      5.    Deallocations
DEALLOCATE(xrat_out,xron_out,xlat_ind,xlon_ind,ijd,zfieldin)  
!
IF (lhook) CALL dr_hook('HOR_INTERPOL_ROTLATLON',1,zhook_handle)
!
CONTAINS
!
SUBROUTINE regrot(PXREG,PYREG,PXROT,PYROT,PXCEN,PYCEN,KCALL)  
!
USE modd_csts, ONLY : xpi
!
USE modi_abor1_sfx
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!*    CONVERSION BETWEEN REGULAR AND ROTATED SPHERICAL COORDINATES.
!*
!*    PXREG     LONGITUDES OF THE REGULAR COORDINATES
!*    PYREG     LATITUDES OF THE REGULAR COORDINATES
!*    PXROT     LONGITUDES OF THE ROTATED COORDINATES
!*    PYROT     LATITUDES OF THE ROTATED COORDINATES
!*              ALL COORDINATES GIVEN IN DEGREES N (NEGATIVE FOR S)
!*              AND DEGREES E (NEGATIVE VALUES FOR W)
!*    KXDIM     DIMENSION OF THE GRIDPOINT FIELDS IN THE X-DIRECTION
!*    KYDIM     DIMENSION OF THE GRIDPOINT FIELDS IN THE Y-DIRECTION
!*    KX        NUMBER OF GRIDPOINT IN THE X-DIRECTION
!*    KY        NUMBER OF GRIDPOINTS IN THE Y-DIRECTION
!*    PXCEN     REGULAR LONGITUDE OF THE SOUTH POLE OF THE ROTATED GRID
!*    PYCEN     REGULAR LATITUDE OF THE SOUTH POLE OF THE ROTATED GRID
!*
!*    KCALL=-1: FIND REGULAR AS FUNCTIONS OF ROTATED COORDINATES.
!*    KCALL= 1: FIND ROTATED AS FUNCTIONS OF REGULAR COORDINATES.
!*
!*    J.E. HAUGEN   HIRLAM   JUNE -92
!
!-----------------------------------------------------------------------
!
INTEGER, INTENT(IN) :: KCALL
REAL, INTENT(IN) :: PXCEN,PYCEN 
REAL, DIMENSION(:), INTENT(INOUT) :: PXREG, PYREG
REAL, DIMENSION(:), INTENT(INOUT) :: PXROT, PYROT         
!
!-----------------------------------------------------------------------
!
REAL :: ZRAD,ZSYCEN,ZCYCEN,ZXMXC,ZSXMXC,ZCXMXC,ZSYREG,ZCYREG, &
        ZSYROT,ZCYROT,ZCXROT,ZSXROT,ZRADI  
INTEGER :: JI, ISIZE
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-----------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('REGROT',0,zhook_handle)
!
isize = SIZE(pxreg)
!
zrad = xpi/180.
zradi = 1./zrad
zsycen = sin(zrad*(pycen+90.))
zcycen = cos(zrad*(pycen+90.))
!
IF (kcall.EQ.1) THEN
  !
  DO ji = 1,isize
    !
    zxmxc  = zrad*(pxreg(ji) - pxcen)
    zsxmxc = sin(zxmxc)
    zcxmxc = cos(zxmxc)
    !
    zsyreg = sin(zrad*pyreg(ji))
    zcyreg = cos(zrad*pyreg(ji))
    !
    zsyrot = zcycen*zsyreg - zsycen*zcyreg*zcxmxc
    zsyrot = min(max(zsyrot,-1.0),+1.0)
    !
    pyrot(ji) = asin(zsyrot)*zradi
    !
    zcyrot = cos(zrad*pyrot(ji))
    zcxrot = (zcycen*zcyreg*zcxmxc + zsycen*zsyreg)/zcyrot  
    zcxrot = min(max(zcxrot,-1.0),+1.0)
    zsxrot = zcyreg*zsxmxc/zcyrot
    !
    pxrot(ji) = acos(zcxrot)*zradi
    !
    IF (zsxrot.LT.0.0) pxrot(ji) = -pxrot(ji)
    !
  ENDDO
  !
ELSEIF (kcall.EQ.-1) THEN
  !
  DO ji = 1,isize
    !
    zsxrot = sin(zrad*pxrot(ji))
    zcxrot = cos(zrad*pxrot(ji))
    zsyrot = sin(zrad*pyrot(ji))
    zcyrot = cos(zrad*pyrot(ji))
    !
    zsyreg = zcycen*zsyrot + zsycen*zcyrot*zcxrot
    zsyreg = max(zsyreg,-1.0)
    zsyreg = min(zsyreg,+1.0)
    !
    pyreg(ji) = asin(zsyreg)*zradi
    !
    zcyreg = cos(pyreg(ji)*zrad)
    !
    zcxmxc = (zcycen*zcyrot*zcxrot - zsycen*zsyrot)/zcyreg  
    zcxmxc = max(zcxmxc,-1.0)
    zcxmxc = min(zcxmxc,+1.0)
    zsxmxc = zcyrot*zsxrot/zcyreg
    zxmxc  = acos(zcxmxc)*zradi
    IF (zsxmxc.LT.0.0) zxmxc = -zxmxc
    !
    pxreg(ji) = zxmxc + pxcen
    !
  ENDDO
  !
ELSE
  !
  WRITE(6,'(1X,''INVALID KCALL IN REGROT'')')
  CALL abor1_sfx('HOR_INTERPOL_ROTLATON:REGROT:KCALL MUST BE 1 OR -1')
  !
ENDIF
!   
IF (lhook) CALL dr_hook('REGROT',1,zhook_handle)
!
END SUBROUTINE regrot
!
!-------------------------------------------------------------------------------------
END SUBROUTINE hor_interpol_rotlatlon