arpege_stretch_a.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 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
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
  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
!-------------------------------------------------------------------------------