! $Id: atm2geo.F90 5144 2024-07-29 21:01:04Z abarral $

SUBROUTINE atm2geo(im, jm, pte, ptn, plon, plat, pxx, pyy, pzz)
  USE dimphy
  USE lmdz_phys_para
  USE lmdz_grid_phy, ONLY: grid_type, unstructured, regular_lonlat
  USE lmdz_abort_physic, ONLY: abort_physic
  USE lmdz_yomcst

  IMPLICIT NONE

  CHARACTER (len = 6) :: clmodnam
  CHARACTER (len = 20) :: modname = 'atm2geo'
  CHARACTER (len = 80) :: abort_message

  ! Change wind local atmospheric coordinates to geocentric

  ! Geocentric :
  ! axe x is eastward : crosses (0 N, 0 E) point.
  ! axe y crosses (0 N, 90 E) point.
  ! axe z is 'up' : crosses north pole
  INTEGER, INTENT (IN) :: im, jm
  REAL, DIMENSION (im, jm), INTENT (IN) :: pte  ! Eastward vector component
  REAL, DIMENSION (im, jm), INTENT (IN) :: ptn  ! Northward vector component
  REAL, DIMENSION (im, jm), INTENT (IN) :: plon, plat
  REAL, DIMENSION (im, jm), INTENT(OUT) :: pxx, pyy, pzz  ! Component in the geocentric referential
  REAL :: rad, reps

  rad = rpi / 180.0E0
  reps = 1.0e-5

  pxx(:, :) = &
          - pte(:, :) * SIN(rad * plon(:, :)) &
                  - ptn(:, :) * SIN(rad * plat(:, :)) * COS(rad * plon(:, :))

  pyy(:, :) = &
          + pte(:, :) * COS(rad * plon(:, :)) &
                  - ptn(:, :) * SIN(rad * plat(:, :)) * SIN(rad * plon(:, :))

  pzz(:, :) = &
          + ptn(:, :) * COS(rad * plat (:, :))

  IF (grid_type==regular_lonlat) THEN
    ! Value at North Pole
    IF (is_north_pole_dyn) THEN
      pxx(:, 1) = - pte (1, 1)
      pyy(:, 1) = - ptn (1, 1)
      pzz(:, 1) = pzz(1, 1) ! => 0
    ENDIF

    ! Value at South Pole
    IF (is_south_pole_dyn) THEN
      pxx(:, jm) = pxx(1, jm)
      pyy(:, jm) = pyy(1, jm)
      pzz(:, jm) = pzz(1, jm) ! => 0
    ENDIF

  ELSE IF (grid_type==unstructured) THEN
    ! Pole nord pour Dynamico
    WHERE (plat(:, :) >= 90.0d+0 - reps)
      pxx (:, :) = -ptn (:, :)
      pyy (:, :) = pte (:, :)
      pzz (:, :) = 0.0e0
    END WHERE

  ELSE
    abort_message = 'Problem: unknown grid type'
    CALL abort_physic(modname, abort_message, 1)
  END IF

END SUBROUTINE atm2geo