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
!!      P. Samuelsson 10/2014 Moved subroutine REGROT to separate routine
!!------------------------------------------------------------------
!
!
!
USE modd_prep,       ONLY : xlat_out, xlon_out
USE modd_grid_rotlatlon
USE modd_surf_par,   ONLY : xundef
USE modd_grid_grib,  ONLY : nni
!
!RJ: missing modi
USE modi_regrot_lonlat_rot
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      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
!

 INTEGER, ALLOCATABLE :: ii(:),jj(:)

 REAL,    ALLOCATABLE :: xlat_ind(:),xlon_ind(:),  &
                           xrat_out(:),xron_out(:),  &
                           w00(:),w01(:),            &
                           w10(:),w11(:)  

 LOGICAL, ALLOCATABLE :: lmask(:)

INTEGER :: i,j,k,l,ij,ij00,ij01,ij10,ij11,ino,jl
REAL    :: wx,wy,wsum
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)

!
!*      1.    Allocations
!
ALLOCATE(xrat_out(ino),       &
           xron_out(ino),       &
           xlat_ind(ino),       &
           xlon_ind(ino),       &
                 ii(ino),       &
                 jj(ino),       &
                w00(ino),       &
                w01(ino),       &
                w10(ino),       &
                w11(ino))  

ALLOCATE(lmask(nni))
!
!*  Transformation of latitudes/longitudes into rotated coordinates

    WRITE(kluout,*)'XLAT_OUT',xlat_out(10:10)
    WRITE(kluout,*)'XLON_OUT',xlon_out(10:10)

    CALL regrot_lonlat_rot(xlon_out,xlat_out,   &
                  xron_out,xrat_out,            &
                  ino,1,ino,1,                  &
                  xrlop,xrlap,1                 )  

    WRITE(kluout,*)'XRAT_OUT',xrat_out(10:10)
    WRITE(kluout,*)'XRON_OUT',xron_out(10:10)

    DO ij=1,ino
       xlat_ind(ij) = ( xrat_out(ij) - xrila1) / xrdy + 1.
       xlon_ind(ij) = ( xron_out(ij) - xrilo1) / xrdx + 1.
    ENDDO

    pfieldout(:,:) = xundef

    DO jl=1,SIZE(pfieldin,2)

    lmask= .true.
    WHERE ( abs(pfieldin(:,jl)-xundef) < 1.e-6 ) lmask = .false.

    DO ij=1,ino

         ii(ij)  = int(xlon_ind(ij))
         jj(ij)  = int(xlat_ind(ij))

         wx  = xlon_ind(ij) - float(ii(ij))
         wy  = xlat_ind(ij) - float(jj(ij))

         w00(ij) = (1.-wx)*(1.-wy)
         w01(ij) = (1.-wx)*    wy
         w10(ij) =     wx *(1.-wy)
         w11(ij) =     wx *    wy

         k    = ii(ij)
         l    = jj(ij)
         ij00 = k   + nrx*(l   -1)
         ij01 = k   + nrx*(l+1 -1)
         ij10 = k+1 + nrx*(l   -1)
         ij11 = k+1 + nrx*(l+1 -1)

         IF (.NOT. lmask(ij00)) w00(ij) = 0.
         IF (.NOT. lmask(ij01)) w01(ij) = 0.
         IF (.NOT. lmask(ij10)) w10(ij) = 0.
         IF (.NOT. lmask(ij11)) w11(ij) = 0.

            wsum = w00(ij) + w01(ij) + &
                     w10(ij) + w11(ij)  

            IF ( abs(wsum) < 1.e-6 ) cycle

            w00(ij) = w00(ij) / wsum
            w01(ij) = w01(ij) / wsum
            w10(ij) = w10(ij) / wsum
            w11(ij) = w11(ij) / wsum

    ENDDO

    !
    ! Bi linear
    !

    WRITE(kluout,*)'NRX,NRY',nrx,nry

       DO ij=1,ino

          k    = ii(ij)
          l    = jj(ij)
          ij00 = k   + nrx*(l   -1)
          ij01 = k   + nrx*(l+1 -1)
          ij10 = k+1 + nrx*(l   -1)
          ij11 = k+1 + nrx*(l+1 -1)

          WRITE(kluout,*)pfieldin(ij00,jl)

          pfieldout(ij,jl) = w00(ij)*pfieldin(ij00,jl) +       &
                               w01(ij)*pfieldin(ij01,jl) +       &
                               w10(ij)*pfieldin(ij10,jl) +       &
                               w11(ij)*pfieldin(ij11,jl)  
 
       ENDDO
    ENDDO


!
!*      5.    Deallocations
!
DEALLOCATE(xrat_out,xron_out,    &
             xlat_ind,xlon_ind,    &
             ii,jj,                &
             w00,w01,w10,w11,      &
             lmask)  
!
IF (lhook) CALL dr_hook('HOR_INTERPOL_ROTLATLON',1,zhook_handle)
!-------------------------------------------------------------------------------------
END SUBROUTINE hor_interpol_rotlatlon