mode_geo_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.
MODULE mode_geo_gauss
  !!
  !!     AUTHOR
  !!     ------
  !!      J-D Gril
  !!
  !!     MODIFICATION
  !!     ------------
  !!       Original  10/2005
  !!
  !--------------------------------------------------------------------------- 
  USE eggangles, ONLY : lola
  !
!
  USE yomhook   ,ONLY : lhook,   dr_hook
  USE parkind1  ,ONLY : jprb
!
  IMPLICIT NONE
  !
  INTERFACE gauss_tr
     MODULE PROCEDURE gauss_tr_v, gauss_tr_s
  END INTERFACE
  !
  INTERFACE gauss_rt
     MODULE PROCEDURE gauss_rt_v, gauss_rt_s
  END INTERFACE
  !
  INTERFACE etir
     MODULE PROCEDURE etir_v, etir_s
  END INTERFACE
  !
  INTERFACE retre
     MODULE PROCEDURE retre_v, retre_s
  END INTERFACE
  !
  INTERFACE rotate_g
     MODULE PROCEDURE rotate_g_v, rotate_g_s
  END INTERFACE
  !
  INTERFACE anti_rotate_g
     MODULE PROCEDURE anti_rotate_g_v, anti_rotate_g_s
  END INTERFACE
  ! 
  INTERFACE map_fac
     MODULE PROCEDURE map_fac_v, map_fac_s
  END INTERFACE
  !
 CONTAINS
  !==========================GAUSS_TR=========================================
  TYPE(lola) function gauss_tr_s (pt_reel,pt_pole,pcodil) result (pt_g)
    TYPE (lola), INTENT(IN)          :: pt_reel
    TYPE (lola), INTENT(IN)          :: pt_pole
    REAL, INTENT(IN)                 :: pcodil
    REAL(KIND=JPRB) :: zhook_handle

    IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:GAUSS_TR_S',0,zhook_handle)
    pt_g=etir(anti_rotate_g(pt_reel,pt_pole),pcodil)
  IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:GAUSS_TR_S',1,zhook_handle)

  END FUNCTION gauss_tr_s

  FUNCTION gauss_tr_v (PT_REEL,PT_POLE,PCODIL) RESULT (PT_G)
    TYPE (lola), DIMENSION(:), INTENT(IN)          :: pt_reel
    TYPE (lola), INTENT(IN)          :: pt_pole
    REAL, INTENT(IN)                 :: pcodil
    TYPE (lola), DIMENSION(SIZE(PT_REEL))       :: pt_g
    REAL(KIND=JPRB) :: zhook_handle

    IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:GAUSS_TR_V',0,zhook_handle)
    pt_g=etir(anti_rotate_g(pt_reel,pt_pole),pcodil)
  IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:GAUSS_TR_V',1,zhook_handle)

  END FUNCTION gauss_tr_v
  !==========================GAUSS_RT=========================================
  TYPE(lola) function gauss_rt_s (pt_g,pt_pole,pcodil) result (pt_reel)
    TYPE (lola), INTENT(IN)          :: pt_g
    TYPE (lola), INTENT(IN)          :: pt_pole
    REAL, INTENT(IN)                 :: pcodil
    REAL(KIND=JPRB) :: zhook_handle

    IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:GAUSS_RT_S',0,zhook_handle)
    pt_reel=rotate_g(retre(pt_g,pcodil),pt_pole)
  IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:GAUSS_RT_S',1,zhook_handle)

  END FUNCTION gauss_rt_s

  FUNCTION gauss_rt_v (PT_G,PT_POLE,PCODIL) RESULT (PT_REEL)
    TYPE (lola), DIMENSION(:), INTENT(IN)          :: pt_g
    TYPE (lola), INTENT(IN)          :: pt_pole
    REAL, INTENT(IN)                 :: pcodil
    TYPE (lola), DIMENSION(SIZE(PT_G))     :: pt_reel
    REAL(KIND=JPRB) :: zhook_handle

    IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:GAUSS_RT_V',0,zhook_handle)
    pt_reel=rotate_g(retre(pt_g,pcodil),pt_pole)
  IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:GAUSS_RT_V',1,zhook_handle)

  END FUNCTION gauss_rt_v
  !==========================ROTATE_G=========================================
  TYPE(lola) function rotate_g_s (pt_tr,pt_pole) result (pt_reel)
    USE eggangles, ONLY : cosin_to_angle, p_asin
    TYPE (lola), INTENT(IN)          :: pt_tr
    TYPE (lola), INTENT(IN)          :: pt_pole

    REAL(KIND=JPRB)    :: zsin, zcos
    REAL(KIND=JPRB) :: zhook_handle

    IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:ROTATE_G_S',0,zhook_handle)
    pt_reel%LAT = p_asin(sin(pt_pole%LAT)*sin(pt_tr%LAT)+cos(pt_pole%LAT)*cos(pt_tr%LAT)*cos(pt_tr%LON)) 
    
    IF (cos(pt_reel%LAT) /= 0.0) THEN     
      zsin = -1.0*cos(pt_tr%LAT)*sin(pt_tr%LON)/cos(pt_reel%LAT)
      zcos = (cos(pt_pole%LAT)*sin(pt_tr%LAT)-sin(pt_pole%LAT)*cos(pt_tr%LAT)*cos(pt_tr%LON))/cos(pt_reel%LAT)
      pt_reel%LON = cosin_to_angle(zcos,zsin) + pt_pole%LON
    ELSE
      pt_reel%LON = 0.0
    ENDIF
  IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:ROTATE_G_S',1,zhook_handle)

  END FUNCTION rotate_g_s

  FUNCTION rotate_g_v (PT_TR,PT_POLE) RESULT (PT_REEL)
    USE eggangles, ONLY : cosin_to_angle, p_asin
    TYPE (lola), DIMENSION(:), INTENT(IN)          :: pt_tr
    TYPE (lola), INTENT(IN)                        :: pt_pole
    TYPE (lola), DIMENSION(SIZE(PT_TR))            :: pt_reel

    REAL(KIND=JPRB), DIMENSION(SIZE(PT_TR))                   :: zsin, zcos
    REAL(KIND=JPRB) :: zhook_handle

    IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:ROTATE_G_V',0,zhook_handle)
    pt_reel%LAT = p_asin(sin(pt_pole%LAT)*sin(pt_tr%LAT)+cos(pt_pole%LAT)*cos(pt_tr%LAT)*cos(pt_tr%LON)) 
    
    pt_reel%LON = 0.0
    WHERE (cos(pt_reel%LAT) /= 0.0)
      zsin = -1.0*cos(pt_tr%LAT)*sin(pt_tr%LON)/cos(pt_reel%LAT)
      zcos = (cos(pt_pole%LAT)*sin(pt_tr%LAT)-sin(pt_pole%LAT)*cos(pt_tr%LAT)*cos(pt_tr%LON))/cos(pt_reel%LAT)
      pt_reel%LON = cosin_to_angle(zcos,zsin) + pt_pole%LON
    ENDWHERE
  IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:ROTATE_G_V',1,zhook_handle)


