C Copyright 1981-2007 ECMWF
C 
C Licensed under the GNU Lesser General Public License which
C incorporates the terms and conditions of version 3 of the GNU
C General Public License.
C See LICENSE and gpl-3.0.txt for details.
C

      INTEGER FUNCTION INTUVPH( KFIELD1, KFIELD2, INLEN,
     X                          KUGRIBO, KVGRIBO,
     X                          OUTLENU, OUTLENV)
C
C---->
C**** INTUVPH
C
C     Purpose
C     -------
C
C     Interpolate GRIB format input spectral vorticity and divergence,
C     or spectral U and V, fields to GRIB format U and V fields.
C
C
C     Interface
C     ---------
C
C     IRET = INTUVPH(KFIELD1,KFIELD2,INLEN,
C                    KUGRIBO,KVGRIBO,OUTLENU,OUTLENV)
C
C     Input
C     -----
C
C     KFIELD1 - Input vorticity or U field  (spectral, GRIB format).
C     KFIELD2 - Input divergence or V field (spectral, GRIB format).
C     INLEN  - Input field length (words).
C
C
C     Output
C     ------
C
C     KUGRIBO - Output U field (GRIB format).
C     KVGRIBO - Output V field (GRIB format).
C     OUTLENU - Output U field length (words).
C     OUTLENV - Output V field length (words).
C
C
C     Method
C     ------
C
C     Convert spectral vorticity/divergence to spectral U/V and then
C     interpolate U and V to output fields.
C
C     Note that a common block is used in intf.h to hold the U/V
C     fields before interpolation.
C
C     Externals
C     ---------
C
C     IBASINI - Ensure basic interpolation setup is done.
C     INTUVDH - Encodes/decodes data into/from GRIB code.
C     INTUVXH - Interpolate U or V component spectral field to grid point.
C     JVOD2UV - Converts spectral vorticity/divergence to spectral U/V.
C     JMEMHAN - Allocate scratch memory.
C     GRIBEX  - GRIB decoding/encoding.
C     ISCRSZ  - Calculate number of values in generated field.
C     FIXAREA - Fixup area definition to correspond to grid definitions
C     INTLOG  - Log error message.
C     RESET_C - Reset interpolation handling options using GRIB product.
C     INSANE  - Ensure no outrageous values given for interpolation.
C     GRSMKP  - P factor calculation switch for routine GRIBEX.
C
C
C     Author
C     ------
C
C     J.D.Chambers     ECMWF     February 2001
C
C
C----<
C     -----------------------------------------------------------------|
C
      IMPLICIT NONE
C
C     Function arguments
C
      INTEGER KFIELD1(*), KFIELD2(*), INLEN
      INTEGER KUGRIBO(*), KVGRIBO(*), OUTLENU, OUTLENV
C
#include "parim.h"
#include "nifld.common"
#include "nofld.common"
#include "grfixed.h"
#include "intf.h"
C
C     Parameters
C
      INTEGER JPROUTINE, JPALLOC, JPDEALL, JPSCR3, JPSCR4,JPSCR5
      PARAMETER (JPROUTINE = 40160 )
      PARAMETER (JPALLOC = 1) 
      PARAMETER (JPDEALL = 0) 
      PARAMETER (JPSCR3 = 3) 
      PARAMETER (JPSCR4 = 4) 
      PARAMETER (JPSCR5 = 5) 
C
C     Local variables
C
      CHARACTER*1 HFUNC
      REAL EW, NS
      LOGICAL LOLDWIND
      INTEGER IERR,KPR,ISZVD,ISZUV,IWORD,ISIZE,ILENF,ISAME,IPARAM
      INTEGER NEXT, LOOP, MTRUNC, NTRUNC, NTROLD, NTROLD2, NPARAM
      INTEGER MIRESO, MORESO, NOLD, NLEN
C
      LOGICAL LFIRST, LNEWUV
      CHARACTER*3 EXTRA
      DATA LFIRST/.TRUE./, LNEWUV/.TRUE./, EXTRA/'NO '/
      SAVE LFIRST, LNEWUV
C
      DATA NTROLD/-1/, NTROLD2/-1/
      SAVE NTROLD, NTROLD2
      INTEGER IPVORT, IPDIV, IP_U, IP_V, IDIVOFF
