source: LMDZ4/branches/V3_test/libf/phylmd/atm2geo.F @ 3919

!
! $Header$
!
C
      SUBROUTINE atm2geo ( im, jm, pte, ptn, plon, plat, pxx, pyy, pzz )
      USE dimphy
      IMPLICIT NONE
      include 'dimensions90.h'
cc
cc Change wind local atmospheric coordinates to
cc geocentric
cc
cxxx      INCLUDE 'param.h'
c
      INTEGER, INTENT (in)              :: im, jm
      REAL, DIMENSION (im,jm), INTENT (in) :: pte, ptn
      REAL, DIMENSION (im,jm), INTENT (in) :: plon, plat
      REAL, DIMENSION (im,jm), INTENT(out) :: pxx, pyy, pzz
c
      REAL, PARAMETER :: rpi = 3.141592653E0
      REAL, PARAMETER :: rad = rpi / 180.0E0
c
      REAL, DIMENSION (im,jm) :: zsinlon, zcoslon
      REAL, DIMENSION (im,jm) :: zsinlat, zcoslat
c
      LOGICAL, SAVE :: linit = .FALSE.
c$OMP THREADPRIVATE(linit)
c
cxxx      IF ( .NOT. linit ) THEN 
          zsinlon = SIN (rad * plon)
          zcoslon = COS (rad * plon)
          zsinlat = SIN (rad * plat)
          zcoslat = COS (rad * plat)
          linit = .TRUE.
cxxx      ENDIF 
c
      pxx = - zsinlon * pte - zsinlat * zcoslon * ptn
      pyy =   zcoslon * pte - zsinlat * zsinlon * ptn
      pzz =   zcoslat * ptn

c
c Value at North Pole
      IF (jjphy_begin==1) THEN
        pxx ( :,  1) = - ptn ( 1, 1)
        pyy ( :,  1) = - pte ( 1, 1)
        pzz ( :,  1) = 0.0
      ENDIF
c Value at South Pole
      
      IF (jjphy_end==jjm+1) THEN
        pxx ( :, jm) = + ptn ( 1, jm)
        pyy ( :, jm) = + pte ( 1, jm)
        pzz ( :, jm) = 0.0
      ENDIF
 
      RETURN 
      END SUBROUTINE atm2geo