hor_interpol_gauss.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_gauss(KLUOUT,PFIELDIN,PFIELDOUT)
!     #################################################################################
!
!!****  *HOR_INTERPOL_GAUSS* - Interpolation from a gaussian grid
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!     V. Masson 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!      M. Jidane   Dec 2013 : initialize NNI if not already done
!!------------------------------------------------------------------
!
USE modd_prep,       ONLY : xlat_out, xlon_out, linterp
USE modd_grid_gauss, ONLY : xila1, xilo1, xila2, xilo2, ninla, ninlo, nilen, lrotpole, &
                              xlap, xlop, xcoef  
USE modd_grid_grib,  ONLY : nni
USE modd_surf_par,   ONLY : xundef
!
USE modi_horibl_surf
!
!
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
!
REAL, DIMENSION(:), ALLOCATABLE :: zlat     ! latitudes
REAL, DIMENSION(:), ALLOCATABLE :: zlon     ! longitudes
INTEGER, DIMENSION(:), ALLOCATABLE :: imaskin  ! input mask
INTEGER, DIMENSION(:), ALLOCATABLE :: imaskout ! output mask
INTEGER                         :: ino      ! output number of points
INTEGER                         :: jl       ! loop counter
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------------
!
!*      1.    Allocations
!
IF (lhook) CALL dr_hook('HOR_INTERPOL_GAUSS',0,zhook_handle)
ino = SIZE(xlat_out)
!
IF (nni==0) nni=nilen
ALLOCATE(imaskin(nni))
!
ALLOCATE(zlat(ino))
ALLOCATE(zlon(ino))
ALLOCATE(imaskout(ino))
imaskout = 1
!
!*      2.     Is pole rotated?
!
IF (lrotpole) THEN
!* transformation of output latitudes, longitudes into rotated coordinates
  CALL arpege_stretch_a(ino,xlap,xlop,xcoef,xlat_out,xlon_out,zlat,zlon)
ELSE
  zlat = xlat_out 
  zlon = xlon_out 
END IF
!
!*      3.    Input mask
!
DO jl=1,SIZE(pfieldin,2)
!
imaskin(:) = 1.
WHERE(pfieldin(:,jl)==xundef) imaskin = 0.
!
!
!*      4.    Interpolation with horibl
!
  CALL horibl_surf(xila1,xilo1,xila2,xilo2,ninla,ninlo,nni,pfieldin(:,jl),ino,zlon,zlat,pfieldout(:,jl), &
                .false.,kluout,linterp,imaskin,imaskout)  
END DO
!
!
!*      5.    Deallocations
!
DEALLOCATE(zlat)
DEALLOCATE(zlon)
DEALLOCATE(imaskin )
DEALLOCATE(imaskout)
!
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('HOR_INTERPOL_GAUSS',1,zhook_handle)
 CONTAINS
!
!     ##########################################################################
SUBROUTINE arpege_stretch_a(KN,PLAP,PLOP,PCOEF,PLAR,PLOR,PLAC,PLOC)
!     ##########################################################################
!!****  *ARPEGE_STRETCH_A* - Projection to Arpege stretched grid
!!
!!   PURPOSE
!!   -------
!!
!!   Projection from standard Lat,Lon grid to Arpege stretched grid
!!
!!   METHOD
!!   ------
!!
!!   The projection is defined in two steps :
!!    1. A rotation to place the stretching pole at the north pole
!!    2. The stretching
!!   This routine is a basic implementation of the informations founded in 
!!     'Note de travail Arpege n#3'
!!     'Transformation de coordonnees'
!!     J.F.Geleyn 1988
!!   This document describes a slightly different transformation in 3 steps. Only the
!!   two first steps are to be taken in account (at the time of writing this paper has
!!   not been updated).
!!
!!   EXTERNAL
!!   --------
!!
!!
!!   IMPLICIT ARGUMENTS
!!   ------------------
!!
!!   REFERENCE
!!   ---------
!!
!!   This routine is based on : 
!!     'Note de travail ARPEGE' number 3
!!     by J.F. GELEYN (may 1988)
!!
!!   AUTHOR
!!   ------
!!
!!   V.Bousquet
!!
!!   MODIFICATIONS
!!   -------------
!!
!!   Original       07/01/1999
!!   V. Masson      01/2004    Externalization of surface
!!
!
! 0. DECLARATIONS
! ---------------
!
  USE modd_csts,ONLY : xpi
!
  IMPLICIT NONE