#ifdef POINTER_64
      INTEGER*8 IZNFLDO
#endif
      REAL ZNFLDO
      POINTER ( IZNFLDO, ZNFLDO )
      DIMENSION ZNFLDO( 1 )
#ifdef POINTER_64
      INTEGER*8 IUV, IVD
#endif
      REAL UV, VD
      POINTER ( IUV, UV )
      POINTER ( IVD, VD )
      DIMENSION UV( 1 ), VD( 1 )
C
C     Externals
C
      INTEGER RESET_C,ISCRSZ,FIXAREA,AURESOL
      INTEGER IBASINI, INSANE, INTUVDH, INTUVXH
C
C     -----------------------------------------------------------------|
C*    Section 1.   Initialise
C     -----------------------------------------------------------------|
C
  100 CONTINUE
C
      INTUVPH = 0
      IERR    = 0
      KPR     = 0
C
      IF( LFIRST ) THEN
        CALL GETENV('IGNORE_UV_EXTRA_MODE', EXTRA)
        IF((EXTRA(1:1).EQ.'Y').OR.(EXTRA(1:1).EQ.'y')) LNEWUV = .FALSE.
        IF( LNEWUV ) THEN
          CALL INTLOG(JP_DEBUG,
     X      'INTUVPH: IGNORE_UV_EXTRA_MODE not turned on',JPQUIET)
        ELSE
          CALL INTLOG(JP_DEBUG,
     X      'INTUVPH: IGNORE_UV_EXTRA_MODE turned on',JPQUIET)
        ENDIF
        LFIRST = .FALSE.
      ENDIF
C
      MIRESO = NIRESO
      MORESO = NORESO
C
      LOLDWIND = LWINDSET
C
C     Ensure that basic initialisation has been done
C
      IERR = IBASINI(0)
      IF( IERR.NE.0 ) THEN
        CALL INTLOG(JP_ERROR,'INTUVPH: basic initialise failed',JPQUIET)
        INTUVPH = IERR
        GOTO 900
      ENDIF
C
C     -----------------------------------------------------------------|
C*    Section 2.   Unpack the input fields.
C     -----------------------------------------------------------------|
C
  200 CONTINUE
C
C     Need to establish input truncation, so unpack GRIB sections 1
C     and 2.
C
      IPARAM = 0
      ISZVD = 1
      IERR = INTUVDH(VD(IPVORT),ISZVD,KFIELD1,INLEN,'I',IPARAM)
      IF( IERR.NE.0 ) THEN
        CALL INTLOG(JP_ERROR,'INTUVPH: GRIB header decode failed',IERR)
        INTUVPH = IERR
        GOTO 900
      ENDIF
C
C     Check that the input is an ECMWF spectral field
C
      IF( (ISEC1(1).NE.128).OR.(ISEC2(1).NE.50) ) THEN
        CALL INTLOG(JP_ERROR,
     X    'INTUVPH: Input is not ECMWF spectral field.',JPQUIET)
        INTUVPH = JPROUTINE + 2
        GOTO 900
      ENDIF
C
C     Get scratch memory for input unpacked fields.
C     Unpacked field memory areas are adjacent.
C
      NIRESO = ISEC2(2)
      ISZVD  = (NIRESO+1)*(NIRESO+2)
      IPVORT = 1
      IPDIV  = 1 + ISZVD
      CALL JMEMHAN( JPSCR4, IVD, ISZVD*2, JPALLOC, IERR)
      IF( IERR.NE.0 ) THEN
        CALL INTLOG(JP_ERROR,
     X    'INTUVPH: Scratch memory type 4 allocation failed.',JPQUIET)
        INTUVPH = IERR
        GOTO 900
      ENDIF
C
C     Decode input fields..
C
      IERR = INTUVDH(VD(IPVORT),ISZVD,KFIELD1,INLEN,'D',IPARAM)
      IF( IERR.NE.0 ) THEN
        CALL INTLOG(JP_ERROR,'INTUVPH: Vorticity decoding failed',IERR)
        INTUVPH = IERR
        GOTO 900
      ENDIF
