find_variable.f

Go to the documentation of this file.

!                   ************************
                    SUBROUTINE find_variable
!                   ************************
!
     &(fformat,fid,var_name,res,n,ierr,time,eps_time,record,
     & time_record, offset)
!
!***********************************************************************
! HERMES   V7P2
!***********************************************************************
!
!brief    Returns the value for each point of a given variable,
!+        for a given time or a given record.
!+        This subroutine can interpolate from two time steps
!+        if the time does not exist in the file with an accuracy of EPS.
!+
!+        If no optional parameter is given, the last record is chosen.
!+
!+        Accepted formats so far: 'SERAFIN ', 'SERAFIND', 'MED     '
!
!history  Y. AUDOUIN (EDF LAB, LNHE)
!+        28/07/2015
!+        V7P1
!+   First version.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        10/06/2016
!+        V7P2
!+   Just a message corrected.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| EPS_TIME       |-->| EPSILON TO DETERMINE IF TWO TIMES ARE EQUAL
!|                |   | IF NOT PRESENT, CONSIDERED TO BE 0.D0
!| IERR           |<--| 0 IF NO ERROR DURING THE EXECUTION
!| FFORMAT        |-->| FORMAT OF THE FILE
!| FID            |-->| FILE LOGICAL UNIT
!| RECORD         |-->| RECORD OF RESULTS IN THE FILE THAT SHOULD BE
!|                |   | TAKEN. RECORD = -1 TAKES THE LAST RECORD.
!| RES            |<->| VALUE FOR EACH POINT AT TIME STEP RECORD
!|                |   | FOR THE VARIABLE VAR_NAME
!| TIME           |-->| TIME TO INTERPOLATE FROM THE FILE
!| TIME_RECORD    |<--| TIME OF THE RECORD THAT HAS BEEN CHOSEN
!| VAR_NAME       |-->| VARIABLE FOR WHICH WE NEED THE VALUE
!| N              |-->| SIZE OF EXPECTED RES_VALUE
!| OFFSET         |-->| OFFSET FOR TIME
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ONLY : date_mjd2sec
      USE interface_hermes
      USE declarations_special
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      CHARACTER(LEN=8),  INTENT(IN)    :: FFORMAT
      INTEGER,           INTENT(IN)    :: FID
      CHARACTER(LEN=16), INTENT(IN)    :: VAR_NAME
      INTEGER,           INTENT(IN)    :: N
      DOUBLE PRECISION,  INTENT(INOUT) :: RES(n)
      INTEGER,           INTENT(INOUT) :: IERR
      DOUBLE PRECISION,  INTENT(IN)   , OPTIONAL :: EPS_TIME,TIME
      INTEGER,           INTENT(IN)   , OPTIONAL :: RECORD
      DOUBLE PRECISION,  INTENT(INOUT), OPTIONAL :: TIME_RECORD
      DOUBLE PRECISION,  INTENT(IN), OPTIONAL :: OFFSET
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER RECORD1,RECORD2
      DOUBLE PRECISION :: TIME1,TIME2,COEF,EPS
      DOUBLE PRECISION,ALLOCATABLE :: RES2(:)
      INTEGER NTIMESTEP,I
      INTEGER :: RRECORD
      CHARACTER(LEN=16), ALLOCATABLE :: VAR_LIST(:)
      CHARACTER(LEN=16), ALLOCATABLE :: UNIT_LIST(:)
      INTEGER :: NVAR
      DOUBLE PRECISION :: TEL_OFFSET1
      INTEGER, DIMENSION(6) :: DATE
      DOUBLE PRECISION :: OFFSET_FILE
      DOUBLE PRECISION :: TIME_FULL
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      IF(PRESENT(offset)) THEN
        tel_offset1 = offset
      ELSE
        tel_offset1 = 0.d0
      ENDIF

      IF(PRESENT(time)) THEN
!
!       TIME ASKED
!
!       EPSILON TO BE CONSIDERED
        IF(PRESENT(eps_time))THEN
          eps=eps_time
        ELSE
          eps=0.d0
        ENDIF
!
!       GET THE NUMBER OF TIME STEPS
!
        CALL get_data_ntimestep(fformat,fid,ntimestep,ierr)
        CALL check_call(ierr,'FIND_VARIABLE:GET_DATA_TIME:REC1')
!
!       IDENTIFY THE TWO RECORDS IN BETWEEN OUR TIME
!
        record1 = 0
        record2 = ntimestep-1
        CALL get_data_time(fformat,fid,record1,time1,ierr)
        CALL check_call(ierr,'FIND_VARIABLE:GET_DATA_TIME:REC1')
        CALL get_data_time(fformat,fid,record2,time2,ierr)
        CALL check_call(ierr,'FIND_VARIABLE:GET_DATA_TIME:REC2')

        IF (tel_offset1.GE.1.d-16) THEN
          CALL get_mesh_date(fformat, fid, date, ierr)
          IF (any(date.NE.0)) THEN
            offset_file = date_mjd2sec(date(1:3),date(4:6))
            time_full = time + tel_offset1
          ELSE
            time_full = time
            offset_file = 0.d0
          ENDIF
        ELSE
          time_full = time
          offset_file = 0.d0
        ENDIF

        time1 = time1+offset_file
        time2 = time2+offset_file