!
! 0.1. Declaration of arguments
! -----------------------------

  INTEGER,             INTENT(IN)  :: kn            ! Number of points to convert
  REAL,                INTENT(IN)  :: plap          ! Latitude of stretching pole
  REAL,                INTENT(IN)  :: plop          ! Longitude of stretching pole
  REAL,                INTENT(IN)  :: pcoef         ! Stretching coefficient
  REAL, DIMENSION(KN), INTENT(IN)  :: plar          ! Lat. of points
  REAL, DIMENSION(KN), INTENT(IN)  :: plor          ! Lon. of points
  REAL, DIMENSION(KN), INTENT(OUT) :: plac          ! Computed pseudo-lat. of points
  REAL, DIMENSION(KN), INTENT(OUT) :: ploc          ! Computed pseudo-lon. of points
!
  REAL                             :: zsinstretchla ! Sine of stretching point lat.
  REAL                             :: zsinstretchlo ! Sine of stretching point lon.
  REAL                             :: zcosstretchla ! Cosine of stretching point lat.
  REAL                             :: zcosstretchlo ! Cosine of stretching point lon.
  REAL                             :: zsinla        ! Sine of computed point latitude
  REAL                             :: zsinlo        ! Sine of computed point longitude
  REAL                             :: zcosla        ! Cosine of computed point latitude
  REAL                             :: zcoslo        ! Cosine of computed point longitude
  REAL                             :: zsinlas       ! Sine of point's pseudo-latitude
  REAL                             :: zsinlos       ! Sine of point's pseudo-longitude
  REAL                             :: zcoslos       ! Cosine of point's pseudo-lon.
  REAL                             :: za,zb,zd      ! Dummy variables used for 
  REAL                             :: zx,zy         ! computations
!  
  INTEGER                          :: jloop1        ! Dummy loop counter
  REAL(KIND=JPRB) :: zhook_handle
!
  IF (lhook) CALL dr_hook('ARPEGE_STRETCH_A',0,zhook_handle)
  zsinstretchla = sin(plap*xpi/180.)
  zcosstretchla = cos(plap*xpi/180.)
  zsinstretchlo = sin(plop*xpi/180.)
  zcosstretchlo = cos(plop*xpi/180.)
  ! L = longitude (0 = Greenwich, + toward east)
  ! l = latitude (90 = N.P., -90 = S.P.)
  ! p stands for stretching pole
  plac(:) = plar(:) * xpi / 180.
  ploc(:) = plor(:) * xpi / 180.
  ! A = 1 + c.c
  za = 1. + pcoef*pcoef
  ! B = 1 - c.c
  zb = 1. - pcoef*pcoef
  DO jloop1=1, kn
    zsinla = sin(plac(jloop1))
    zcosla = cos(plac(jloop1))
    zsinlo = sin(ploc(jloop1))
    zcoslo = cos(ploc(jloop1))
    ! X = cos(Lp-L)
    zx = zcoslo*zcosstretchlo + zsinlo*zsinstretchlo
    ! Y = sin(Lp-L)
    zy = zsinstretchlo*zcoslo - zsinlo*zcosstretchlo
    ! D = (1+c.c) + (1-c.c)(sin lp.sin l + cos lp.cos l.cos(Lp-L))
    zd = za + zb*(zsinstretchla*zsinla+zcosstretchla*zcosla*zx)
    !          (1-c.c)+(1+c.c)((sin lp.sin l + cos lp.cos l.cos(Lp-L))
    ! sin lr = -------------------------------------------------------
    !                                  D
    zsinlas = (zb + za*(zsinstretchla*zsinla+zcosstretchla*zcosla*zx)) / zd
    ! D' = D * cos lr
    zd = zd * (amax1(1e-6,sqrt(1.-zsinlas*zsinlas)))
    !          2.c.(cos lp.sin l - sin lp.cos l.cos(Lp-L))
    ! cos Lr = -------------------------------------------
    !                              D'
    zcoslos = 2.*pcoef*(zcosstretchla*zsinla-zsinstretchla*zcosla*zx) / zd
    !          2.c.cos l.cos(Lp-L)
    ! sin Lr = -------------------
    !                  D'
    zsinlos = 2.*pcoef*(zcosla*zy) / zd
    ! saturations (corrects calculation errors)
    zsinlas = max(zsinlas,-1.)
    zsinlas = min(zsinlas, 1.)
    zcoslos = max(zcoslos,-1.)
    zcoslos = min(zcoslos, 1.)
    ! back from sine & cosine
    plac(jloop1) = asin(zsinlas)
    IF (zsinlos>0) THEN
      ploc(jloop1) =  acos(zcoslos)
    ELSE
      ploc(jloop1) = -acos(zcoslos)
    ENDIF
  ENDDO
  ploc(:) = ploc(:) * 180. / xpi
  plac(:) = plac(:) * 180. / xpi
IF (lhook) CALL dr_hook('ARPEGE_STRETCH_A',1,zhook_handle)
END SUBROUTINE arpege_stretch_a
!-------------------------------------------------------------------------------
END SUBROUTINE hor_interpol_gauss