C
      IERR = INTUVDH(VD(IPDIV),ISZVD,KFIELD2,INLEN,'D',IPARAM)
      IF( IERR.NE.0 ) THEN
        CALL INTLOG(JP_ERROR,'INTUVPH: Divergence decoding failed',IERR)
        INTUVPH = IERR
        GOTO 900
      ENDIF
C
C     Setup interpolation options from input GRIB characteristics.
C
      IERR = RESET_C( ISEC1, ISEC2, ZSEC2, ISEC4)
      IF( IERR.NE.0 ) THEN
        CALL INTLOG(JP_ERROR,
     X    'INTUVPH: Setup interp. options from GRIB failed.',JPQUIET)
        INTUVPH = IERR
        GOTO 900
      ENDIF
C
C     Check that no outrageous values given for interpolation
C
      ISAME = INSANE()
      IF( (ISAME.GT.0).AND.(ISAME.NE.27261) ) THEN
        CALL INTLOG(JP_ERROR,
     X    'INTUVPH: Interpolation cannot use given values.',JPQUIET)
        INTUVPH = ISAME
        GOTO 900
      ENDIF
C
C     -----------------------------------------------------------------|
C*    Section 3.   Convert spectral vorticity/divergence
C                  to spectral U/V
C     -----------------------------------------------------------------|
C
  300 CONTINUE
C
C     Spectral U and V for Tn are to be generated from vorticity
C     and divergence spectral T(n-1)
C
C     Get scratch memory for U and V spectral fields.
C     The memory areas are adjacent.
C
      IF( ISEC1(6).EQ.155 ) THEN
C
C       Truncate vorticity and divergence to correspond to U/V
C
        NTRUNC = NORESO - 1
C
        IF( LNEWUV ) THEN
          MTRUNC = NTRUNC + 1
        ELSE
          NTRUNC = NTRUNC + 1
          MTRUNC = NTRUNC
        ENDIF
C
C       Check whether user requested a truncation before interpolation
C
        IF( .NOT.LNORESO ) THEN
C
C         If user did not supply a truncation value, see whether or
C         not the 'autoresol' flag is set. If not, use the input
C         truncation.
C
          IF( LARESOL ) THEN
            IF( (NOREPR.EQ.JPREGULAR).OR.(NOREPR.EQ.JPREGROT) ) THEN
              EW = FLOAT(NOGRID(1))/PPMULT
              NS = FLOAT(NOGRID(2))/PPMULT
            ELSE
              EW = 90.0/FLOAT(NOGAUSS)
              NS = EW
            ENDIF
            NTRUNC = AURESOL(EW,NS)
            IF( NTRUNC.NE.NTROLD ) THEN
              NTROLD = NTRUNC
              CALL INTLOG(JP_WARN,
     X          'INTUVPH: Resolution automatically set to ', NTRUNC)
            ENDIF
          ELSE
            NTRUNC = NIRESO
          ENDIF
C
        ENDIF
C
C       Check whether the output resolution is greater than the input
C
        IF( NTRUNC.GT.NIRESO ) THEN
C
C         Issue warning if the output resolution was user-supplied
C
          IF( .NOT.LARESOL ) THEN
C
C           Revert to the input truncation
C
            IF( NIRESO.NE.NTROLD2 ) THEN
              CALL INTLOG(JP_WARN,
     X          'INTUVPH: spectral -> grid point interpolation',JPQUIET)
              CALL INTLOG(JP_WARN,
     X          'INTUVPH: User supplied resolution = ',NTRUNC)
              CALL INTLOG(JP_WARN,
     X          'INTUVPH: Input field resolution   = ',NIRESO)
              CALL INTLOG(JP_WARN,
     X          'INTUVPH: User supplied resolution ignored',JPQUIET)
              CALL INTLOG(JP_WARN,
     X          'INTUVPH: Input field resolution has been used',JPQUIET)
              NTROLD2 = NIRESO
            ENDIF
            NTRUNC = NIRESO
C
          ELSE
C
C           Revert to the input truncation
C
            NTRUNC = NIRESO
            IF( NTRUNC.NE.NTROLD2 ) THEN
              NTROLD2 = NTRUNC
              CALL INTLOG(JP_WARN,
     X          'INTUVPH: Auto-resolution selection too high',JPQUIET)
              CALL INTLOG(JP_WARN,
     X          'INTUVPH: Resolution set to input resolution: ',NTRUNC)
            ENDIF
          ENDIF
        ENDIF