!    DO I=1,SIZE(PT_TR)
!      PT_REEL(I) = ROTATE_G(PT_TR(I),PT_POLE)
!    END DO

  END FUNCTION rotate_g_v
  !==========================ANTI_ROTATE_G====================================
  TYPE(lola) function anti_rotate_g_s (pt_reel,pt_pole) result (pt_tr)
    USE eggangles, ONLY : cosin_to_angle, p_asin
    TYPE (lola), INTENT(IN)          :: pt_reel
    TYPE (lola), INTENT(IN)          :: pt_pole

    REAL(KIND=JPRB) :: zsin, zcos
    REAL(KIND=JPRB) :: zhook_handle
   
    IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:ANTI_ROTATE_G_S',0,zhook_handle)
    pt_tr%LAT = p_asin(sin(pt_pole%LAT)*sin(pt_reel%LAT)+cos(pt_pole%LAT)*cos(pt_reel%LAT)*cos(pt_reel%LON-pt_pole%LON))
 
   IF (cos(pt_tr%LAT) /= 0.0) THEN     
      zsin = -1.0*cos(pt_reel%LAT)*sin(pt_reel%LON-pt_pole%LON)/cos(pt_tr%LAT)
      zcos = (cos(pt_pole%LAT)*sin(pt_reel%LAT)-sin(pt_pole%LAT)*cos(pt_reel%LAT)*cos(pt_reel%LON-pt_pole%LON))/cos(pt_tr%LAT)
      pt_tr%LON = cosin_to_angle(zcos,zsin)
    ELSE
      pt_tr%LON = 0.0
    ENDIF
  IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:ANTI_ROTATE_G_S',1,zhook_handle)
  END FUNCTION anti_rotate_g_s

  FUNCTION anti_rotate_g_v (PT_REEL,PT_POLE) RESULT (PT_TR)
    USE eggangles, ONLY : cosin_to_angle, p_asin
    TYPE (lola), DIMENSION(:), INTENT(IN)          :: pt_reel
    TYPE (lola), INTENT(IN)                        :: pt_pole
    TYPE (lola), DIMENSION(SIZE(PT_REEL))          :: pt_tr  

    REAL(KIND=JPRB), DIMENSION(SIZE(PT_REEL)) :: zsin, zcos
    REAL(KIND=JPRB) :: zhook_handle

   IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:ANTI_ROTATE_G_V',0,zhook_handle)
   pt_tr%LAT = p_asin(sin(pt_pole%LAT)*sin(pt_reel%LAT)+cos(pt_pole%LAT)*cos(pt_reel%LAT)*cos(pt_reel%LON-pt_pole%LON))
 
    pt_tr%LON = 0.0
    WHERE (cos(pt_tr%LAT) /= 0.0)
      zsin = -1.0*cos(pt_reel%LAT)*sin(pt_reel%LON-pt_pole%LON)/cos(pt_tr%LAT)
      zcos = (cos(pt_pole%LAT)*sin(pt_reel%LAT)-sin(pt_pole%LAT)*cos(pt_reel%LAT)*cos(pt_reel%LON-pt_pole%LON))/cos(pt_tr%LAT)
      pt_tr%LON = cosin_to_angle(zcos,zsin)
    ENDWHERE
  IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:ANTI_ROTATE_G_V',1,zhook_handle)

