read_mesh_coord.f

Go to the documentation of this file.

!                   **************************
                    SUBROUTINE read_mesh_coord
!                   **************************
!
     &(fformat,nfic,x,y,npoin,projection,lati0,longi0,z)
!
!***********************************************************************
! HERMES   V7P1
!***********************************************************************
!
!brief    Reads the coordinates in the geometry file.
!+        Latitude-longitude coordinates transformed into mercator.
!
!history  Y. AUDOUIN (EDF LAB, LNHE)
!+        06/05/2015
!+        V7P1
!+   First version.
!
!history  M.S.TURNBULL (HRW)
!+        18/09/2015
!+        V7P1
!+   Correction to the conversion from degree to Mercator projection.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| FFORMAT        |-->| FORMAT FOR GEOMETRY FILE
!| NFIC           |-->| LOGICAL UNIT FOR GEOMETRY FILE
!| X              |<--| X COORDINATES
!| Y              |<--| Y COORDINATES
!| PROJECTION     |<--| TYPE OF PROJECTION
!| LATI0          |<--| LATITUDE
!| LONGI0         |<--| LONGITUDE
!| Z              |<--| Z COORDINATES
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_hermes
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)           :: NPOIN,NFIC
      DOUBLE PRECISION, INTENT(OUT) :: X(npoin),Y(npoin)
      DOUBLE PRECISION, INTENT(OUT), OPTIONAL :: Z(npoin)
      INTEGER, INTENT(INOUT)        :: PROJECTION
      DOUBLE PRECISION, INTENT(IN)  :: LATI0,LONGI0
      CHARACTER(LEN=8), INTENT(IN)  :: FFORMAT
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER :: IERR, I, NDIM
!
      DOUBLE PRECISION, PARAMETER :: R=6.37d6
      DOUBLE PRECISION PS4,TAN1,TAN2,LONGIRAD
!
      INTRINSIC log,tan,atan
!
!-----------------------------------------------------------------------
!
!     total number of dimension
      ndim = 2
      IF (PRESENT(z)) ndim = 3

      CALL get_mesh_coord(fformat,nfic,1,ndim,npoin,x,ierr)
      IF(ierr.NE.0) THEN
        WRITE(lu,*) 'ERROR WHILE READING X ARRAY'
        CALL plante(1)
      ENDIF
      CALL get_mesh_coord(fformat,nfic,2,ndim,npoin,y,ierr)
      IF(ierr.NE.0) THEN
        WRITE(lu,*) 'ERROR WHILE READING Y ARRAY'
        CALL plante(1)
      ENDIF
      IF(PRESENT(z)) THEN
        CALL get_mesh_coord(fformat,nfic,3,ndim,npoin,z,ierr)
        IF(ierr.NE.0) THEN
          WRITE(lu,*) 'ERROR WHILE READING Z ARRAY'
          CALL plante(1)
        ENDIF
      ENDIF
!
      CALL corrxy(x,y,npoin)
!
!     LATITUDE-LONGITUDE COORDINATES (DEGREES)
!     CHANGED INTO MERCATOR PROJECTION
!
      IF(projection.EQ.3) THEN
        ps4=atan(1.d0)
        longirad=longi0*ps4/45.d0
        tan2=tan(lati0*ps4/90.d0+ps4)
        IF(tan2.LT.0.d0) THEN
          WRITE(lu,*) 'LATI0=',lati0,' IS PROBABLY WRONG'
          CALL plante(1)
          stop
        ENDIF
        DO i=1,npoin
          x(i)=r*(x(i)*ps4/45.d0-longirad)
          tan1=tan(y(i)*ps4/90.d0+ps4)
          IF(tan1.LT.0.d0) THEN
            WRITE(lu,*) 'LATITUDE MUST BE GIVEN IN DEGREES'
            WRITE(lu,*) 'HERE Y(I)=',y(i),' FOR I=',i
            WRITE(lu,*) 'USE CORRXY (BIEF) FOR CONVERSION'
            CALL plante(1)
            stop
          ENDIF
          y(i)=r*(log(tan1)-log(tan2))
        ENDDO
!       NOW IT IS MERCATOR
        projection=2
        WRITE(lu,*) ' '
        WRITE(lu,*) 'READ_MESH_COORD :'
        WRITE(lu,*) 'COORDINATES CHANGED INTO MERCATOR PROJECTION'
        WRITE(lu,*) 'WITH LATITUDE OF ORIGIN POINT = ',lati0
        WRITE(lu,*) 'AND LONGITUDE OF ORIGIN POINT = ',longi0
        WRITE(lu,*) ' '
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END