C
        MTRUNC = NTRUNC
        IF( LNEWUV ) NTRUNC = MTRUNC - 1
C
        CALL INTLOG(JP_DEBUG,'INTUVPH: vo/div truncation = ', NTRUNC)
        CALL INTLOG(JP_DEBUG,'INTUVPH: U/V truncation    = ', MTRUNC)
C
        ISIZE =  (NTRUNC+1)*(NTRUNC+2)
        CALL JMEMHAN( JPSCR5, IZNFLDO, ISIZE*2, JPALLOC, IERR)
        IF( IERR.NE.0 ) THEN
          CALL INTLOG(JP_ERROR,
     X      'INTUVP: Scratch memory type 5 allocation failed.',JPQUIET)
          INTUVPH = JPROUTINE + 4
          GOTO 900
        ENDIF
C
        CALL SH2SH( VD(IPVORT), NIRESO, ZNFLDO, NTRUNC )
C
        IDIVOFF = IPVORT + (NTRUNC+1)*(NTRUNC+2)
        CALL SH2SH( VD(IPDIV), NIRESO, ZNFLDO(IDIVOFF), NTRUNC )
C
        NEXT = 0
        DO LOOP = NTRUNC, 0, -1
          NEXT = NEXT + LOOP + 1
          ZNFLDO(IDIVOFF+NEXT*2-2) = 0
          ZNFLDO(IDIVOFF+NEXT*2-1) = 0
          ZNFLDO(IPVORT +NEXT*2-2) = 0
          ZNFLDO(IPVORT +NEXT*2-1) = 0
        ENDDO
C
        ISZUV = (MTRUNC+1)*(MTRUNC+2)
        IP_U  = 1
        IP_V  = 1 + ISZUV
        CALL JMEMHAN( JPSCR3, IUV, ISZUV*2, JPALLOC, IERR)
        IF( IERR.NE.0 ) THEN
          CALL INTLOG(JP_ERROR,
     X      'INTUVPH: Scratch memory type 3 allocation failed.',JPQUIET)
            INTUVPH = IERR
          GOTO 900
        ENDIF
C
C       Generate U and V from vorticity and divergence,
C
        CALL JVOD2UV(ZNFLDO(IPVORT),ZNFLDO(IDIVOFF),NTRUNC,
     X               UV(IP_U),UV(IP_V),MTRUNC)
C
      ELSE
C
C       Get scratch memory for U and V spectral fields.
C       The memory areas are adjacent.
C
        NTRUNC = NIRESO
        ISZUV = (NIRESO+1)*(NIRESO+2)
        IP_U  = 1
        IP_V  = 1 + ISZUV
        CALL JMEMHAN( JPSCR3, IUV, ISZUV*2, JPALLOC, IERR)
        IF( IERR.NE.0 ) THEN
          CALL INTLOG(JP_ERROR,
     X      'INTUVPH: Scratch memory type 3 allocation failed.',JPQUIET)
          INTUVPH = IERR
          GOTO 900
        ENDIF
C
C       Transfer input spectral U and V to the memory areas.
C
        DO LOOP = 0, ISZUV-1
          UV(IP_U + LOOP) = VD(IPVORT + LOOP)
          UV(IP_V + LOOP) = VD(IPDIV  + LOOP)
        ENDDO
C
      ENDIF
C
C     -----------------------------------------------------------------|
C*    Section 4. Handle spectral output.
C     -----------------------------------------------------------------|
C
  400 CONTINUE
C
      IF( (LNORESO)            .AND.
     X    ((NOREPR.EQ.JPSPHERE).OR.(NOREPR.EQ.JPSPHROT)) ) THEN
C
        CALL INTLOG(JP_DEBUG,
     X    'INTUVPH: Produce spectral output with truncation',NORESO)
C
C       Set GRIBEX flag to force recalculation of complex packing factor
C
        CALL GRSMKP(1)