!    DO I=1,SIZE(PT_REEL)
!      PT_TR(I) = ANTI_ROTATE_G(PT_REEL(I),PT_POLE)
!    END DO

  END FUNCTION anti_rotate_g_v
  !==========================ETIR=============================================
  TYPE(lola) function etir_s (pt_coord,pcodil) result (pt_et)
    USE eggangles, ONLY : cosin_to_angle
    TYPE (lola), INTENT(IN)          :: pt_coord
    REAL, INTENT(IN)                 :: pcodil

    REAL(KIND=JPRB) :: zpc2, zsin, zcos
    REAL(KIND=JPRB) :: zhook_handle

    IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:ETIR_S',0,zhook_handle)
    pt_et%LON = pt_coord%LON
    zpc2  = pcodil*pcodil
    zsin = ((1.0-zpc2)+(1.0+zpc2)*sin(pt_coord%LAT))/((1.0+zpc2)+(1.0-zpc2)*sin(pt_coord%LAT))
    zcos = 2.0*pcodil*cos(pt_coord%LAT)/((1.0+zpc2)+(1.0-zpc2)*sin(pt_coord%LAT))

    pt_et%LAT = cosin_to_angle(zcos,zsin)
  IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:ETIR_S',1,zhook_handle)

  END FUNCTION etir_s

  FUNCTION etir_v (PT_COORD,PCODIL) RESULT (PT_ET)
    USE eggangles, ONLY : cosin_to_angle
    TYPE (lola), DIMENSION(:), INTENT(IN)   :: pt_coord
    REAL, INTENT(IN)                        :: pcodil
    TYPE (lola), DIMENSION(SIZE(PT_COORD))  :: pt_et

    REAL(KIND=JPRB), DIMENSION(SIZE(PT_COORD))         :: zpc2, zsin, zcos
    REAL(KIND=JPRB) :: zhook_handle

    IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:ETIR_V',0,zhook_handle)
    pt_et%LON = pt_coord%LON
    zpc2  = pcodil*pcodil
    zsin = ((1.0-zpc2)+(1.0+zpc2)*sin(pt_coord%LAT))/((1.0+zpc2)+(1.0-zpc2)*sin(pt_coord%LAT))
    zcos = 2.0*pcodil*cos(pt_coord%LAT)/((1.0+zpc2)+(1.0-zpc2)*sin(pt_coord%LAT))

    pt_et%LAT = cosin_to_angle(zcos,zsin)
  IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:ETIR_V',1,zhook_handle)

  END FUNCTION etir_v
  !==========================RETRE============================================
  TYPE(lola) function retre_s (pt_coord,pcodil) result (pt_re)
    USE eggangles, ONLY : cosin_to_angle
    TYPE (lola), INTENT(IN)          :: pt_coord
    REAL, INTENT(IN)                 :: pcodil

    REAL(KIND=JPRB)     :: zpc2, zsin, zcos
    REAL(KIND=JPRB) :: zhook_handle

    IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:RETRE_S',0,zhook_handle)
    pt_re%LON = pt_coord%LON
    zpc2  = pcodil*pcodil
    zsin = (-1.0*(1.0-zpc2)+(1.0+zpc2)*sin(pt_coord%LAT))/((1.0+zpc2)-(1.0-zpc2)*sin(pt_coord%LAT))
    zcos = 2.0*pcodil*cos(pt_coord%LAT)/((1.0+zpc2)-(1.0-zpc2)*sin(pt_coord%LAT))

    pt_re%LAT = cosin_to_angle(zcos,zsin)
  IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:RETRE_S',1,zhook_handle)

  END FUNCTION retre_s

  FUNCTION retre_v (PT_COORD,PCODIL) RESULT (PT_RE)
    USE eggangles, ONLY : cosin_to_angle
    TYPE (lola), DIMENSION(:), INTENT(IN)   :: pt_coord
    REAL, INTENT(IN)                        :: pcodil
    TYPE (lola), DIMENSION(SIZE(PT_COORD))  :: pt_re

    REAL(KIND=JPRB), DIMENSION(SIZE(PT_COORD))         :: zpc2, zsin, zcos
    REAL(KIND=JPRB) :: zhook_handle

    IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:RETRE_V',0,zhook_handle)
    pt_re%LON = pt_coord%LON
    zpc2  = pcodil*pcodil
    zsin = (-1.0*(1.0-zpc2)+(1.0+zpc2)*sin(pt_coord%LAT))/((1.0+zpc2)-(1.0-zpc2)*sin(pt_coord%LAT))
    zcos = 2.0*pcodil*cos(pt_coord%LAT)/((1.0+zpc2)-(1.0-zpc2)*sin(pt_coord%LAT))

    pt_re%LAT = cosin_to_angle(zcos,zsin)
  IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:RETRE_V',1,zhook_handle)

