longitu.f

Go to the documentation of this file.

!                   ******************
                    SUBROUTINE longitu
!                   ******************
!
     &(xel,coslat,ikle,nelmax,nelem)
!
!***********************************************************************
! BIEF   V7P1
!***********************************************************************
!
!brief    CORRECT THE ARRAYS THAT DEPEND ON THE LONGITUDE
!+                OF THE GIVEN ELEMENT.
!
!history  M.S.TURNBULL (HRW) and S.E.BOURBAN (HRW)
!+        01/10/2015
!+        V7P1
!+   Correction to the computation of the XEL in spherical coordinates.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| XEL            |<->| VECTOR DEFINED PER ELEMENT
!| COSLAT         |-->| COSINUS OF LATITUDE OF ORIGIN POINT
!| IKLE           |-->| CONNECTIVITY TABLE
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NELEM,NELMAX
      INTEGER         , INTENT(IN)    :: IKLE(nelmax,3)
      DOUBLE PRECISION, INTENT(IN)    :: COSLAT(*)
      DOUBLE PRECISION, INTENT(INOUT) :: XEL(nelmax,3)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELEM
!
      DOUBLE PRECISION R,PIR,SCAEL
!
      INTRINSIC atan
!
!-----------------------------------------------------------------------
!
! EARTH RADIUS
!
      r = 6370000.d0
!
!-----------------------------------------------------------------------
!
!     FILTER XEL ON THE DATELINE, AND THEN ONLY SCALE
!
!     WARNING: XEL ON EITHER SIDE OF THE SAME SEGMENT MIGHT BE DIFFERENT
!     ANOTHER SOLUTION WHERE SCAEL is SCAEDGE MIGHT BE MORE APPROPRIATE
!
      pir = 4.d0 * atan(1.d0) * r
!
      DO ielem = 1,nelem
!
        IF(xel(ielem,2).GT.pir) THEN
          xel(ielem,2) = xel(ielem,2) - 2.d0 * pir
        ELSEIF(xel(ielem,2).LT.-pir) THEN
          xel(ielem,2) = xel(ielem,2) + 2.d0 * pir
        ENDIF
!
        IF(xel(ielem,3).GT.pir) THEN
          xel(ielem,3) = xel(ielem,3) - 2.d0 * pir
        ELSEIF(xel(ielem,3).LT.-pir) THEN
          xel(ielem,3) = xel(ielem,3) + 2.d0 * pir
        ENDIF
!
        scael = ( coslat(ikle(ielem,1))
     &          + coslat(ikle(ielem,2)) +
     &            coslat(ikle(ielem,3)) ) / 3.d0
        xel(ielem,2) = xel(ielem,2) * scael
        xel(ielem,3) = xel(ielem,3) * scael
!
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END