C
        ISIZE =  (NORESO+1)*(NORESO+2)
        CALL JMEMHAN( JPSCR5, IZNFLDO, ISIZE*2, JPALLOC, IERR)
        IF( IERR.NE.0 ) THEN
          CALL INTLOG(JP_ERROR,
     X      'INTUVPH: Scratch memory type 5 allocation failed.',JPQUIET)
          INTUVPH = JPROUTINE + 4
          GOTO 900
        ENDIF
C
        IP_U = 1
        CALL SH2SH( UV(1), NIRESO, ZNFLDO(IP_U), NORESO )
C
        IP_V = 1 + (NORESO+1)*(NORESO+2)
        CALL SH2SH( UV(1+ISZUV), NIRESO, ZNFLDO(IP_V), NORESO )
C
        NIRESO = NORESO
C
      ENDIF
C
C     Has all processing been done (ie is the output spectral)?
C
      IF( (NOREPR.EQ.JPSPHERE).OR.(NOREPR.EQ.JPSPHROT) ) THEN
C
C       Code U into GRIB
C
        IERR = INTUVDH(ZNFLDO(IP_U),ISZUV,KUGRIBO,OUTLENU,'C',JP_U)
        IF( IERR.NE.0 ) THEN
          CALL INTLOG(JP_ERROR,
     X      'INTUVPH: U encoding into GRIB failed.',IERR)
          INTUVPH = JPROUTINE + 4
          GOTO 900
        ENDIF
C
C       Code V into GRIB
C
        IERR = INTUVDH(ZNFLDO(IP_V),ISZUV,KVGRIBO,OUTLENV,'C',JP_V)
        IF( IERR.NE.0 ) THEN
          CALL INTLOG(JP_ERROR,
     X      'INTUVPH: V encoding into GRIB failed.',IERR)
          INTUVPH = JPROUTINE + 4
          GOTO 490
        ENDIF
C
  490   CONTINUE
C
C       Turn off GRIBEX flag which forces recalculation of complex
C       packing factor
C
        CALL GRSMKP(0)
C
        GOTO 900
C
      ENDIF
C
C     -----------------------------------------------------------------|
C*    Section 6.   Generate interpolated GRIB format U and V fields.
C     -----------------------------------------------------------------|
C
  600 CONTINUE
C
C     Get scratch space for interpolation
C
      IERR = FIXAREA()
      IF( IERR.NE.0 ) THEN
        CALL INTLOG(JP_ERROR,'INTUVPH: Fixarea failed.',JPQUIET)
        INTUVPH = IERR
        GOTO 900
      ENDIF
C
      ISIZE = ISCRSZ()*2
      IF( ISIZE.LE.0 ) THEN
        CALL INTLOG(JP_ERROR,'INTUVPH: Get scratch space failed.',ISIZE)
        INTUVPH = JPROUTINE + 5
        GOTO 900
      ENDIF
      CALL JMEMHAN( JPSCR5, IZNFLDO, ISIZE, JPALLOC, IERR)
      IF( IERR.NE.0 ) THEN
        CALL INTLOG(JP_ERROR,
     X    'INTUVPH: Scratch memory(5) allocation failed.',JPQUIET)
        INTUVPH = JPROUTINE + 5
        GOTO 900
      ENDIF
C
      LWIND = .TRUE.
      LWINDSET = .TRUE.
C
C     Interpolate U and V
C
      NOLD = NIRESO
      NIRESO = MTRUNC
      IERR = INTUVXH(UV,ISZUV,ZNFLDO,KUGRIBO,KVGRIBO,OUTLENU,OUTLENV)
      IF( IERR.NE.0 ) THEN
        CALL INTLOG(JP_ERROR,
     X    'INTUVPH: U/V interpolation failed.',JPQUIET)
        INTUVPH = JPROUTINE + 6
        GOTO 900
      ENDIF
C
      NIRESO = NOLD
C
C     -----------------------------------------------------------------|
C*    Section 9.   Return
C     -----------------------------------------------------------------|
C
  900 CONTINUE
C
C     Clear change flags for next product processing and reset wind flag
C
      LCHANGE = .FALSE.
      LSMCHNG = .FALSE.
      LWINDSET  = LOLDWIND
      LWIND = .FALSE.
C
      NIRESO = MIRESO
      NORESO = MORESO
C
      RETURN
      END