!
!       Quick check to see if the time is indeed within the first
!       and last time step of the file
!
        IF((time1.GT.time_full.AND.abs(time1-time_full).GT.eps)
     & .OR.(time2.LT.time_full.AND.abs(time2-time_full).GT.eps)) THEN
          WRITE(lu,*) 'TIME: ',time_full,
     &                'IS NOT WITHIN THE RANGE',
     &                ' OF THE FILE I.E.[',
     &                time1,'-',time2,']'
          CALL plante(1)
          stop
        ENDIF
!
        DO WHILE ((time1.LT.time_full).AND.
     &       (abs(time1-time_full).GT.eps) )
          record1 = record1 + 1
          CALL get_data_time(fformat,fid,record1,time1,ierr)
          time1 = time1+offset_file
          CALL check_call(ierr,'FIND_VARIABLE:GET_DATA_TIME:REC1')
        ENDDO
!       If the time is on the file return the value for that record
!       otherwise we interpolate between the two time step
        IF(abs(time1-time_full).LE.eps) THEN
          CALL get_data_value(fformat,fid,record1,var_name,res,n,ierr)
          CALL check_call(ierr,'FIND_VARIABLE:GET_DATA_VALUE:REC1')
        ELSE
!         Get the time and values for the record before and after time
          record2 = record1
          record1 = record1 -1
          CALL get_data_time(fformat,fid,record1,time1,ierr)
          CALL check_call(ierr,'FIND_VARIABLE:GET_DATA_TIME:REC1')
          CALL get_data_time(fformat,fid,record2,time2,ierr)
          time1 = time1+offset_file
          time2 = time2+offset_file
          CALL check_call(ierr,'FIND_VARIABLE:GET_DATA_TIME:REC2')
!         Get the value for each time step
          CALL get_data_value(fformat,fid,record1,var_name,res,n,ierr)
          CALL check_call(ierr,'FIND_VARIABLE:GET_DATA_VALUE:REC1')
          ALLOCATE(res2(n),stat=ierr)
          CALL check_allocate(ierr,'FIND_VARIABLE:RES2')
          CALL get_data_value(fformat,fid,record2,var_name,res2,n,ierr)
          CALL check_call(ierr,'FIND_VARIABLE:GET_DATA_VALUE:REC1')
!
!         Interpolates in time
          coef=(time_full-time1)/(time2-time1)
!
          DO i=1,n
            res(i)=(res2(i)-res(i))*coef+res(i)
          ENDDO
          DEALLOCATE(res2)
        ENDIF
!
      ELSE
        CALL get_data_ntimestep(fformat,fid,ntimestep,ierr)
        CALL check_call(ierr, 'READ_DATA:GET_DATA_NTIMESTEP')
        ! Check if record was given if not 0 is taken
        IF (PRESENT(record)) THEN
          ! Taking last time step
          IF(record.EQ.-1) THEN
            rrecord = ntimestep - 1
          ELSE
            rrecord = record
          ENDIF
        ELSE
          rrecord = 0
        ENDIF
        ! WE CHECK IF THE RECORD IS IN THE FILE BY GETTING THE NUMBER OF
        ! TIMESTEPS AND CHECKING THAT THE RECORD IS INDEED
        ! BETWEEN 0 AND NTIMESTEP - 1
        IF((rrecord.LT.0).AND.(rrecord.GE.ntimestep)) THEN
          ierr = hermes_record_unknown_err
          RETURN
        ENDIF
        ! CHECKING THAT THE VARIABLE IS INDEED IN THE FILE
        !  BY LOOPING ON THE LIST OF VARIABLE
        CALL get_data_nvar(fformat,fid,nvar,ierr)
        CALL check_call(ierr, 'FIND_VARIABLE:GET_DATA_NVAR')
        ALLOCATE(var_list(nvar),stat=ierr)
        ALLOCATE(unit_list(nvar),stat=ierr)
        CALL check_allocate(ierr, 'VAR_LIST')
        CALL get_data_var_list(fformat,fid,nvar,var_list,
     &                         unit_list,ierr)
        CALL check_call(ierr, 'FIND_VARIABLE:GET_DATA_VAR_LIST')
        i=1
        DO WHILE(i.LE.nvar)
          IF (var_list(i).EQ.var_name) EXIT
          i = i +1
        ENDDO
        IF(i.EQ.nvar+1) THEN
          ierr = hermes_var_unknown_err
          RETURN
        ENDIF
        DEALLOCATE(var_list)
        DEALLOCATE(unit_list)

        ! Getting the data value by record number
!
        ! If asked for giving the time value associated with the record
        IF(PRESENT(time_record)) THEN
          CALL get_data_time(fformat,fid,rrecord,time_record,ierr)
          IF(ierr.NE.0) THEN
            WRITE(lu,*) 'ERROR WHILE READING TIME VALUE ',
     &                  'FOR RECORD:',rrecord
            CALL plante(1)
            stop
          ENDIF
        ENDIF
        CALL get_data_value(fformat,fid,rrecord,var_name,res,
     &                      n,ierr)
        IF(ierr.NE.0) THEN
          WRITE(lu,*) 'ERROR WHILE READING VALUE ',
     &                'FOR VARIABLE:',var_name,
     &                'FOR RECORD:',rrecord
          CALL plante(1)
          stop
        ENDIF
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END