!    DO I=1,SIZE(PT_COORD)
!       PT_RE(I) = RETRE(PT_COORD(I),PCODIL)
!    END DO

  END FUNCTION retre_v
  !==========================MAP FACTOR=======================================
  REAL FUNCTION map_fac_s (POLE,PCODIL,PT_COORD) RESULT (PMF)
    TYPE (lola), INTENT(IN)          :: pole, pt_coord
    REAL, INTENT(IN)                 :: pcodil
    
    REAL(KIND=JPRB) :: zpc2
    REAL(KIND=JPRB) :: zhook_handle
    IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:MAP_FAC_S',0,zhook_handle)
    zpc2  = pcodil*pcodil
    pmf = (2.0*pcodil)/((1+zpc2)+(1-zpc2)* &
            (sin(pt_coord%LAT)*sin(pole%LAT)+cos(pt_coord%LAT)*cos(pole%LAT)*cos(pt_coord%LON-pole%LON)))      
  IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:MAP_FAC_S',1,zhook_handle)
   
  END FUNCTION map_fac_s

  FUNCTION map_fac_v (POLE,PCODIL,PT_COORD) RESULT (PMF)
    TYPE (lola), INTENT(IN)                 :: pole
    TYPE (lola), DIMENSION(:), INTENT(IN)   :: pt_coord
    REAL, INTENT(IN)                        :: pcodil
    REAL, DIMENSION(SIZE(PT_COORD))         :: pmf

    REAL(KIND=JPRB), DIMENSION(SIZE(PT_COORD))         :: zpc2
    REAL(KIND=JPRB) :: zhook_handle

    IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:MAP_FAC_V',0,zhook_handle)
    zpc2  = pcodil*pcodil
    pmf = (2.0*pcodil)/((1+zpc2)+(1-zpc2)* &
            (sin(pt_coord%LAT)*sin(pole%LAT)+cos(pt_coord%LAT)*cos(pole%LAT)*cos(pt_coord%LON-pole%LON)))      
  IF (lhook) CALL dr_hook('MODE_GEO_GAUSS:MAP_FAC_V',1,zhook_handle)

!    DO I=1,SIZE(PT_COORD)
!      PMF(I) = MAP_FAC(POLE,PCODIL,PT_COORD(I))
!    END DO    
    
  END FUNCTION map_fac_v
END MODULE mode_geo_gauss