#include "w3macros.h"
!/ ------------------------------------------------------------------- /
      PROGRAM W3PRNC
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           M. Accensi              |
!/                  |           F. Ardhuin              |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         14-Oct-2015 |
!/                  +-----------------------------------+
!/
!/    01-Jan-2011 : Creation                            ( version 4.01 )   
!/    17-Nov-2011 : Fix bug on latitudes                ( version 4.04 )
!/    30-Sep-2012 : Implement tidal analysis            ( version 4.08 )
!/    29-Oct-2012 : Parallelization of tidal analysis   ( version 4.08 )
!/     4-Mar-2012 : allows any NetCDF dimensions names  ( version 4.09 ) 
!/    13-Mar-2012 : Makes compatible with NC3           ( version 4.10 ) 
!/    18-Oct-2013 : Debug compile issue with TIDE switch( version 4.12 )
!/    18-Oct-2013 : Initialize interpolation weights    ( version 4.12 )
!/    20-Dec-2013 : Allow scale factor and offset in    ( version 4.16 )
!/                  NetCDF variables (S. Zieger)
!/    24-Oct-2014 : Allows "As Is" curvilinear grids    ( version 5.02 )
!/    14-Oct-2015 : Add a check for latitude reversed   ( version 5.11 )
!/    20-Jan-2017 : Update to new W3GSRUMD APIs         ( version 6.02 )
!/
!/    Copyright 2009 National Weather Service (NWS),
!/       National Oceanic and Atmospheric Administration.  All rights
!/       reserved.  WAVEWATCH III is a trademark of the NWS. 
!/       No unauthorized use without permission.
!/
!  1. Purpose :
!
!     Pre-processing of the input water level, current, wind and ice
!     fields as well as assimilation data ... from NetCDF input       
!
!  2. Method :
!
!     See documented input file.
!
!  3. Parameters :
!
!     Local parameters.
!     ----------------------------------------------------------------
!       NDSI    Int.  Input unit number ("ww3_prnc.inp").
!       NDSLL   Int.  Unit number(s) of long-lat file(s)
!       NDSF    I.A.  Unit number(s) of input file(s).
!       NDSDAT  Int.  Unit number for output data file.
!       IFLD    Int.  Integer input type.
!       ITYPE   Int.  Integer input 'format' type.
!       NFCOMP  Int.  Number of partial input to be processed.
!       FLTIME  Log.  Time flag for input fields, if false, single
!                     field, time read from NDSI.
!       IDLALL  Int.  Layout indicator used by INA2R. +
!       IDFMLL  Int.  Id. FORMAT indicator.           |
!       FORMLL  C*16  Id. FORMAT.                     | Long-lat
!       FROMLL  C*4   'UNIT' / 'NAME' indicator       |    file(s)
!       NAMELL  C*20  Name of long-lat file(s)        +
!       IDLAF   I.A.   +
!       IDFMF   I.A.   |
!       FORMF   C.A.   | Idem. fields file(s)
!       FROMF   C*4    |
!       NAMEF   C*20   +
!       FORMT   C.A.  Format or time in field.
!       XC      R.A.  Components of input vector field or first
!                     input scalar field
!       XCFAC   Real  Scale factor for input scalar field
!       XCOFF   Real  Offset for input scalar field
!       YC      R.A.  Components of input vector field or second
!                     input scalar field
!       YCFAC   Real  Scale factor for input scalar field
!       YCOFF   Real  Offset for input scalar field
!       FX,FY   R.A.  Output fields.
!       ACC     Real  Required interpolation accuracy.
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!      W3NMOD    Subr. W3GDATMD Set number of model.
!      W3SETG    Subr.   Id.    Point to selected model.
!      W3NDAT    Subr. W3WDATMD Set number of model for wave data.
!      W3SETW    Subr.   Id.    Point to selected model for wave data.
!      W3NOUT    Subr. W3ODATMD Set number of model for output.
!      W3SETO    Subr.   Id.    Point to selected model for output.
!      ITRACE    Subr. W3SERVMD Subroutine tracing initialization.
!      STRACE    Subr.   Id.    Subroutine tracing.
!      NEXTLN    Subr.   Id.    Get next line from input filw
!      EXTCDE    Subr.   Id.    Abort program as graceful as possible.
!      STME21    Subr. W3TIMEMD Convert time to string.
!      INAR2R    Subr. W3ARRYMD Read in an REAL array.
!      INAR2I    Subr.   Id.    Read in an INTEGER array.
!      PRTBLK    Subr.   Id.    Print plot of array.
!      W3IOGR    Subr. W3IOGRMD Reading/writing model definition file.
!      W3FLDO    Subr. W3FLDSMD Opening of WAVEWATCH III generic shell
!                               data file.
!      W3FLDP    Subr.   Id.    Prepare interp. from arbitrary grid.
!      W3FLDG    Subr.   Id.    Reading/writing shell input data.
!      W3FLDD    Subr.   Id.    Reading/writing shell assim. data.
!      W3GSUC    Func. W3GSRUMD Create grid-search-utility object
!      W3GSUD    Subr. W3GSRUMD Destroy grid-search-utility object
!      W3GRMP    Func. W3GSRUMD Compute interpolation weights
!     ----------------------------------------------------------------
!
!  5. Called by :
!
!     None, stand-alone program.
!
!  6. Error messages :
!
!     - Checks on files and reading from file.
!     - Checks on validity of input parameters.
!
!  7. Remarks :
!
!     - Input fields need to be continuous in longitude and latitude.
!     - Program attempts to detect closure type using longitudes of the
!       grid. Thus, it does not allow the user to specify the closure 
!       type, and so tripole closure is not supported.
!
!  8. Structure :
!
!     ----------------------------------------------------
!        1.a  Number of models.
!                   ( W3NMOD , W3NOUT , W3SETG , W3SETO )
!          b  I-O setup.
!          c  Print heading(s).
!        2.   Read model definition file.      ( W3IOGR )
!        3.a  Read major types from input file.
!          b  Check major types.
!          c  Additional input format types and time.
!        4.   Prepare interpolation.
!          a  Longitude - latitude grid
!          b  Grid(s) from file.               ( W3FLDP )
!          c  Initialize fields.
!          d  Input location and format.
!        5    Prepare input and output files.
!          a  Open input file
!          b  Open and prepare output file     ( W3FLDO )
!        6    Until end of file
!          a  Read new time and fields
!          b  Interpolate fields
!          c  Write fields                     ( W3FLDG )
!     ----------------------------------------------------
!
!  9. Switches :
!
!     !/WNT0  = !/WNT1
!     !/WNT1  Correct wind speeds to (approximately) conserve the wind
!             speed over the interpolation box.
!     !/WNT2  Id. energy (USE ONLY ONE !)
!     !/CRT1  Like !/WNT1 for currents.
!     !/CRT2  Like !/WNT2 for currents.
!     !/MPI   Parallel processing is used for tidal analysis.
!
!     !/O3    Additional output in fields processing loop.
!     !/O15   Generate file with the times of the processed fields.
!
!     !/S     Enable subroutine tracing.
!     !/T     Enable test output,
!     !/T1    Full interpolation data.
!     !/T1a   Echo of lat-long data in type Fn
!     !/T2    Full input data.
!     !/T3    Print-plot of output data.
!
!     !/NCO   NCEP NCO modifications for operational implementation.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE CONSTANTS
!/
!     USE W3GDATMD, ONLY: W3NMOD, W3SETG
!/NL1      USE W3ADATMD,ONLY: W3NAUX, W3SETA
      USE W3ODATMD, ONLY: W3NOUT, W3SETO
      USE W3ODATMD, ONLY: IAPROC, NAPROC, NAPERR, NAPOUT
      USE W3SERVMD, ONLY : ITRACE, NEXTLN, EXTCDE, STRSPLIT
!/S      USE W3SERVMD, ONLY : STRACE
      USE W3ARRYMD, ONLY : INA2R, INA2I
!/T2      USE W3ARRYMD, ONLY : PRTBLK
!/T3      USE W3ARRYMD, ONLY : PRTBLK
      USE W3IOGRMD, ONLY: W3IOGR
      USE W3FLDSMD, ONLY: W3FLDO, W3FLDP, W3FLDG, W3FLDD,  &
                          W3FLDTIDE1, W3FLDTIDE2
!/
      USE W3GDATMD
      USE W3GSRUMD
      USE W3ODATMD, ONLY: NDSE, NDST, NDSO, FNMPRE

!/TIDE      USE W3TIDEMD
      USE W3TIMEMD
      USE NETCDF
!
      IMPLICIT NONE
!
!/MPI      INCLUDE "mpif.h"
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
!/MPI      INTEGER                 :: IERR_MPI
      INTEGER                 :: TIDESIZ
      REAL, ALLOCATABLE       :: TIDALSTUFFA(:)
      INTEGER                 :: TIDALSTUFF
      INTEGER, ALLOCATABLE    :: NXL(:)
      INTEGER                 :: NTI, NDSEN, IND
      CHARACTER(LEN=1024)     :: STRFIELDSNAME
      CHARACTER(LEN=100)      :: FIELDSNAME(4)
      CHARACTER(LEN=1024)     :: STRDIMSNAME
      CHARACTER(LEN=100)      :: DIMSNAME(2)
      INTEGER                 :: NIDIMS, NFIELDS
      INTEGER                 :: JPC
      INTEGER(KIND=4)         :: IX, IY, JX, NXLM
      INTEGER                 :: NDSI, NDSM, NDSDAT, NDSTRC, NTRACE,  &
                                 IERR, IFLD, ITYPE, J, NFCOMP,        &
                                 TIME(2), TIMESTART(2),               &
                                 NXI, NYI, NXJ(2), NYJ(2),            &
                                 NDSLL, IDLALL, IDFMLL, NDSF(2),      &
                                 IDLAF(2), IDFMF(2), TIME2(2),        &
                                 MXM, MYM, DATTYP, RECLDT, IDAT,      &
                                 NDAT, NDATMX, JJ, IS(4), JS(4)
      INTEGER                 :: NYB
      INTEGER                 :: ILAND = -999
!/O15      INTEGER                 :: NDSTIME
      INTEGER, ALLOCATABLE    :: IX21(:,:), IX22(:,:),                &
                                 IY21(:,:), IY22(:,:),                &
                                 JX21(:,:), JX22(:,:),                &
                                 JY21(:,:), JY22(:,:), MAPOVR(:,:)
      INTEGER, ALLOCATABLE    :: MASK(:,:)
      TYPE(T_GSU)             :: GSI
!/S      INTEGER, SAVE           :: IENT = 0
!/T2      INTEGER                 :: IXP0, IXPN, IXPWDT = 60
!/T3      INTEGER                 :: IX0, IXN, IXWDT = 60
!/T3      INTEGER, ALLOCATABLE    :: MAPOUT(:,:)
      REAL                    :: X0I, XNI, Y0I, YNI, SXI, SYI, XDIF,  &
                                 ADD, X, Y, FACTOR, EFAC, NODATA, RW(4)
      REAL                    :: ACC = 0.05
      REAL, ALLOCATABLE       :: RD11(:,:), RD21(:,:),                &
                                 RD12(:,:), RD22(:,:),                &
                                 XD11(:,:), XD21(:,:),                &
                                 XD12(:,:), XD22(:,:),                &
                                 FX(:,:), FY(:,:), FA(:,:),           &
                                 A1(:,:), A2(:,:), A3(:,:)
      REAL, POINTER           :: ALA(:,:), ALO(:,:)
      REAL, ALLOCATABLE       :: XC(:,:), YC(:,:), AC(:,:), DATA(:,:)
      REAL                    :: XCFAC, XCOFF, YCFAC, YCOFF
      LOGICAL                 :: INGRID
      LOGICAL                 :: FLSTAB, FLBERG, CLO(2), FLTIME, FLHDR
      INTEGER                 :: ICLO
!/T      LOGICAL                 :: FLMOD
      CHARACTER               :: COMSTR*1, IDFLD*3, IDTYPE*2,         &
                                 IDTIME*23, FROMLL*4, FORMLL*16,      &
                                 NAMELL*80, FROMF*4, NAMEF*80
      CHARACTER*14            :: IDSTR1(-7:5)
      CHARACTER*15            :: IDSTR3(3)
      CHARACTER*32            :: FORMT(2), FORMF(2)
      CHARACTER*20            :: IDSTR2(6)
      INTEGER                 :: GTYPEDUM = 0
!
      INTEGER               :: ncid, iret, ndims, parents, nvars,     &
                               NDIMSGRID, NDIMSVAR, natts, varidtmp, varidf(50)
      INTEGER               :: numDimsLon, numAttsLon
      INTEGER               :: DIMSVAR(4)
       INTEGER               :: dimln(5),dimid(5) 
      CHARACTER*20          :: dimname(5)
      CHARACTER*50          :: timeunits, calendar
      DOUBLE PRECISION      :: dateref, datecur
      REAL                  :: scur, mncur, hcur, fillvalue, TIMESHIFT
      INTEGER               :: I,K,L,ITIME, dref, mref, yref, dcur, mcur, ycur
      INTEGER               :: hhref,mmref,ssref
!
! Variables used in tidal analysis
      REAL(KIND=8), ALLOCATABLE :: ALLTIMES(:)
      INTEGER                   :: TIDE_NDEF, TIDE_ITREND
      REAL(KIND=8), ALLOCATABLE ::       SDEV0(:), SDEV(:), RMSR(:), &
                                         RMSR0(:), RMSRP(:), RESMAX(:)
      REAL ,        ALLOCATABLE ::       TIDE_DATA_ALL(:,:,:)
      REAL ,        ALLOCATABLE ::       SSQ(:), RES(:)
      INTEGER ,        ALLOCATABLE ::       IMAX(:)
      REAL                      :: TIDE_LAT, TIDE_LON
      CHARACTER*256             :: TIDE_file
      INTEGER                   :: TIDEFLAG
      CHARACTER                 :: TIDECONSTNAMES*1024
      CHARACTER(LEN=100)         :: LIST(70)
      LOGICAL                   :: TIDESWITCH, DISTSWITCH 

!/TIDET      INTEGER, PARAMETER      :: LRB = 4
!/TIDET      INTEGER(KIND=8)         :: RPOS
!/TIDET      INTEGER                 ::   LRECL, NREC
!/TIDET      REAL(KIND=LRB), ALLOCATABLE :: NULLBUFF(:)
!/TIDET      CHARACTER(LEN=21)       :: FNAMETXT
!/MPI      INTEGER                   :: STAT_MPI(MPI_STATUS_SIZE)
!
      EQUIVALENCE              ( NXI , NXJ(1) ) , ( NYI , NYJ(1) )
!/
!/ ------------------------------------------------------------------- /
!/
      DATA IDSTR1 / 'ice thickness ' , 'ice viscosity' ,              &
                    'ice density   ' , 'ice modulus  ' ,              &
                    'ice flow diam.' , 'mud density  ' ,              &
                    'mud thickness ', 'mud viscosity ' ,              &
                    'ice conc.     ' , 'water levels ' ,              &
                    'winds         ' , 'currents     ' ,              &
                    'data          ' /
      DATA IDSTR2 / 'pre-processed file  ' , 'long.-lat. grid     ' , &
                    'grid from file (1)  ' , 'grid from file (2)  ' , &
                    'data (assimilation) ' , 'pre-pro. file + tide'  /
      DATA IDSTR3 / 'mean parameters', '1D spectra     ',             &
                    '2D spectra     ' /
!
!/NCO/!     CALL W3TAGB('WAVEPREP',1998,0007,0050,'NP21   ')
!
!--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 1.a  Set number of models
!
      CALL W3NMOD ( 1, 6, 6 )
      CALL W3SETG ( 1, 6, 6 )
!/NL1      CALL W3NAUX (    6, 6 )
!/NL1      CALL W3SETA ( 1, 6, 6 )
      CALL W3NOUT (    6, 6 )
      CALL W3SETO ( 1, 6, 6 )
!
! 1.b  IO set-up.
!
      NDSI   = 10
      NDSO   =  6
      NDSE   =  6
      NDST   =  6
      NDSM   = 11
      NDSDAT = 12
!/O15      NDSTIME = 13
!
      NDSTRC =  6
      NTRACE = 10
      CALL ITRACE ( NDSTRC, NTRACE )
!
!/NCO/!
!/NCO/! Redo according to NCO
!/NCO/!
!/NCO      NDSI   = 11
!/NCO      NDSO   =  6
!/NCO      NDSE   = NDSO
!/NCO      NDST   = NDSO
!/NCO      NDSM   = 12
!/NCO      NDSDAT = 51
!/NCO      NDSTRC = NDSO
!
!/S      CALL STRACE (IENT, 'W3PRNC')
!
!
! 1.c MPP initializations
!
!/SHRD      NAPROC = 1
!/SHRD      IAPROC = 1
!
!/MPI      CALL MPI_INIT      ( IERR_MPI )
!/MPI      CALL MPI_COMM_SIZE ( MPI_COMM_WORLD, NAPROC, IERR_MPI )
!/MPI      CALL MPI_COMM_RANK ( MPI_COMM_WORLD, IAPROC, IERR_MPI )
!/MPI      IAPROC = IAPROC + 1  ! this is to have IAPROC between 1 and NAPROC
!
      IF ( IAPROC .EQ. NAPERR ) THEN
          NDSEN  = NDSE
        ELSE
          NDSEN  = -1
        END IF
!
      IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,900)
!
      J      = LEN_TRIM(FNMPRE)
      OPEN (NDSI,FILE=FNMPRE(:J)//'ww3_prnc.inp',STATUS='OLD',ERR=800,IOSTAT=IERR)
      REWIND (NDSI)
      READ (NDSI,'(A)',END=801,ERR=802,IOSTAT=IERR) COMSTR
      IF (COMSTR.EQ.' ') COMSTR = '$'
      IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,901) COMSTR
!
!--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 2.  Read model definition file.
!
      CALL W3IOGR ( 'READ', NDSM )
      IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,902) GNAME
      ALLOCATE ( IX21(NX,NY), IX22(NX,NY), IY21(NX,NY), IY22(NX,NY),  &
                 JX21(NX,NY), JX22(NX,NY), JY21(NX,NY), JY22(NX,NY),  &
                 MAPOVR(NX,NY) )
      ALLOCATE ( RD11(NX,NY), RD21(NX,NY), RD12(NX,NY), RD22(NX,NY),  &
                 XD11(NX,NY), XD21(NX,NY), XD12(NX,NY), XD22(NX,NY),  &
                 FX(NX,NY), FY(NX,NY), FA(NX,NY),                     &
                 A1(NX,NY), A2(NX,NY), A3(NX,NY) )
!
!--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 3.a Read types and variables from input file.
!
      CALL NEXTLN ( COMSTR , NDSI , NDSE )
      READ (NDSI,*,END=801,ERR=802,IOSTAT=IERR) IDFLD, IDTYPE, FLTIME, FLHDR
!
! 3.b Check types.
!
      FLSTAB = IDFLD .EQ. 'WNS'
      FLBERG = IDFLD .EQ. 'ISI'
      IF        ( IDFLD.EQ.'IC1' ) THEN
          IFLD    = -7
        ELSE IF ( IDFLD.EQ.'IC2' ) THEN
          IFLD    = -6
        ELSE IF ( IDFLD.EQ.'IC3' ) THEN
          IFLD    = -5
        ELSE IF ( IDFLD.EQ.'IC4' ) THEN
          IFLD    = -4
        ELSE IF ( IDFLD.EQ.'IC5' ) THEN
          IFLD    = -3
        ELSE IF ( IDFLD.EQ.'MDN' ) THEN
          IFLD    = -2
        ELSE IF ( IDFLD.EQ.'MTH' ) THEN
          IFLD    = -1
        ELSE IF ( IDFLD.EQ.'MVS' ) THEN
          IFLD    = 0
        ELSE IF ( IDFLD.EQ.'ICE' .OR. FLBERG ) THEN
          IFLD    = 1
        ELSE IF ( IDFLD.EQ.'LEV' ) THEN
          IFLD    = 2
        ELSE IF ( IDFLD.EQ.'WND' .OR. FLSTAB ) THEN
          IFLD    = 3
        ELSE IF ( IDFLD.EQ.'CUR' ) THEN
          IFLD    = 4
        ELSE IF ( IDFLD.EQ.'DAT' ) THEN
          IFLD    = 5
        ELSE
          WRITE (NDSE,1030) IDFLD
          CALL EXTCDE ( 1 )
        END IF
!
      NFCOMP = 1
      TIDEFLAG = 0
      IF (IDFLD.EQ.'DAT') THEN
          ITYPE   = 5
        ELSE IF (IDTYPE.EQ.'AI') THEN
          ITYPE   = 1
        ELSE IF (IDTYPE.EQ.'AT') THEN
          ITYPE   = 6
          TIDEFLAG= 1
          CALL NEXTLN ( COMSTR , NDSI , NDSE )
          READ (NDSI,'(A)',END=801,ERR=803,IOSTAT=IERR) TIDECONSTNAMES
          LIST(:)=''
          CALL STRSPLIT(TIDECONSTNAMES,LIST)
        ELSE IF (IDTYPE.EQ.'LL') THEN
          ITYPE   = 2
        ELSE IF (IDTYPE.EQ.'F1') THEN
          ITYPE   = 3
        ELSE IF (IDTYPE.EQ.'F2') THEN
          ITYPE   = 4
          NFCOMP = 2
        ELSE
          WRITE (NDSE,1031) IDTYPE
          CALL EXTCDE ( 2 )
        END IF
!
!/T      IF (ITYPE.NE.1 .AND. ITYPE.NE.5) WRITE (NDST,9000) ACC
!
      IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,930) IDSTR1(IFLD), IDSTR2(ITYPE)
      IF ( ITYPE.NE.1 .AND. ITYPE.NE.6 ) THEN
!/WNT0          IF ( IAPROC .EQ. NAPOUT .AND.IFLD.EQ.3) WRITE (NDSO,1930)
!/WNT1          IF ( IAPROC .EQ. NAPOUT .AND.IFLD.EQ.3) WRITE (NDSO,1930)
!/WNT2          IF ( IAPROC .EQ. NAPOUT .AND.IFLD.EQ.3) WRITE (NDSO,2930)
!/CRT1          IF ( IAPROC .EQ. NAPOUT .AND.IFLD.EQ.4) WRITE (NDSO,1930)
!/CRT2          IF ( IAPROC .EQ. NAPOUT .AND.IFLD.EQ.4) WRITE (NDSO,2930)
        END IF
      IF ( IAPROC .EQ. NAPOUT .AND.FLBERG ) WRITE (NDSO,938)
      IF ( IAPROC .EQ. NAPOUT .AND.FLSTAB ) WRITE (NDSO,939)
!
      CALL NEXTLN ( COMSTR , NDSI , NDSE )
      READ (NDSI,'(A)',END=801,ERR=802,IOSTAT=IERR) STRDIMSNAME
!
      FIELDSNAME(:)=''
      DIMSNAME(:)=''
      CALL STRSPLIT(STRDIMSNAME,DIMSNAME)
      ! Counts the number of dimensions
      NIDIMS=0
      DO I=1,2
        IF (LEN_TRIM(DIMSNAME(I)).NE.0) NIDIMS=NIDIMS+1
      ENDDO
!
      CALL NEXTLN ( COMSTR , NDSI , NDSE )
      READ (NDSI,'(A)',END=801,ERR=802,IOSTAT=IERR) STRFIELDSNAME
!
      FIELDSNAME(:)=''
      CALL STRSPLIT(STRFIELDSNAME,FIELDSNAME)
      ! Counts the number of variables
      NFIELDS=0
      DO WHILE (len_trim(FIELDSNAME(NFIELDS+1)).NE.0)
        NFIELDS=NFIELDS+1
      ENDDO
!
! 3.c Additional input for format types and time
! ... time
!
      IF (.NOT. FLTIME) THEN
        CALL NEXTLN ( COMSTR , NDSI , NDSE )
        READ (NDSI,*,END=801,ERR=802,IOSTAT=IERR) TIMESTART
        IF (TIMESTART(1).LT.10000000) THEN
          WRITE (NDSE,1035) TIME
          CALL EXTCDE ( 4 )
        END IF
        CALL STME21 ( TIMESTART , IDTIME )
        IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,931) IDTIME
      END IF
!
      J = 1
      IF ( FLAGLL ) THEN
        FACTOR = 1.
      ELSE
        FACTOR = 1.E-3
      END IF

!
! 3.d Read Input netcdf file
!
      CALL NEXTLN ( COMSTR , NDSI , NDSE )
      READ (NDSI,*,END=801,ERR=802,IOSTAT=IERR) NAMEF
      IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,967) NAMEF

      J = LEN_TRIM(FNMPRE)
      iret=nf90_open(path=FNMPRE(:J)//NAMEF,mode=nf90_noWRITE,ncid=ncid)
      CALL check_err(iret)
!
! gets requested variables ids and dimensions
!
      DO I=1,NFIELDS
        iret = nf90_inq_varid(ncid,TRIM(FIELDSNAME(I)),varidf(I))
        CALL check_err(iret)
        iret = nf90_inquire_variable(ncid, varidf(I), ndims=NDIMSVAR,  &
                                     natts = nAtts)
        CALL check_err(iret)
        iret = nf90_inquire_variable(ncid, varidf(I),                  &
                                     dimids=DIMSVAR(:NDIMSVAR))
        CALL check_err(iret)
! Now the dimensions
        DO J=1,NDIMSVAR
          iret=nf90_inquire_dimension(ncid,DIMSVAR(J),name=dimname(J), &
                                      len=dimln(J))
          CALL check_err(iret)
          END DO
        END DO
!
! instanciates time
!
      yref=0.
      mref=0.
      dref=0.
      hhref=0.
      mmref=0.
      ssref=0.
      iret=nf90_inq_varid(ncid,"time",varidtmp)
      IF ( iret/=nf90_noerr ) iret=nf90_inq_varid(ncid,"MT",varidtmp)
      CALL check_err(iret)
      iret=nf90_get_att(ncid,varidtmp,"calendar",calendar)
      IF ( iret/=nf90_noerr ) THEN 
        WRITE(NDSE,1028)
      ELSE IF ((INDEX(calendar, "julian").EQ.0) .AND. &
               (INDEX(calendar, "gregorian").EQ.0)) THEN
        WRITE(NDSE,1029)
      END IF
!
      iret=nf90_get_att(ncid,varidtmp,"units",timeunits)
      IF (INDEX(timeunits, "seconds").NE.0) THEN
        read(timeunits(15:18),'(I4.4)',END=804,ERR=805,IOSTAT=IERR) yref
        read(timeunits(20:21),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) mref
        read(timeunits(23:24),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) dref
        read(timeunits(26:27),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) hhref
        read(timeunits(29:30),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) mmref
        read(timeunits(32:33),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) ssref
      ELSE IF (INDEX(timeunits, "days").NE.0) THEN
        read(timeunits(12:15),'(I4.4)',END=804,ERR=805,IOSTAT=IERR) yref
        read(timeunits(17:18),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) mref
        read(timeunits(20:21),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) dref
        read(timeunits(23:24),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) hhref
        read(timeunits(26:27),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) mmref
        read(timeunits(29:30),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) ssref
      ELSE IF (INDEX(timeunits, "hours").NE.0) THEN
        read(timeunits(13:16),'(I4.4)',END=804,ERR=805,IOSTAT=IERR) yref
        read(timeunits(18:19),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) mref
        read(timeunits(21:22),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) dref
        read(timeunits(24:25),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) hhref
        read(timeunits(27:28),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) mmref
        read(timeunits(30:31),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) ssref
      ELSE IF (INDEX(timeunits, "minutes").NE.0) THEN
        read(timeunits(15:18),'(I4.4)',END=804,ERR=805,IOSTAT=IERR) yref
        read(timeunits(20:21),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) mref
        read(timeunits(23:24),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) dref
        read(timeunits(26:27),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) hhref
        read(timeunits(29:30),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) mmref
        read(timeunits(32:33),'(I2.2)',END=804,ERR=805,IOSTAT=IERR) ssref
      END IF
      dateref=JULDAY(dref,mref,yref) + DBLE(hhref)/24. + DBLE(mmref)/1440. + DBLE(ssref)/86400.
!
! instanciates generic variables dimensions
!
      NXI=0
      NYI=0 
      NDIMSGRID=2
      DO i=1,NDIMSVAR
        IF (dimname(i) .EQ. "time".OR.dimname(i) .EQ."MT") NTI = dimln(i)
        IF (dimname(i) .EQ. DIMSNAME(1)) NXI = dimln(i)
        IF (dimname(i) .EQ. DIMSNAME(1).AND.NIDIMS.EQ.1) THEN
          NDIMSGRID=1
          NYI = 1
          END IF
        IF (NIDIMS.GE.2) THEN 
          IF (dimname(i) .EQ. DIMSNAME(2)) NYI = dimln(i)
          END IF
        END DO
!
      IF (NXI*NYI.EQ.0) GOTO 864
!
      IF (ITYPE.NE.1.AND.ITYPE.NE.6) THEN
        ALLOCATE (ALA(NXI,NYI))
        ALLOCATE (ALO(NXI,NYI))
        iret=nf90_inq_varid(ncid,"longitude",varidtmp)
        IF ( iret/=nf90_noerr ) iret=nf90_inq_varid(ncid,"lon",varidtmp)
        IF ( iret/=nf90_noerr ) iret=nf90_inq_varid(ncid,"Longitude",varidtmp)
        IF ( iret/=nf90_noerr ) iret=nf90_inq_varid(ncid,"x",varidtmp)
        IF ( iret/=nf90_noerr ) iret=nf90_inq_varid(ncid,"X",varidtmp)
        IRET = nf90_inquire_variable(ncid, varidtmp, ndims = numDimsLon, natts = numAttsLon)
        call check_err(iret)
        IF (numDimsLon.EQ.1) THEN
          iret=nf90_get_var(ncid,varidtmp,X0I,start=(/1/))
          call check_err(iret)
          iret=nf90_get_var(ncid,varidtmp,XNI,start=(/NXI/))
          call check_err(iret)
          iret=nf90_get_var(ncid,varidtmp,ALO(:,1))
          call check_err(iret)
          DO i=1,NYI
            ALO(:,i)=ALO(:,1)
            END DO
        ELSE 
          iret=nf90_get_var(ncid,varidtmp,X0I,start=(/1,1/))
          call check_err(iret)
          iret=nf90_get_var(ncid,varidtmp,XNI,start=(/NXI,1/))
          call check_err(iret)
          iret=nf90_get_var(ncid,varidtmp,ALO(:,:))
          call check_err(iret)
          END IF
!
        iret=nf90_inq_varid(ncid,"latitude",varidtmp)
        IF ( iret/=nf90_noerr ) iret=nf90_inq_varid(ncid,"lat",varidtmp)
        IF ( iret/=nf90_noerr ) iret=nf90_inq_varid(ncid,"Latitude",varidtmp)
        IF ( iret/=nf90_noerr ) iret=nf90_inq_varid(ncid,"y",varidtmp)
        IF ( iret/=nf90_noerr ) iret=nf90_inq_varid(ncid,"Y",varidtmp)
        iret = nf90_inquire_variable(ncid, varidtmp, ndims = numDimsLon, natts = numAttsLon)
        CALL check_err(iret)
        iret=nf90_get_var(ncid,varidtmp,Y0I, start=(/1/))
        CALL check_err(iret)
        IF (numDimsLon.EQ.1) THEN
          iret=nf90_get_var(ncid,varidtmp,ALA(1,:))
          CALL check_err(iret)
          YNI=ALA(1,NYI)
          DO i=1,NXI
            ALA(i,:)=ALA(1,:)
            END DO
        ELSE 
          iret=nf90_get_var(ncid,varidtmp,ALA(:,:))
          CALL check_err(iret)
          YNI=ALA(1,NYI)
          END IF
        END IF
        IF ((ITYPE.EQ.2) .AND. (GTYPE .EQ. RLGTYPE) .AND. (Y0I.GT.YNI)) THEN
          WRITE (NDSE,1032)
          CALL EXTCDE ( 2 )
        END IF

!
! gets FillValue
!
      iret=nf90_inq_varid(ncid,FIELDSNAME(1),varidtmp)
      CALL check_err(iret)
      iret=nf90_get_att(ncid,varidtmp,"_FillValue", fillvalue)
!
! ... type 1: "As Is" (AI)
!
      IF (ITYPE.EQ.1.OR.ITYPE.EQ.6) THEN
!
        NXI    = NX
        NYI    = NY
        ALLOCATE ( MASK(NXI,NYI) )
        MASK   = 1
        IF(GTYPE .EQ. UNGTYPE) THEN
!
!  X0, Y0 are the coordinates of the lower-left point in mesh
!
          RW(1) = FACTOR*X0        ; RW(2) = FACTOR*MAXX
          RW(3) = FACTOR*Y0        ; RW(4) = FACTOR*MAXY
        ELSE
          RW(1) = FACTOR*XGRD(1,1) ; RW(2) = FACTOR*XGRD(NY,NX)
          RW(3) = FACTOR*YGRD(1,1) ; RW(4) = FACTOR*YGRD(NY,NX)
        END IF
        IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,932) NXI, NYI
        IF ( FLAGLL ) THEN
          IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,933) RW(1),RW(2),RW(3),RW(4)
        ELSE
          IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,733) RW(1),RW(2),RW(3),RW(4)
        END IF
!
! ... type 2
!
      ELSE IF (ITYPE.EQ.2) THEN
!
        IF (NXI.LT.2 .OR. NYI.LT.2) THEN
          WRITE (NDSE,1036) NXI, NYI
          CALL EXTCDE ( 5 )
        END IF
        ALLOCATE ( MASK(NXI,NYI) )
        MASK   = 1
        IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,932) NXI, NYI

        IF ( FLAGLL ) THEN
          IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,933) FACTOR*X0I, FACTOR*XNI,                &
                           FACTOR*Y0I, FACTOR*YNI
        ELSE
          IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,733) FACTOR*X0I, FACTOR*XNI,                &
                           FACTOR*Y0I, FACTOR*YNI
        END IF
!
! ... type 5
!
      ELSE IF (ITYPE.EQ.5) THEN
        CALL NEXTLN ( COMSTR , NDSI , NDSE )
        READ (NDSI,*,END=801,ERR=802,IOSTAT=IERR)                   &
              DATTYP, RECLDT, NODATA
        IF (DATTYP.LT.0 .OR. DATTYP.GT.2) THEN
          WRITE (NDSE,1033) DATTYP
          CALL EXTCDE ( 6 )
        END IF
        IF (RECLDT.LE.0) THEN
          WRITE (NDSE,1034) RECLDT
          CALL EXTCDE ( 7 )
        END IF
        IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,934) IDSTR3(DATTYP+1), RECLDT, NODATA
        WRITE (IDFLD,935) DATTYP
        DEALLOCATE ( IX21, IX22, IY21, IY22, JX21, JX22, JY21, JY22, &
                     MAPOVR )
        DEALLOCATE ( RD11, RD21, RD12, RD22, XD11, XD21, XD12, XD22, &
                     FX, FY, FA, A1, A2, A3 )
!
! ... types 3 and 4 ... in preprocessing loop ....
!
      END IF
!
!--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 4   Prepare interpolation.
!
      IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,940)
!
      IF (ITYPE.NE.1 .AND. ITYPE.NE.5 .AND. ITYPE.NE.6 ) THEN
!
! 4.a Longitude - latitude grid
!
        IF (ITYPE.EQ.2) THEN
          IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,941)
!
! ... setup coordinates
!
          SXI    = (XNI-X0I)/REAL(NXI-1)
          SYI    = (YNI-Y0I)/REAL(NYI-1)
          ICLO   = ICLOSE_NONE
          IF ( FLAGLL ) THEN
            IF ( ABS(ABS(REAL(NXI)*SXI)-360.) .LT. 0.1*ABS(SXI) ) THEN
              ICLO = ICLOSE_SMPL
            END IF
          END IF
!
! ... create grid search utility
!
          GSI = W3GSUC( .TRUE., FLAGLL, ICLO, ALO, ALA )
!
! ... construct Interpolation data
!
!/T1              WRITE (NDST,9045)
          IF (GTYPE .NE. UNGTYPE) THEN
            DO IY=1,NY
              DO IX=1,NX
                INGRID = W3GRMP( GSI, XGRD(IY,IX), YGRD(IY,IX),      &
                                 IS, JS, RW )
                IF ( .NOT.INGRID ) THEN
                    IF ( IAPROC .EQ. NAPOUT ) WRITE(NDSO,1042) IX, IY, XGRD(IY,IX), YGRD(IY,IX)
                  IX21(IX,IY) = 1
                  IX22(IX,IY) = 1
                  IY21(IX,IY) = 1
                  IY22(IX,IY) = 1
                  RD11(IX,IY) = 0.
                  RD21(IX,IY) = 0.
                  RD12(IX,IY) = 0.
                  RD22(IX,IY) = 0.

                  CYCLE
                  END IF
                IX21(IX,IY) = IS(1)
                IX22(IX,IY) = IS(2)
                IY21(IX,IY) = JS(1)
                IY22(IX,IY) = JS(4)
                RD11(IX,IY) = RW(1)
                RD21(IX,IY) = RW(2)
                RD12(IX,IY) = RW(4)
                RD22(IX,IY) = RW(3)
!/T1                  WRITE (NDST,9046) IX, IY,                           &
!/T1                    IX21(IX,IY),IX22(IX,IY),IY21(IX,IY),IY22(IX,IY),  &
!/T1                    RD11(IX,IY),RD12(IX,IY),RD21(IX,IY),RD22(IX,IY)
                END DO
              END DO
          ELSE ! GTYPE .NE. UNGTYPE
            DO IX=1, NX
              X = XYB(IX,1)
              Y = XYB(IX,2)
              IX21(IX,1) =   1 + INT(MOD(360.+(X-X0I),360.)/SXI)
!
! Manages the simple closure of the grid
!
              IF (ICLO.EQ.ICLOSE_NONE) THEN 
                IF (IX21(IX,1).LT.1.OR.IX21(IX,1).GT.NXI-1) WRITE(NDSO,1042) IX, IY, X, Y 
                IX21(IX,1) =   MAX ( 1 , MIN(IX21(IX,1),NXI-1) )
                IX22(IX,1) =   IX21(IX,1) + 1
              ELSE
                IX21(IX,1) =   MAX ( 1 , MIN(IX21(IX,1),NXI) )
                IX22(IX,1) =   MOD(IX21(IX,1),NXI)+1
                END IF 
              IY21(IX,1) =   1 + INT((Y-Y0I)/SYI)
              IF (IY21(IX,1).LT.1.OR.IY21(IX,1).GT.NYI-1) WRITE(NDSO,1042) IX, IY, X, Y
              IY21(IX,1) =   MAX ( 1 , MIN(IY21(IX,1),NYI-1) )
              IY22(IX,1) =   IY21(IX,1) + 1
!
              RW(1) =  MOD(360.+(X-X0I),360.)/SXI - REAL(IX21(IX,1)-1)
              RW(2) =  (Y-Y0I)/SYI - REAL(IY21(IX,1)-1)
!       
              IF (IX21(IX,1).LE.1 .AND. RW(1).LT.ACC) THEN
                IF (RW(1).LT.0.) THEN
                  RW(1) = 0.
                  WRITE (NDSO,1043) X
!/T                  FLMOD  = .TRUE.
                  END IF
                END IF
!
              IF (IX21(IX,1).GE.(NXI-1) .AND. RW(1).GT.1.-ACC) THEN
                IF (RW(1).GT.1.) THEN
                  WRITE (NDSO,1043) X
                  RW(1) = 1.
!/T                       FLMOD  = .TRUE.
                  END IF
                END IF
!
              IF (IY21(IX,1).LE.1 .AND. RW(2).LT.ACC) THEN
                IF (RW(2).LT.0.) THEN
                  WRITE (NDSO,1044) Y
                  RW(2) = 0.
!/T                  FLMOD  = .TRUE.
                END IF
              END IF
!
              IF (IY21(IX,1).GE.NYI .AND. RW(2).GT.1.-ACC) THEN
                IF (RW(2).GT.1) THEN
                  WRITE (NDSO,1044) Y
                  RW(2) = 1.
!/T                       FLMOD  = .TRUE.
                END IF
              END IF
!       
              EFAC = SQRT ( MAX(0.,ABS(RW(1)-0.5)-0.5)**2 +  &
                            MAX(0.,ABS(RW(2)-0.5)-0.5)**2 )
              EFAC        = 1. / ( 1. + 0.25*EFAC**2 )
                
              RD11(IX,1) = EFAC * (1.-RW(1)) * (1.-RW(2))
              RD21(IX,1) = EFAC *    RW(1)   * (1.-RW(2))
              RD12(IX,1) = EFAC * (1.-RW(1)) *    RW(2)
              RD22(IX,1) = EFAC *    RW(1)   *    RW(2)
            END DO ! IX=1, NX
          END IF ! GTYPE .NE. UNGTYPE
!
          CALL W3GSUD( GSI )

!
! 4.b Grid(s) from file
!
        ELSE ! ITYPE.EQ.2
          IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,942)
!
! ... prepare overlay map
!
          DO IY=1, NY
            DO IX=1, NX
              IF ( MAPSTA(IY,IX) .EQ. 0 ) THEN
                  MAPOVR(IX,IY) = ILAND
                ELSE
                  MAPOVR(IX,IY) = 0
              END IF
            END DO
          END DO
!
! ... loop over fields
!
          DO J=1, NFCOMP
!
            IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,943) J
!
! ... file info lat-long file
!
            CALL NEXTLN ( COMSTR , NDSI , NDSE )
            READ (NDSI,*,END=801,ERR=802,IOSTAT=IERR)             &
                  NXJ(J), NYJ(J), CLO(J)
            IF (NXJ(J).LT.2 .OR. NYJ(J).LT.2) THEN
                WRITE (NDSE,1036) NXJ(J), NYJ(J)
                CALL EXTCDE ( 10 )
              END IF
            IF ( ALLOCATED(MASK) ) DEALLOCATE (MASK)
            ALLOCATE ( MASK(NXJ(J),NYJ(J)) )
            MASK   = 1
            IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,944) NXJ(J), NYJ(J), CLO(J)
!
            CALL NEXTLN ( COMSTR , NDSI , NDSE )
            READ (NDSI,*,END=801,ERR=802,IOSTAT=IERR)             &
                  FROMLL, IDLALL, IDFMLL, FORMLL
            IF (IDLALL.LT.1 .OR. IDLALL.GT.4) IDLALL   = 1
            IF (IDFMLL.LT.1 .OR. IDFMLL.GT.3) IDFMLL   = 1
            IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,945) IDLALL, IDFMLL
            IF (IDFMLL.EQ.2) WRITE (NDSO,946) FORMLL
!
            CALL NEXTLN ( COMSTR , NDSI , NDSE )
            READ (NDSI,*,END=801,ERR=802,IOSTAT=IERR) NDSLL, NAMELL
!/NCO                NDSLL  = 20 + NFCOMP
            IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,947) NDSLL
            IF ( IAPROC .EQ. NAPOUT.AND.FROMLL.EQ.'NAME') WRITE (NDSO,948) NAMELL
            IF (NDSLL.EQ.NDSI) THEN
                WRITE (NDSE,1038)
                CALL NEXTLN ( COMSTR , NDSI , NDSE )
            ELSE
!
! ... open lat-long file
!
              IF ( IDFMLL .EQ. 3 ) THEN
                IF (FROMLL.EQ.'NAME') THEN
                  JJ = LEN_TRIM(FNMPRE)
                  OPEN (NDSLL,FILE=FNMPRE(:JJ)//NAMELL,     &
                              FORM='UNFORMATTED',STATUS='OLD',    &
                              ERR=845,IOSTAT=IERR)
                ELSE
                  OPEN (NDSLL, FORM='UNFORMATTED',          &
                               STATUS='OLD',ERR=845,IOSTAT=IERR)
                END IF
              ELSE
                IF (FROMLL.EQ.'NAME') THEN
                  JJ = LEN_TRIM(FNMPRE)
                  OPEN (NDSLL,FILE=FNMPRE(:JJ)//NAMELL,     &
                              STATUS='OLD',ERR=845,IOSTAT=IERR)
                ELSE
                  OPEN (NDSLL,                              &
                              STATUS='OLD',ERR=845,IOSTAT=IERR)
                END IF
              END IF
!
            END IF
!
! ... read lat-lon data
!
            IF ( ASSOCIATED(ALA) ) THEN
                DEALLOCATE ( ALA, ALO )
                NULLIFY ( ALA, ALO )
            END IF
            ALLOCATE ( ALA(NXJ(J),NYJ(J)), ALO(NXJ(J),NYJ(J)) )
            CALL INA2R (ALA, NXJ(J), NYJ(J), 1, NXJ(J), 1, NYJ(J),&
                 NDSLL, NDST, NDSE, IDFMLL, FORMLL, IDLALL, 1., 0.)
            CALL INA2R (ALO, NXJ(J), NYJ(J), 1, NXJ(J), 1, NYJ(J),&
                 NDSLL, NDST, NDSE, IDFMLL, FORMLL, IDLALL, 1., 0.)
!
            IF ( NDSLL .NE. NDSI ) CLOSE (NDSLL)
!
! ... file info mask file
!
            IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,949)
!
            CALL NEXTLN ( COMSTR , NDSI , NDSE )
            READ (NDSI,*,END=801,ERR=802,IOSTAT=IERR)             &
                  FROMLL, IDLALL, IDFMLL, FORMLL
            IF (IDLALL.LT.1 .OR. IDLALL.GT.4) IDLALL   = 1
            IF (IDFMLL.LT.1 .OR. IDFMLL.GT.3) IDFMLL   = 1
            IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,945) IDLALL, IDFMLL
            IF (IDFMLL.EQ.2) WRITE (NDSO,946) FORMLL
!
            CALL NEXTLN ( COMSTR , NDSI , NDSE )
            READ (NDSI,*,END=801,ERR=802,IOSTAT=IERR) NDSLL, NAMELL
!/NCO                NDSLL  = 22 + NFCOMP
            IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,947) NDSLL
            IF (FROMLL.EQ.'NAME') WRITE (NDSO,948) NAMELL
            IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,*) ' '
            IF (NDSLL.EQ.NDSI) THEN
                WRITE (NDSE,1038)
                CALL NEXTLN ( COMSTR , NDSI , NDSE )
            ELSE
!
! ... open mask file
!
              IF ( IDFMLL .EQ. 3 ) THEN
                IF (FROMLL.EQ.'NAME') THEN
                  JJ = LEN_TRIM(FNMPRE)
                  OPEN (NDSLL,FILE=FNMPRE(:JJ)//NAMELL,     &
                            FORM='UNFORMATTED',STATUS='OLD',    &
                            ERR=846,IOSTAT=IERR)
                ELSE
                  OPEN (NDSLL,FORM='UNFORMATTED',           &
                              STATUS='OLD',ERR=846,IOSTAT=IERR)
                END IF
              ELSE
                IF (FROMLL.EQ.'NAME') THEN
                  JJ = LEN_TRIM(FNMPRE)
                  OPEN (NDSLL,FILE=FNMPRE(:JJ)//NAMELL,     &
                              STATUS='OLD',ERR=846,IOSTAT=IERR)
                ELSE
                  OPEN (NDSLL,                              &
                              STATUS='OLD',ERR=846,IOSTAT=IERR)
                END IF
              END IF
!
            END IF
!
! ... read mask data
!
            CALL INA2I (MASK, NXJ(J), NYJ(J), 1,NXJ(J), 1,NYJ(J), &
                 NDSLL, NDST, NDSE, IDFMLL, FORMLL, IDLALL, 1, 0)
            IF ( NDSLL .NE. NDSI ) CLOSE (NDSLL)
!
!/T1a                WRITE (NDST,9050)
!/T1a                DO IY=1, NYJ(J)
!/T1a                  DO IX=1,NXJ(J)
!/T1a                    WRITE (NDST,9051) IX, IY, ALA(IX,IY),             &
!/T1a                           ALO(IX,IY), MASK(IX,IY)
!/T1a                    END DO
!/T1a                  END DO
!
! ... generate interpolation data
!
            IF ( J .EQ. 1 ) THEN
              CALL W3FLDP ( NDSO, NDST, NDSE, IERR, FLAGLL,       &
                NX, NY, NX, NY, YGRD, XGRD, MAPOVR, ILAND,        &
                NXJ(J), NYJ(J), NXJ(J), NYJ(J), CLO(J), ALA, ALO, &
                MASK, RD11, RD21, RD12, RD22, IX21, IX22, IY21,   &
                IY22 )
            ELSE
              CALL W3FLDP ( NDSO, NDST, NDSE, IERR, FLAGLL,       &
                NX, NY, NX, NY, YGRD, XGRD, MAPOVR, ILAND,        &
                NXJ(J), NYJ(J), NXJ(J), NYJ(J), CLO(J), ALA, ALO, &
                MASK, XD11, XD21, XD12, XD22, JX21, JX22, JY21,   &
                JY22 )
            END IF ! J .EQ. 1 
!
          END DO ! J=1, NFCOMP
!
! ... average two fields !
!
          IF ( NFCOMP .EQ. 2) THEN
            DO IX=1, NX
              DO IY=1, NY
                IF ( MAPOVR(IX,IY) .GE. 2) THEN
                  FACTOR = 1. / REAL(MAPOVR(IX,IY))
                  RD11(IX,IY) = FACTOR * RD11(IX,IY)
                  RD12(IX,IY) = FACTOR * RD12(IX,IY)
                  RD21(IX,IY) = FACTOR * RD21(IX,IY)
                  RD22(IX,IY) = FACTOR * RD22(IX,IY)
                  XD11(IX,IY) = FACTOR * XD11(IX,IY)
                  XD12(IX,IY) = FACTOR * XD12(IX,IY)
                  XD21(IX,IY) = FACTOR * XD21(IX,IY)
                  XD22(IX,IY) = FACTOR * XD22(IX,IY)
                END IF
              END DO
            END DO
          END IF ! NFCOMP .EQ. 2
!
        END IF ! ITYPE.EQ.2
      END IF ! ITYPE.NE.1 .AND. ITYPE.NE.5
!
! 4.c Input location and format
!
      DO J=1, NFCOMP
!
        IF ( ITYPE .EQ. 5 ) THEN
            IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,960)
          ELSE
            IF (ITYPE.LE.3) THEN
                IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,961) NXJ(J), NYJ(J)
              ELSE
                IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,962) J, NXJ(J), NYJ(J)
              END IF
          END IF ! ITYPE .EQ. 5
!
      END DO
!
!--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 5   Prepare files
!
      IF ( NFCOMP .EQ. 1 ) THEN
          NXJ  (2) = NXJ  (1)
          NYJ  (2) = NYJ  (1)
          NDSF (2) = NDSF (1)
          IDLAF(2) = IDLAF(1)
          IDFMF(2) = IDFMF(1)
          FORMT(2) = FORMT(1)
          FORMF(2) = FORMF(1)
        END IF

! 5.b Open and prepare output file
!
      IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,971)
      J      = LEN_TRIM(FNMPRE)
      IF ( ITYPE .LE. 4 .OR. ITYPE.EQ.6 ) THEN
!/TIDE          IF (ITYPE.EQ.6) THEN 
!/TIDE            CALL VUF_SET_PARAMETERS
!/TIDE            TIDE_NDEF = NFIELDS
!/TIDE            CALL TIDE_FIND_INDICES_ANALYSIS(LIST)
!/TIDE            END IF
            IF ( IAPROC .EQ. NAPOUT )  &
               CALL W3FLDO ( 'WRITE', IDFLD, NDSDAT, NDST, NDSE,      &
               NX, NY, GTYPE, IERR, FPRE=FNMPRE(:J),                  &
               FHDR=FLHDR, TIDEFLAGIN=TIDEFLAG)
        ELSE
          IF ( IAPROC .EQ. NAPOUT )  &
          CALL W3FLDO ( 'WRITE', IDFLD, NDSDAT, NDST, NDSE,           &
               RECLDT, 0, GTYPEDUM, IERR, FPRE=FNMPRE(:J) )
        END IF
!/TIDET   IF (TIDEFLAG.GT.0) THEN  
!/TIDET   LRECL  = TIDE_MF*LRB*NFIELDS*2
!/TIDET   NREC = LRECL / LRB
!/TIDET   ALLOCATE(NULLBUFF(NREC))
!/TIDET   NULLBUFF(1:NREC) = 0.

!/TIDET   OPEN (990,FILE='tidana.dat',FORM='UNFORMATTED',       &
!/TIDET                 ACCESS='STREAM')
!/TIDET   FNAMETXT  = 'tidanaNNN.txt'
!/TIDET          WRITE (FNAMETXT(7:9),'(I3.3)') IAPROC
!/TIDET   OPEN (989,FILE=FNAMETXT,status='unknown')
!/TIDET   ENDIF 
!
! 5.c Initialize fields
!
      IF ( ITYPE .NE. 5 ) THEN
          FX = 0.
          FY = 0.
          FA = 0.
          MXM = MAX ( NXJ(1), NXJ(2) )
          MYM = MAX ( NYJ(1), NYJ(2) )
          IF (ITYPE.EQ.1.AND.GTYPE.EQ.UNGTYPE) THEN 
            ALLOCATE ( XC(MXM,1), YC(MXM,1), AC(MXM,1) )
          ELSE
            ALLOCATE ( XC(MXM,MYM), YC(MXM,MYM), AC(MXM,MYM) )
          END IF
          XC = 0.
          YC = 0.
          AC = 0.
        END IF
!
! For ITYPE.LE.5 time stemps are read one by one. 
! Otherwise points are read one by one for tidal analysis
!
      TIDESWITCH=.FALSE.
!/TIDE      TIDESWITCH=.TRUE.
      IF (ITYPE.GE.6.AND.TIDESWITCH) THEN
! 
! Reads in the full time vector
!
        IF (NX*NY.GT.4000) THEN 
!/MPI     IF ((NX*NY)/NAPROC.LT.4000) THEN 
!/MPI        IF (IAPROC.EQ.NAPOUT) WRITE(NDSE,*) 'Starting tidal analysis ... '
!/MPI     ELSE
!/MPI       IF (IAPROC.EQ.NAPOUT) WRITE(NDSE,*) 'Starting tidal analysis for ',NX*NY,     &
!/MPI                         ' points. This can take hours ...'
!/MPI     ENDIF 
!/MPI     IF (NX*NY.LT.4000) THEN 
          WRITE(NDSE,*) 'Starting tidal analysis for ',NX*NY, ' points.'
          IF (NAPROC.EQ.1) WRITE(NDSE,*) 'This can take hours ...Consider running this with MPI '
          END IF
!/MPI     END IF
        iret=nf90_inq_varid(ncid,"time",varidtmp)
        IF ( iret/=nf90_noerr ) iret=nf90_inq_varid(ncid,"MT",varidtmp)
        CALL check_err(iret)
        ALLOCATE(ALLTIMES(NTI))
        iret=nf90_get_var(ncid,varidtmp,ALLTIMES,start=(/1/))
        CALL check_err(iret)
        IF (INDEX(timeunits, "seconds").NE.0)   ALLTIMES=ALLTIMES/86400.
        IF (INDEX(timeunits, "minutes").NE.0)   ALLTIMES=ALLTIMES/1440.
        IF (INDEX(timeunits, "hours").NE.0)   ALLTIMES=ALLTIMES/24.
        ALLTIMES=dateref+ALLTIMES
        CALL CALDAT(INT(ALLTIMES(1)),dcur,mcur,ycur)
        hcur=(ALLTIMES(1)-INT(ALLTIMES(1)))*24.
        IF (ABS(hcur-NINT(hcur)).LT.1.E-5) hcur=NINT(hcur)
        mncur=(hcur-INT(hcur))*60.
        IF (ABS(mncur-NINT(mncur)).LT.1.E-5) mncur=NINT(mncur)
        scur=(mncur-INT(mncur))*60.
        IF (ABS(scur-NINT(scur)).LT.1.E-5) scur=NINT(scur)
! 
! Performs tidal analysis
!
!/TIDE        TIDE_NTI = NTI 
!/TIDE        TIDE_NDEF = NFIELDS
!/TIDE        ALLOCATE(SDEV0(TIDE_NDEF),SDEV(TIDE_NDEF), RMSR(TIDE_NDEF), &
!/TIDE                   RES(TIDE_NDEF), SSQ(TIDE_NDEF),RMSR0(TIDE_NDEF), &
!/TIDE                 RMSRP(TIDE_NDEF), IMAX(TIDE_NDEF), RESMAX(TIDE_NDEF))

!/TIDE        ALLOCATE( TIDE_DATA(TIDE_NTI,TIDE_NDEF) )  
!/TIDE        ALLOCATE( TIDE_DAYS(TIDE_NTI), TIDE_SECS(TIDE_NTI), TIDE_HOURS(TIDE_NTI) )  
!/TIDE        ALLOCATE(V_ARG(170,TIDE_NTI),F_ARG(170,TIDE_NTI),U_ARG(170,TIDE_NTI))
!/TIDE        TIDE_NX=NX
!/TIDE        TIDE_NY=NY
!/TIDE        ALLOCATE(TIDAL_CONST(NX,NY,TIDE_MF,2,2))
!/TIDE        TIDAL_CONST(:,:,:,:,:)=0.
!/TIDE        DO I=1,NFIELDS
!/TIDE          iret=nf90_inq_varid(ncid,FIELDSNAME(I),varidf(I))
!/TIDE        END DO
!/TIDE        iret=nf90_get_att(ncid,varidtmp,"_FillValue", fillvalue)

!/TIDE       IRET = NF90_INQUIRE_ATTRIBUTE(NCID,VARIDF(1),'scale_factor')
!/TIDE             IF ( IRET .EQ. NF90_NOERR ) THEN
!/TIDE          IRET = NF90_GET_ATT(NCID,VARIDF(1),'scale_factor',XCFAC)
!/TIDE        ENDIF
!/TIDE              IRET = NF90_INQUIRE_ATTRIBUTE(NCID,VARIDF(1),'add_offset')
!/TIDE              IF ( IRET .EQ. NF90_NOERR) THEN
!/TIDE                 IRET = NF90_GET_ATT(NCID,VARIDF(1),'add_offset',XCOFF)
!/TIDE              END IF
!/TIDE              IF ( NFCOMP.EQ.2 .OR. IFLD.GE.3 .OR. FLBERG ) THEN
!/TIDE                IRET = NF90_INQUIRE_ATTRIBUTE(NCID,VARIDF(2),'scale_factor')
!/TIDE                IF (IRET .EQ. NF90_NOERR) THEN
!/TIDE                   IRET = NF90_GET_ATT(NCID,VARIDF(2),'scale_factor',YCFAC)
!/TIDE                ENDIF
!/TIDE                IRET = NF90_INQUIRE_ATTRIBUTE(NCID,VARIDF(2),'add_offset')
!/TIDE                IF (IRET .EQ. NF90_NOERR) THEN
!/TIDE         IRET = NF90_GET_ATT(NCID,VARIDF(2),'add_offset',YCOFF)
!/TIDE                END IF
!/TIDE              END IF

!
! Loop on grid points: reads all points for same IY to limit I/O time
! this may require too much memory
! Modified on 12-Feb-2013: now reads one block from IX to IX+NXLM-1
!       NSEAL  = 1 + (NSEA-IAPROC)/NAPROC
        ALLOCATE(NXL(NAPROC))
        DO JPC=1,NAPROC
          NXL(JPC)  = 1 + (NX-JPC)/NAPROC
          END DO
        NXLM=NXL(1)   ! this is the max value of NXL
        NXL(1:NAPROC-1) = NXLM
        NXL(NAPROC)=NX-(NAPROC-1)*NXLM

        ALLOCATE(TIDE_DATA_ALL(NXL(IAPROC),NTI,NFIELDS))


        WRITE(NDSE,*) "Number of points for this processor :", IAPROC,NXL(IAPROC)
              DISTSWITCH = .FALSE.
!/DIST        DISTSWITCH = .TRUE.
!/TIDE        IF (DISTSWITCH) THEN 
!/TIDE          TIDESIZ = MAXVAL(NXL) * TIDE_MF * NFIELDS * 2
!/MPI          CALL MPI_TYPE_VECTOR ( TIDESIZ, 1, 1, MPI_REAL, TIDALSTUFF, IERR_MPI )
!/MPI          CALL MPI_TYPE_COMMIT ( TIDALSTUFF, IERR_MPI )
!/TIDE          ALLOCATE(TIDALSTUFFA(TIDESIZ))
!/TIDE          TIDALSTUFFA(:)=0.
!/TIDE          END IF
!
! Loops on Y dimension
!
        DO IY=1,NY
          IND=0
!
! big data read split in small data read ... 
          IF (NDIMSGRID.EQ.1) THEN 
            DO I=1,NFIELDS
                iret=nf90_get_var(ncid,varidf(I),TIDE_DATA_ALL(:,:,I), &
                                start=(/1+(IAPROC-1)*NXLM,1/),count=(/NXL(IAPROC),NTI/))
              CALL check_err(iret)
           if (i.eq.1) TIDE_DATA_ALL(:,:,I)=TIDE_DATA_ALL(:,:,I)*XCFAC+XCOFF
           if (i.eq.2) TIDE_DATA_ALL(:,:,I)=TIDE_DATA_ALL(:,:,I)*YCFAC+YCOFF
              END DO
          ELSE 
            DO I=1,NFIELDS
              iret=nf90_get_var(ncid,varidf(I),TIDE_DATA_ALL(:,:,I), &
                                start=(/1+(IAPROC-1)*NXLM,IY,1/),count=(/NXL(IAPROC),1,NTI/))
              CALL check_err(iret)
              if (i.eq.1) TIDE_DATA_ALL(:,:,I)=TIDE_DATA_ALL(:,:,I)*XCFAC+XCOFF
              if (i.eq.2) TIDE_DATA_ALL(:,:,I)=TIDE_DATA_ALL(:,:,I)*YCFAC+YCOFF
              END DO
          END IF
!
          DO JX=1,NXL(IAPROC)  ! This loop is parallelized with MPI 
!/DIST      IX = JX + (IAPROC-1)*NXLM ! IAPROC + (JX-1)*NAPROC
!/SHRD      IX = JX
!
!/TIDE            TIDE_NTI=0
!/TIDE            DO I=1,NTI
!
! Defines usable timesteps ... criteria could be improved
! remove the times when the point IX,IY is dry ... 
! and redefine TIDE_NTI based on wet times only
!
!/TIDE              IF (TIDE_DATA_ALL(JX,I,1).NE.fillvalue     &
!/TIDE                  .AND.TIDE_DATA_ALL(JX,I,NFIELDS).NE.fillvalue                   &
!/TIDE                  .AND.TIDE_DATA_ALL(JX,I,1).NE.0.0) THEN 
!/TIDE                TIDE_NTI=TIDE_NTI+1            
!/TIDE                TIDE_DATA(TIDE_NTI,:)=TIDE_DATA_ALL(JX,I,:)
!/TIDE                TIDE_DAYS(TIDE_NTI)=INT(ALLTIMES(I))
!/TIDE                TIDE_SECS(TIDE_NTI)=(ALLTIMES(I)-TIDE_DAYS(TIDE_NTI))*86400
!/TIDE              END IF
!/TIDE            END DO
!
!/TIDE            TIDE_HOURS(1:TIDE_NTI)=24.d0*dfloat(TIDE_DAYS(1:TIDE_NTI)) &
!/TIDE                                   +dfloat(TIDE_SECS(1:TIDE_NTI))/3600.d0
!
!/TIDE            IF (TIDE_NTI.GT.(TIDE_MF*3)) THEN 
!/TIDE              TIDE_LAT= YGRD(IY,IX)
!/TIDE              IF (ABS(TIDE_LAT).LT.5.) TIDE_LAT=SIGN(5.,TIDE_LAT)
!/TIDE                DO I=1,TIDE_NTI
!/TIDE                  CALL SETVUF(TIDE_HOURS(I),TIDE_LAT,I)
!/TIDE                  END DO
!/TIDE              CALL flex_tidana_webpage(IX,IY,XGRD(IY,IX),TIDE_LAT,TIDE_DAYS(1),TIDE_DAYS(TIDE_NTI), &
!/TIDE                                   TIDE_NDEF, TIDE_ITREND, RES, SSQ, RMSR0,          &
!/TIDE                                   SDEV0, RMSR, RESMAX, IMAX, 0)
                
!/TIDE              TIDAL_CONST(IX,IY,1:TIDE_MF,1:NFIELDS,1)=TIDE_AMPC(1:TIDE_MF,1:NFIELDS)
!/TIDE              TIDAL_CONST(IX,IY,1:TIDE_MF,1:NFIELDS,2)=TIDE_PHG(1:TIDE_MF,1:NFIELDS)
!/TIDET    WRITE (989,'(2I10,X,176F10.3)'),IX,TIDE_NTI,TIDE_AMPC(1:TIDE_MF,1:NFIELDS) !,TIDE_PHG(1:TIDE_MF,1:NFIELDS)
!/TIDET    WRITE (989,'(2I10,X,176F10.3)'),IX,TIDE_NTI,TIDE_PHG(1:TIDE_MF,1:NFIELDS)
!/TIDET    RPOS = 1_8 + LRECL*(IX-1_8)
!/TIDET    WRITE (990,POS=RPOS),NULLBUFF(1:NREC)
!/TIDET    WRITE (990,POS=RPOS),TIDE_AMPC(1:TIDE_MF,1:NFIELDS),TIDE_PHG(1:TIDE_MF,1:NFIELDS)
!              
! Only keeps constituents with high enough SNR ? (T-TEST > 2 : not a very good test ...)
!
!              DO J=1,TIDE_MF
!                DO I=1,NFIELDS
!                  WRITE(990+IAPROC,*) IX,IY,J,I,TIDE_TTEST(J,I),TIDE_AMPC(J,I),TIDE_PHG(J,I)
!                  IF (TIDE_TTEST(J,I).GT.1.) THEN 
!                    TIDAL_CONST(IX,IY,J,I,1)=TIDE_AMPC(J,I)
!                    TIDAL_CONST(IX,IY,J,I,2)=TIDE_PHG(J,I)
!                    END IF
!                  END DO
!                END DO

! Fills 1-D array for MPI communication
            IF (DISTSWITCH) THEN 
!/TIDE         DO J=1,TIDE_MF
!/TIDE           DO K=1,NFIELDS
!/TIDE             DO L=1,2
!/TIDE              IND=IND+1
!/TIDE               TIDALSTUFFA(IND)=TIDAL_CONST(IX,IY,J,K,L)
!/TIDE             END DO
!/TIDE           END DO
!/TIDE         END DO
!/TIDE         IF (IAPROC.EQ.NAPOUT) WRITE(6,'(A,I6,A,I6,A,I6)') 'IY, JX = ', &
!/TIDE           IY,',',JX, ' out of  ', NXL(IAPROC)
              END IF
!/TIDE            ELSE 
!/TIDE              WRITE(NDSE,*) 'WARNING NOT ENOUGH DATA AT POINT:',IX,IY, & 
!/TIDE                                                            NTI, TIDE_NTI
!/TIDE              TIDAL_CONST(IX,IY,1:TIDE_MF,1:NFIELDS,1)=0.
!/TIDE              TIDAL_CONST(IX,IY,1:TIDE_MF,1:NFIELDS,2)=0.
!/TIDE              IND=IND+TIDE_MF*NFIELDS*2
!/TIDE            END IF ! end of test on TIDE_NTI
!/TIDET            DO J=1,TIDE_MF
!/TIDET              DO K=1,NFIELDS
!/TIDET                DO L=1,2
!/TIDET                  WRITE(990+IAPROC,'(6I5,F10.4)') IAPROC,NAPOUT,IX,J,K,L,TIDAL_CONST(IX,1,J,K,L)
!/TIDET                END DO
!/TIDET              END DO
!/TIDET            END DO

          END DO  ! JX=1,NXL  - End of parallization
!/TIDET         CLOSE(990+IAPROC)

!
! We need to send and receive the bits of TIDAL_CONST(:,IY,:,:,:)
!
!/TIDE     IF (DISTSWITCH) THEN 
!/TIDE     IF (IAPROC.NE.NAPOUT)  THEN
!/DIST       CALL MPI_SEND ( TIDALSTUFFA, 1, TIDALSTUFF, NAPOUT-1, 0,  &
!/DIST                       MPI_COMM_WORLD, IERR_MPI )
!/TIDE     ELSE
!/TIDET           DO JX=1,NXL(IAPROC)
!/TIDET             IX = JX + (IAPROC-1)*NXLM ! IAPROC + (JX-1)*NAPROC
!/TIDET             DO J=1,TIDE_MF
!/TIDET               DO K=1,NFIELDS
!/TIDET                 DO L=1,2
!/TIDET    WRITE(987,'(6I5,F10.4)') IAPROC,IAPROC,IX,J,K,L,TIDAL_CONST(IX,1,J,K,L)
!/TIDET                 END DO
!/TIDET               END DO
!/TIDET             END DO
!/TIDET           END DO       ! JX=1,NXL
!/TIDE       DO JPC=1,NAPROC
!/TIDE         IF (JPC .NE. NAPOUT )  THEN 
!/DIST           CALL MPI_RECV ( TIDALSTUFFA, 1, TIDALSTUFF, JPC-1, MPI_ANY_TAG,  &
!/DIST                             MPI_COMM_WORLD, STAT_MPI, IERR_MPI )
!/TIDE           ! Now unpacks the 1-D array
!/TIDE           IND=0
!/TIDE           DO JX=1,NXL(JPC)
!/TIDE             IX = JX + (JPC-1)*NXLM ! JPC    + (JX-1)*NAPROC
!/TIDE!
!/TIDE             DO J=1,TIDE_MF
!/TIDE               DO K=1,NFIELDS
!/TIDE                 DO L=1,2
!/TIDE                   IND=IND+1
!/TIDE                   TIDAL_CONST(IX,IY,J,K,L)=TIDALSTUFFA(IND)
!/TIDE                 END DO
!/TIDE               END DO
!/TIDE             END DO
!/TIDE           END DO       ! JX=1,NXL
!/TIDE         END IF         ! JPC .NE. NAPOUT
!/TIDE       END DO           ! JPC=0,NAPROC-1
!/TIDET       WRITE(NDSE,*) "IY = ", IY, "on ", NY
!/TIDE     END IF
!/TIDE     END IF
          
        END DO ! IY=1,NY
!/TIDET    CLOSE (990)
!/TIDET    CLOSE (989)
!/TIDET  IF (IDFLD.EQ.'CUR') WRITE(986,'(F10.3,/)') TIDAL_CONST(:,1,15,1,1) 
!/TIDET  IF (IDFLD.EQ.'CUR') WRITE(986,'(F10.3,/)') TIDAL_CONST(:,1,15,1,1)

!/TIDE     IF (DISTSWITCH) THEN 
!/DIST       CALL MPI_TYPE_FREE(TIDALSTUFF, IERR_MPI)
!/TIDE       DEALLOCATE(TIDALSTUFFA)
!/TIDE       END IF

!/MPI        IF (IAPROC .EQ. NAPOUT )  THEN 
!/MPI          WRITE(NDSE,*) "parallelization done"
!/MPI        ELSE
!/MPI          GOTO 888
!/MPI        END IF
!/TIDE!
!/TIDE! After loop on points, write tidal constituents to file. 
!/TIDE!
!/TIDE        IF ( IAPROC .EQ. NAPOUT.AND.TIDEFLAG.GE.1)  &
!/TIDE          CALL W3FLDTIDE1 ( 'WRITE',  NDSDAT, NDST, NDSE, NX, NY, IDFLD, IERR )
!/TIDE          CALL W3FLDTIDE2 ( 'WRITE',  NDSDAT, NDST, NDSE, NX, NY, IDFLD, 0, IERR )
!/TIDE!
!/TIDE        GOTO 880
           END IF ! end of test   IF (ITYPE.GE.6.AND.TIDESWITCH)

!
!--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 6   Begin loop over input fields
!
!     Read scale factor and offset for input fields
        XCFAC = 1.0
        YCFAC = 1.0
        XCOFF = 0.0
        YCOFF = 0.0
!
        IF ( ITYPE .LE. 4 .OR. ITYPE.EQ.6 ) THEN
          IRET = NF90_INQUIRE_ATTRIBUTE(NCID,VARIDF(1),'scale_factor')
          IF ( IRET .EQ. NF90_NOERR ) THEN
             IRET = NF90_GET_ATT(NCID,VARIDF(1),'scale_factor',XCFAC)
          ENDIF
          IRET = NF90_INQUIRE_ATTRIBUTE(NCID,VARIDF(1),'add_offset')
          IF ( IRET .EQ. NF90_NOERR) THEN
             IRET = NF90_GET_ATT(NCID,VARIDF(1),'add_offset',XCOFF)
          END IF
          IF ( NFCOMP.EQ.2 .OR. IFLD.GE.3 .OR. FLBERG ) THEN
            IRET = NF90_INQUIRE_ATTRIBUTE(NCID,VARIDF(2),'scale_factor')
            IF (IRET .EQ. NF90_NOERR) THEN
               IRET = NF90_GET_ATT(NCID,VARIDF(2),'scale_factor',YCFAC)
            ENDIF
            IRET = NF90_INQUIRE_ATTRIBUTE(NCID,VARIDF(2),'add_offset')
            IF (IRET .EQ. NF90_NOERR) THEN
               IRET = NF90_GET_ATT(NCID,VARIDF(2),'add_offset',YCOFF)
            END IF
          END IF
        END IF
!
!/O15      J      = LEN_TRIM(FNMPRE)
!/O15      OPEN (NDSTIME,FILE=FNMPRE(:J)//'times.'//IDFLD,      &
!/O15            ERR=870,IOSTAT=IERR )
!
        IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,972)
        TIMESHIFT = 0
        DO ITIME=1,NTI
!
! 6.a Read new time and fields
!
          iret=nf90_inq_varid(ncid,"time",varidtmp)
          IF ( iret/=nf90_noerr ) iret=nf90_inq_varid(ncid,"MT",varidtmp)
          CALL check_err(iret)
          iret=nf90_get_var(ncid,varidtmp,datecur,start=(/ITIME/))
          call check_err(iret)
          IF (index(timeunits, "seconds").NE.0)   datecur=datecur/86400.
          IF (index(timeunits, "hours").NE.0)   datecur=datecur/24. 
          datecur=dateref+datecur
!
          CALL caldat(INT(datecur),dcur,mcur,ycur)

          hcur=(datecur-INT(datecur))*24.
          IF (ABS(hcur-NINT(hcur)).LT.1.E-5) hcur=NINT(hcur)
          mncur=(hcur-INT(hcur))*60.
          IF (ABS(mncur-NINT(mncur)).LT.1.E-5) mncur=NINT(mncur)
          scur=(mncur-INT(mncur))*60.
          IF (ABS(scur-NINT(scur)).LT.1.E-5) scur=NINT(scur)
! Conversion to WW3 time integers
          TIME(1)=INT(ycur)*10000+INT(mcur)*100 +INT(dcur) 
          TIME(2)=INT(hcur)*10000+INT(mncur)*100+INT(scur) 
          IF ( .NOT.FLTIME.AND.ITIME.EQ.1) THEN
            TIMESHIFT = DSEC21 ( TIME ,TIMESTART)
            END IF
! Shifts time to start on TIMESTART if FLTIME is false
        CALL TICK21 ( TIME, TIMESHIFT )
!
        CALL STME21 ( TIME , IDTIME )
        IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,973) IDTIME
!/O15        WRITE (NDSTIME, 979, ERR=871,IOSTAT=IERR) TIME
!/O3        IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,974)
!
! ... Input
!
        IF ( ITYPE .LE. 4 .OR. ITYPE.EQ.6) THEN
            IF (NDIMSGRID.EQ.1) THEN 
              iret=nf90_get_var(ncid,varidf(1),XC(:,1),start=(/1,ITIME/),count=(/MXM,1/))
            ELSE
              IF (NDIMSVAR.EQ.3) THEN 
                iret=nf90_get_var(ncid,varidf(1),XC,start=(/1,1,ITIME/),count=(/MXM,MYM,1/))
              ELSE 
                iret=nf90_get_var(ncid,varidf(1),XC,start=(/1,1,1,ITIME/),count=(/MXM,MYM,1,1/))
                END IF
              END IF
! The following line forces to 0 values that are undefine
            WHERE(XC.NE.XC) XC = 0. 
            WHERE(XC.EQ.fillvalue) XC = 0.
            CALL check_err(iret)
            XC  = XC * XCFAC + XCOFF
!
!/T2            WRITE (NDST,9060) 1
!/T2            IXP0   = 1
!/T2            IXPN   = MIN ( IXP0+IXPWDT-1 , NXJ(1) )
!/T2            DO
!/T2              CALL PRTBLK ( NDST, NXJ(1), NYJ(1), MXM, XC, MASK, 0, 0.,&
!/T2                            IXP0, IXPN, 1, 1, NYJ(1), 1, 'Field 1', ' ')
!/T2              IF (IXPN.NE.NXJ(1)) THEN
!/T2                  IXP0    = IXP0 + IXPWDT
!/T2                  IXPN    = MIN ( IXPN+IXPWDT , NXJ(1) )
!/T2                ELSE
!/T2                  EXIT
!/T2                END IF
!/T2              END DO
!
            IF (NFCOMP.EQ.2 .OR. IFLD.GE.3 .OR. FLBERG) THEN

! This is a quick fix that works if the lon,lat,level,time dimensions are in that order
! otherwise, one should check the length of each dimension ... 
                    IF (NDIMSGRID.EQ.1) THEN
                      iret=nf90_get_var(ncid,varidf(2),YC(:,1),start=(/1,ITIME/),count=(/MXM,1/))
                    ELSE
                      IF (NDIMSVAR.EQ.3) THEN
                        iret=nf90_get_var(ncid,varidf(2),YC,start=(/1,1,ITIME/),count=(/MXM,MYM,1/))

                      ELSE
                        iret=nf90_get_var(ncid,varidf(2),YC,start=(/1,1,1,ITIME/),count=(/MXM,MYM,1,1/))
                        END IF
                      END IF
! The following line forces to 0 values that are undefine
                    WHERE(YC.NE.YC) YC = 0.
                    WHERE(YC.EQ.fillvalue) YC = 0.
                    CALL check_err(iret)
                    YC  = YC * YCFAC + YCOFF
!
!/T2                WRITE (NDST,9060) 2
!/T2                IXP0   = 1
!/T2                IXPN   = MIN ( IXP0+IXPWDT-1 , NXJ(2) )
!/T2                DO
!/T2                  CALL PRTBLK ( NDST, NXJ(2), NYJ(2), MXM, YC, MASK, 0, 0., &
!/T2                                IXP0, IXPN, 1, 1, NYJ(2), 1, 'Field 2', ' ')
!/T2                  IF (IXPN.NE.NXJ(2)) THEN
!/T2                      IXP0    = IXP0 + IXPWDT
!/T2                      IXPN    = MIN ( IXPN+IXPWDT , NXJ(2) )
!/T2                    ELSE
!/T2                      EXIT
!/T2                    END IF
!/T2                  END DO
!
                IF (FLSTAB) THEN
! This is a quick fix that works if the lon,lat,level,time dimensions are in that order
! otherwise, one should check the length of each dimension ... 
                    IF (NDIMSGRID.EQ.1) THEN 
                      iret=nf90_get_var(ncid,varidf(3),AC(:,1),start=(/1,ITIME/),count=(/MXM,1/))
                    ELSE 
                      IF (NDIMSVAR.EQ.3) THEN
                        iret=nf90_get_var(ncid,varidf(3),AC,start=(/1,1,ITIME/),count=(/MXM,MYM,1/))
                      ELSE
                        iret=nf90_get_var(ncid,varidf(3),AC,start=(/1,1,1,ITIME/),count=(/MXM,MYM,1,1/))
                        END IF
                      END IF
                    CALL check_err(iret)
                    !AC(:,:)=AC(:,MYM:1:-1)
!
!/T2                    WRITE (NDST,9060) 3
!/T2                    IXP0   = 1
!/T2                    IXPN   = MIN ( IXP0+IXPWDT-1 , NXJ(2) )
!/T2                    DO
!/T2                      CALL PRTBLK ( NDST, NXJ(2), NYJ(2), MXM, AC, MASK, 0,&
!/T2                            0., IXP0, IXPN, 1,1, NYJ(2), 1, 'Field 3', ' ')
!/T2                      IF (IXPN.NE.NXJ(2)) THEN
!/T2                          IXP0    = IXP0 + IXPWDT
!/T2                          IXPN    = MIN ( IXPN+IXPWDT , NXJ(2) )
!/T2                        ELSE
!/T2                          EXIT
!/T2                        END IF
!/T2                      END DO
!
                  END IF
!
              END IF

          ELSE  ! ITYPE .NE. 5
!
            IF ( IAPROC .EQ. NAPOUT ) WRITE(NDSO,*) "ITYPE5 TO DO"
            IF (IDFMF(1).EQ.3) THEN
                READ (NDSF(1),  END=862,ERR=862,IOSTAT=IERR) NDAT
              ELSE
                READ (NDSF(1),*,END=862,ERR=862,IOSTAT=IERR) NDAT
              END IF
!/O3           IF ( IAPROC .EQ. NAPOUT )  WRITE (NDSO,975) NDAT
            IF ( NDAT.GT.0 ) THEN
                ALLOCATE ( DATA(RECLDT,NDAT) )
                DO IDAT=1, NDAT
                  IF (IDFMF(1).EQ.1) THEN
                      READ (NDSF(1),    *   ,END=863,ERR=863,         &
                            IOSTAT=IERR) DATA(:,IDAT)
                    ELSE IF (IDFMF(1).EQ.2) THEN
                      READ (NDSF(1),FORMT(1),END=863,ERR=863,         &
                            IOSTAT=IERR) DATA(:,IDAT)
                    ELSE
                      READ (NDSF(1),         END=863,ERR=863,         &
                            IOSTAT=IERR) DATA(:,IDAT)
                    END IF
                  END DO
              END IF
!
!/T2            WRITE (NDST,9061)
!/T2            DO IDAT=1, NDAT
!/T2              IX     = MIN(6,RECLDT)
!/T2              WRITE (NDST,9062) IDAT, DATA(1:IX,IDAT)
!/T2              IF ( IX.LT.RECLDT ) WRITE (NDST,9063) DATA(IX+1:,:)
!/T2              END DO
!
          END IF
!
! 6.b Interpolate fields
! ... No Interpolation, type AI (should not use array syntax !!!)
!
        IF (ITYPE.EQ.1.OR.ITYPE.EQ.6) THEN
!
            ! change fillvalue
            DO IY=1,NY
              DO IX=1,NX
                IF (XC(IX,IY) .EQ. fillvalue) XC(IX,IY)=0
                IF (YC(IX,IY) .EQ. fillvalue) YC(IX,IY)=0
              END DO
            END DO

            IF (( IFLD.LE.2 ).AND.( .NOT. FLBERG )) THEN
                DO IY=1, NY
                  DO IX=1, NX
                    FA(IX,IY) = XC(IX,IY)
                    END DO
                  END DO
              ELSE
                DO IY=1, NY
                  DO IX=1, NX
                    FX(IX,IY) = XC(IX,IY)
                    FY(IX,IY) = YC(IX,IY)
                    FA(IX,IY) = AC(IX,IY)
                    END DO
                  END DO
              END IF
!
          ELSE IF (ITYPE.NE.5) THEN

            ! change fillvalue
            ! this is done in a previous test
            !DO IY=1,NY
            !  DO IX=1,NX
            !    IF (XC(IX21(IX,IY),IY21(IX,IY)) .EQ. fillvalue) XC(IX21(IX,IY),IY21(IX,IY))=0
            !    IF (XC(IX22(IX,IY),IY21(IX,IY)) .EQ. fillvalue) XC(IX22(IX,IY),IY21(IX,IY))=0
            !    IF (XC(IX21(IX,IY),IY22(IX,IY)) .EQ. fillvalue) XC(IX21(IX,IY),IY22(IX,IY))=0
            !    IF (XC(IX22(IX,IY),IY22(IX,IY)) .EQ. fillvalue) XC(IX22(IX,IY),IY22(IX,IY))=0
            !    IF (YC(IX21(IX,IY),IY21(IX,IY)) .EQ. fillvalue) YC(IX21(IX,IY),IY21(IX,IY))=0
            !    IF (YC(IX22(IX,IY),IY21(IX,IY)) .EQ. fillvalue) YC(IX22(IX,IY),IY21(IX,IY))=0
            !    IF (YC(IX21(IX,IY),IY22(IX,IY)) .EQ. fillvalue) YC(IX21(IX,IY),IY22(IX,IY))=0
            !    IF (YC(IX22(IX,IY),IY22(IX,IY)) .EQ. fillvalue) YC(IX22(IX,IY),IY22(IX,IY))=0
            !  END DO
            !END DO
!
! ... One-component fields
!
!/O3           IF ( IAPROC .EQ. NAPOUT )  WRITE (NDSO,976) ' '
            IF (( IFLD.LE.2 ).AND.( .NOT. FLBERG )) THEN
!
                DO IY=1,NY
                  DO IX=1,NX
                    FA(IX,IY)                                         &
                          = RD11(IX,IY) * XC(IX21(IX,IY),IY21(IX,IY)) &
                          + RD21(IX,IY) * XC(IX22(IX,IY),IY21(IX,IY)) &
                          + RD12(IX,IY) * XC(IX21(IX,IY),IY22(IX,IY)) &
                          + RD22(IX,IY) * XC(IX22(IX,IY),IY22(IX,IY))
                    END DO
                  END DO
!
                IF (NFCOMP.EQ.2) THEN
!/O3             IF ( IAPROC .EQ. NAPOUT )    WRITE (NDSO,976) ' (2) '
                    DO IY=1,NY
                      DO IX=1,NX
                        FA(IX,IY) = FA(IX,IY)                         &
                          + XD11(IX,IY) * YC(JX21(IX,IY),JY21(IX,IY)) &
                          + XD21(IX,IY) * YC(JX22(IX,IY),JY21(IX,IY)) &
                          + XD12(IX,IY) * YC(JX21(IX,IY),JY22(IX,IY)) &
                          + XD22(IX,IY) * YC(JX22(IX,IY),JY22(IX,IY))
                        END DO
                      END DO
                  END IF
!
! ... Two-component fields
!
              ELSE  !so if IFLD.GT.2
!
                DO IY=1,NY
                  DO IX=1,NX
                    FX(IX,IY)                                         &
                          = RD11(IX,IY) * XC(IX21(IX,IY),IY21(IX,IY)) &
                          + RD21(IX,IY) * XC(IX22(IX,IY),IY21(IX,IY)) &
                          + RD12(IX,IY) * XC(IX21(IX,IY),IY22(IX,IY)) &
                          + RD22(IX,IY) * XC(IX22(IX,IY),IY22(IX,IY))
                    FY(IX,IY)                                         &
                          = RD11(IX,IY) * YC(IX21(IX,IY),IY21(IX,IY)) &
                          + RD21(IX,IY) * YC(IX22(IX,IY),IY21(IX,IY)) &
                          + RD12(IX,IY) * YC(IX21(IX,IY),IY22(IX,IY)) &
                          + RD22(IX,IY) * YC(IX22(IX,IY),IY22(IX,IY))
                    FA(IX,IY)                                         &
                          = RD11(IX,IY) * AC(IX21(IX,IY),IY21(IX,IY)) &
                          + RD21(IX,IY) * AC(IX22(IX,IY),IY21(IX,IY)) &
                          + RD12(IX,IY) * AC(IX21(IX,IY),IY22(IX,IY)) &
                          + RD22(IX,IY) * AC(IX22(IX,IY),IY22(IX,IY))
                    A1(IX,IY) = MAX ( 1.E-10 ,                        &
                                  SQRT( FX(IX,IY)**2 + FY(IX,IY)**2 ) )
                    A2(IX,IY)                                         &
                         = RD11(IX,IY) * SQRT(XC(IX21(IX,IY),IY21(IX,IY))**2  &
                                             +YC(IX21(IX,IY),IY21(IX,IY))**2) &
                         + RD21(IX,IY) * SQRT(XC(IX22(IX,IY),IY21(IX,IY))**2  &
                                             +YC(IX22(IX,IY),IY21(IX,IY))**2) &
                         + RD12(IX,IY) * SQRT(XC(IX21(IX,IY),IY22(IX,IY))**2  &
                                             +YC(IX21(IX,IY),IY22(IX,IY))**2) &
                         + RD22(IX,IY) * SQRT(XC(IX22(IX,IY),IY22(IX,IY))**2  &
                                             +YC(IX22(IX,IY),IY22(IX,IY))**2)
                    A3(IX,IY) = SQRT (                                &
                           RD11(IX,IY) * ( XC(IX21(IX,IY),IY21(IX,IY))**2   &
                                         + YC(IX21(IX,IY),IY21(IX,IY))**2 ) &
                         + RD21(IX,IY) * ( XC(IX22(IX,IY),IY21(IX,IY))**2   &
                                         + YC(IX22(IX,IY),IY21(IX,IY))**2 ) &
                         + RD12(IX,IY) * ( XC(IX21(IX,IY),IY22(IX,IY))**2   &
                                         + YC(IX21(IX,IY),IY22(IX,IY))**2 ) &
                         + RD22(IX,IY) * ( XC(IX22(IX,IY),IY22(IX,IY))**2   &
                                         + YC(IX22(IX,IY),IY22(IX,IY))**2 ) )
                    END DO
                  END DO
!
! ... Winds, correct for velocity or energy conservation
!
!/WNT1                IF (IFLD.EQ.3) THEN
!/WNT1                    DO IY=1,NY
!/WNT1                      DO IX=1,NX
!/WNT1                        FACTOR = MIN ( 1.5 , A2(IX,IY)/A1(IX,IY) )
!/WNT1                        FX(IX,IY) = FACTOR * FX(IX,IY)
!/WNT1                        FY(IX,IY) = FACTOR * FY(IX,IY)
!/WNT1                        END DO
!/WNT1                      END DO
!/WNT1                  END IF
!
!/WNT2                IF (IFLD.EQ.3) THEN
!/WNT2                    DO IY=1,NY
!/WNT2                      DO IX=1,NX
!/WNT2                        FACTOR = MIN ( 1.5 , A3(IX,IY)/A1(IX,IY) )
!/WNT2                        FX(IX,IY) = FACTOR * FX(IX,IY)
!/WNT2                        FY(IX,IY) = FACTOR * FY(IX,IY)
!/WNT2                        END DO
!/WNT2                      END DO
!/WNT2                  END IF
!
! ... Currents, correct for velocity or energy conservation
!
!/CRT1                IF (IFLD.EQ.4) THEN
!/CRT1                    DO IY=1,NY
!/CRT1                      DO IX=1,NX
!/CRT1                        FACTOR = MIN ( 1.5 , A2(IX,IY)/A1(IX,IY) )
!/CRT1                        FX(IX,IY) = FACTOR * FX(IX,IY)
!/CRT1                        FY(IX,IY) = FACTOR * FY(IX,IY)
!/CRT1                        END DO
!/CRT1                      END DO
!/CRT1                  END IF
!
!/CRT2                IF (IFLD.EQ.4) THEN
!/CRT2                    DO IY=1,NY
!/CRT2                      DO IX=1,NX
!/CRT2                        FACTOR = MIN ( 1.5 , A3(IX,IY)/A1(IX,IY) )
!/CRT2                        FX(IX,IY) = FACTOR * FX(IX,IY)
!/CRT2                        FY(IX,IY) = FACTOR * FY(IX,IY)
!/CRT2                        END DO
!/CRT2                      END DO
!/CRT2                 END IF
!
              END IF
!
          END IF
!
! ... Test output
!
!/T3        IF ( .NOT. ALLOCATED(MAPOUT) ) ALLOCATE ( MAPOUT(NX,NY) )
!/T3        WRITE (NDST,9065)
!/T3        DO IX=1, NX
!/T3          DO IY=1, NY
!/T3            MAPOUT(IX,IY) = MAPSTA(IY,IX)
!/T3            END DO
!/T3          END DO
!/T3        IX0    = 1
!/T3        IXN    = MIN ( IX0+IXWDT-1 , NX )
!/T3        DO
!/T3          IF (IFLD.EQ.1) THEN
!/T3              CALL PRTBLK (NDSO, NX, NY, NX, FA, MAPOUT, 0, 0.,   &
!/T3                      IX0, IXN, 1, 1, NY, 1, 'Fraction ice', '(-)')
!/T3              IF ( FLBERG )                                       &
!/T3              CALL PRTBLK (NDSO, NX, NY, NX, FA, MAPOUT, 0, 0.,   &
!/T3                      IX0, IXN, 1, 1, NY, 1, 'Iceberg a', '0.1/km')
!/T3            ELSE IF (IFLD.EQ.2) THEN
!/T3              CALL PRTBLK (NDSO, NX, NY, NX, FA, MAPOUT, 0, 0.,   &
!/T3                      IX0, IXN, 1, 1, NY, 1, 'Water level', 'm')
!/T3            ELSE
!/T3              CALL PRTBLK (NDSO, NX, NY, NX, FX, MAPOUT, 0, 0.,   &
!/T3                      IX0, IXN, 1, 1, NY, 1, 'Cart. X-comp', 'm/s')
!/T3              CALL PRTBLK (NDSO, NX, NY, NX, FY, MAPOUT, 0, 0.,   &
!/T3                      IX0, IXN, 1, 1, NY, 1, 'Cart. Y-comp', 'm/s')
!/T3              IF ( FLSTAB )                                       &
!/T3              CALL PRTBLK (NDSO, NX, NY, NX, FA, MAPOUT, 0, 0.,   &
!/T3                      IX0, IXN, 1, 1, NY, 1, 'Tair-Tsea', 'degr')
!/T3            END IF
!/T3          IF (IXN.NE.NX) THEN
!/T3              IX0    = IX0 + IXWDT
!/T3              IXN    = MIN ( IXN+IXWDT , NX )
!/T3            ELSE
!/T3              EXIT
!/T3            END IF
!/T3          END DO
!
! 6.c Write fields
!
        IF ( ITYPE .LE. 4 .OR. ITYPE.EQ.6 ) THEN
!/O3      IF ( IAPROC .EQ. NAPOUT )  WRITE (NDSO,977)
          IF ( IAPROC .EQ. NAPOUT ) CALL W3FLDG ('WRITE', IDFLD, NDSDAT, NDST, NDSE, NX, NY,  &
                       NX, NY, TIME, TIME, TIME, FX, FY, FA, TIME,  &
                       FX, FY, FA, IERR)

        ELSE IF ( ITYPE .EQ. 5 ) THEN
          IF ( NDAT .EQ. 0 ) THEN
!/O3        IF ( IAPROC .EQ. NAPOUT )  WRITE (NDSO,978)
          ELSE
!/O3        IF ( IAPROC .EQ. NAPOUT )  WRITE (NDSO,977)
            IF ( IAPROC .EQ. NAPOUT ) CALL W3FLDD ('WRITE', IDFLD, NDSDAT, NDST, NDSE, TIME,&
                         TIME, RECLDT, NDAT, IDAT, DATA, IERR )
            DEALLOCATE ( DATA )
          END IF
        END IF
        IF (IERR.NE.0) CALL EXTCDE ( 30 )
!
      END DO
!
      DEALLOCATE(XC,YC,AC)  
      IF (ASSOCIATED(ALA)) DEALLOCATE(ALA,ALO)
!
!     End loop over input fields
!--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
  880 CONTINUE
      GOTO 888
!
! Error escape locations
!
  800 CONTINUE
      WRITE (NDSE,1000) IERR
      CALL EXTCDE ( 40 )
!
  801 CONTINUE
      WRITE (NDSE,1001)
      CALL EXTCDE ( 41 )
!
  802 CONTINUE
      WRITE (NDSE,1002) IERR
      CALL EXTCDE ( 42 )
!
  803 CONTINUE
      WRITE (NDSE,1003) IERR
      CALL EXTCDE ( 43 )
!
  804 CONTINUE
      WRITE (NDSE,1004)
      CALL EXTCDE ( 44 )
!
  805 CONTINUE
      WRITE (NDSE,1005) IERR
      CALL EXTCDE ( 45 )
!
  837 CONTINUE
      WRITE (NDSE,1037) NDSI
      CALL EXTCDE ( 46 )
!
  840 CONTINUE
      WRITE (NDSE,1040) X0I, XNI
      CALL EXTCDE ( 47 )
!
  841 CONTINUE
      WRITE (NDSE,1041) Y0I, YNI
      CALL EXTCDE ( 48 )
!
  845 CONTINUE
      WRITE (NDSE,1045) IERR
      CALL EXTCDE ( 49 )
!
  846 CONTINUE
      WRITE (NDSE,1046) IERR
      CALL EXTCDE ( 50 )
!
  850 CONTINUE
      WRITE (NDSE,1050) IERR, NDSF(J), NAMEF
      CALL EXTCDE ( 51 )
!
  860 CONTINUE
      WRITE (NDSE,1060) J, IERR
      CALL EXTCDE ( 52 )
!
  861 CONTINUE
      WRITE (NDSE,1061) TIME, TIME2
      CALL EXTCDE ( 53 )
!
  862 CONTINUE
      WRITE (NDSE,1062) IERR
      CALL EXTCDE ( 54 )
!
  863 CONTINUE
      WRITE (NDSE,1063) IDAT, IERR
      CALL EXTCDE ( 55 )
  864 CONTINUE
      WRITE (NDSE,1064)  TRIM(STRDIMSNAME)
      CALL EXTCDE ( 56 )
!
!/O15  870 CONTINUE
!/O15      WRITE (NDSE,1070) IDFLD, IERR
!/O15      CALL EXTCDE ( 57 )
!
!/O15  871 CONTINUE
!/O15      WRITE (NDSE,1071) IDTIME, IERR
!/O15      CALL EXTCDE ( 58 )
!
  888 CONTINUE
      IF ( IAPROC .EQ. NAPOUT ) WRITE (NDSO,999)
!/MPI      CALL MPI_FINALIZE  ( IERR_MPI )

!
!/NCO/!     CALL W3TAGE('WAVEPREP')


!
! Formats
!
  900 FORMAT (/15X,'  *** WAVEWATCH III  Input pre-processing ***  '/ &
               15X,'==============================================='/)
  901 FORMAT ( '  Comment character is ''',A,''''/)
  902 FORMAT ( '  Grid name : ',A/)
!
  930 FORMAT (/'  Description of inputs'/                             &
               ' --------------------------------------------------'/ &
               '       Input type        : ',A/                       &
               '       Format type       : ',A)
 1930 FORMAT ( '          Field conserves velocity.')
 2930 FORMAT ( '          Field corrected for energy conservation.')
  931 FORMAT (/'       Time shifted. Starting: ',A)
  932 FORMAT (/'       Input grid dim.   :',I9,3X,I5)
  933 FORMAT ( '       Longitude range   :',2F8.2,' (deg)'/           &
               '       Latitude range    :',2F8.2,' (deg)')
  733 FORMAT ( '       X range           :',2F8.2,' (km)'/           &
               '       Y range           :',2F8.2,' (km)')
  934 FORMAT (/'       Data type         : ',A/                       &
               '       Data record length:',I5/                       &
               '       Missing values    :',F8.2)
  935 FORMAT ( 'DT',I1 )
  938 FORMAT ( '          Icebergs included.')
  939 FORMAT ( '          Air-sea temperature differences included.')
!
  940 FORMAT (//'  Preprocessing data'/                               &
               ' --------------------------------------------------')
  941 FORMAT ( '       Interpolation factors ..... '/                 &
               '          (longitude-latitude grid)')
  942 FORMAT ( '       Interpolation factors ..... '/                 &
               '          (grid from file)')
  943 FORMAT (/'           Longitude-latitude file ',I1,' :'/         &
               '          ---------------------------------------')
  944 FORMAT ( '           Input grid dim.   :',I9,3X,I5/             &
               '           Closed longitudes :',L5)
  945 FORMAT ( '           Layout indicator  :',I5/                   &
               '           Format indicator  :',I5)
  946 FORMAT ( '           Format            : ',A)
  947 FORMAT ( '           Unit number       :',I5)
  948 FORMAT ( '           File name         : ',A)
  949 FORMAT (/'           Corresponding map file '/                  &
               '          ---------------------------------------')
!
  960 FORMAT (/'           Data file :'/                              &
               '          ---------------------------------------')
  961 FORMAT (/'           Data file :'/                              &
               '          ---------------------------------------'/   &
               '           Input grid dim.   :',I9,3X,I5)
  962 FORMAT (/'           Data file (',I1,') :'/                     &
               '          ---------------------------------------'/   &
               '           Input grid dim.   :',I9,3X,I5)
  963 FORMAT ( '           Layout indicator  :',I5)
  964 FORMAT ( '           Format indicator  :',I5)
  965 FORMAT ( '           Format for time   : ',A/                   &
               '           Format for data   : ',A)
  966 FORMAT ( '           Unit number       :',I5)
  967 FORMAT ( '           File name         : ',A)
!
  970 FORMAT (/'       Opening input data file .....')
  971 FORMAT (/'       Opening output data file .....')
  972 FORMAT (//'  Processing data'/                                  &
               ' --------------------------------------------------')
  973 FORMAT ( '           Time : ',A)
!/O3  974 FORMAT ( '                  reading ....')
!/O3  975 FORMAT ( '                     number of data records :',I6)
!/O3  976 FORMAT ( '                  interpolating',A,'....')
!/O3  977 FORMAT ( '                  writing ....')
!/O3  978 FORMAT ( '                  skipping ....')
!
!/O15 979 FORMAT (1X,I8.8,1X,I6.6)
!
  999 FORMAT(//'  End of program '/                                   &
               ' ========================================='/          &
               '         WAVEWATCH III Input preprocessing '/)
!
 1000 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ERROR IN OPENING INPUT FILE'/                    &
               '     IOSTAT =',I5/)
!
 1001 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     PREMATURE END OF INPUT FILE'/)
!
 1002 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ERROR IN READING FROM INPUT FILE'/               &
               '     IOSTAT =',I5/)
!
 1003 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ERROR IN READING FROM INPUT FILE'/               &
               '     EXPECTING LIST OF TIDAL CONST. OR FAST OR VFAST'/&
               '     IOSTAT =',I5/)
!
 1004 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     PREMATURE END OF TIME ATTRIBUTE '/               &
               '     FROM INPUT DATA FILE'/)
!
 1005 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ERROR IN READING OF TIME ATTRIBUTE '/            &
               '     FROM INPUT DATA FILE'/                           &
               '     IOSTAT =',I5/)
!
 1028 FORMAT (/' *** WAVEWATCH III WARNING IN W3PRNC : '/             &
               '     CALENDAR ATTRIBUTE NOT FOUND'/                   &
               '     IT MUST RESPECT JULIAN OR GREGORIAN CALENDAR')
 1029 FORMAT (/' *** WAVEWATCH III WARNING IN W3PRNC : '/             &
               '     CALENDAR ATTRIBUTE NOT MATCH'/                   &
               '     IT MUST RESPECT JULIAN OR GREGORIAN CALENDAR')
 1030 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ILLEGAL FIELD ID -->',A,'<--'/)
 1031 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ILLEGAL FORMAT ID -->',A,'<--'/)
 1032 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     LATITUDE VALUES MUST BE REVERSED'/               &
               '     EXAMPLE: ncpdq -h -O -a -lat file.nc'/ )
!
 1033 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ILLEGAL DATA RECORD LENGTH : ',I6/)
 1034 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ILLEGAL DATA TYPE : ',I2/)
!
 1035 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ILLEGAL TIME : ',I8.8,I7.6/)
 1036 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ILLEGAL SIZE OF INPUT GRID : ',I5,1X,I5/)
 1037 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ATTEMPT TO USE INPUT UNIT (',I2,                 &
               ') FOR NAMED FILE'/)
 1038 FORMAT (/' *** WAVEWATCH III WARNING IN W3PRNC : '/             &
               '     DATA READ FROM INPUT FILE')
 1039 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ARRAYS FOR INPUT FIELDS TO SMALL '/              &
               '            GRID SIZE       :',2I5/                   &
               '            ARRAY DIMENSION :',2I5/)
!
 1040 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     LONGITUDE GRID RANGE > 360 DEGR. : ',2F6.1/)
 1041 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ILLEGAL LATITUDE GRID RANGE : ',2F6.1/)
 1042 FORMAT (/' *** WAVEWATCH-III WARNING W3PRNC : '/                &
               '     GRID POINT ',2I6,2F7.2,/                         &
               ' NOT COVERED BY INPUT GRID.'/)
 1043 FORMAT (/' *** WAVEWATCH III WARNING W3PRNC : '/                &
               '     X = ',F10.1,' NOT COVERED BY INPUT GRID.'/)
 1044 FORMAT (/' *** WAVEWATCH III WARNING W3PRNC : '/                &
               '     Y = ',F10.1,' NOT COVERED BY INPUT GRID.'/)
!

!
 1045 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ERROR IN OPENING LAT-LONG DATA FILE'/            &
               '     IOSTAT =',I5/)
!
 1046 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ERROR IN OPENING MASK FILE'/                     &
               '     IOSTAT =',I5/)
!
 1050 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ERROR IN OPENING INPUT DATA FILE'/               &
               '     IOSTAT =',I5/                                    &
               '     NDSF   =',I5/                                    &
               '     NAMEF  = ',A/)
 1051 FORMAT (/' *** WAVEWATCH III ERROR IN W3PREP : '/               &
               '     CANNOT READ UNFORMATTED FROM UNIT',I3/)
!
 1060 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ERROR IN READING TIME FROM FILE (',I1,')'/       &
               '     IOSTAT =',I5/)
 1061 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     INCOMPATIBLE FIELD TIMES '/                      &
               '            FIELD #1 : ',I8.8,I7.6/                   &
               '            FIELD #2 : ',I8.8,I7.6/)
 1062 FORMAT (/' *** WAVEWATCH III ERROR IN W3PREP : '/               &
               '     ERROR IN READING NDAT FROM FILE'/                &
               '     IOSTAT =',I5/)
 1063 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     ERROR IN READING DATA RECORD',I6,' FROM FILE'/   &
               '     IOSTAT =',I5/)
 1064 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/               &
               '     GRID DIMENSIONS ', A,' NOT FOUND... CHECK DIMENSION NAMES')
!
!/O15 1070 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/          &
!/O15               '     ERROR IN CREATING A TIMES FILE FOR ',A/     &
!/O15               '     IOSTAT =',I5/)
!/O15 1071 FORMAT (/' *** WAVEWATCH III ERROR IN W3PRNC : '/          &
!/O15               '     ERROR IN WRITING TIME OUTPUT ',A/           &
!/O15               '     IOSTAT =',I5/)
!
!/T 9000 FORMAT (' TEST W3PRNC : ACC : ',F6.3)
!
!/T 9040 FORMAT (' TEST W3PRNC : INPUT GRID RANGES AND INCR. AFTER CORR.'/ &
!/T              '               LON / X : ',3F10.2,                   &
!/T              ' (GLOBAL=',L1,')'/                                     &
!/T              '               LAT / Y : ',3F10.2)
!/T 9041 FORMAT (' TEST W3PRNC : INTERPOLATION DATA FOR ',A)
!/T 9042 FORMAT ('              ',I4,F8.2,2I4,2F8.2,1X,F6.3,1X,A)
!/T 9043 FORMAT (' TEST W3PRNC : GRID SHIFTED BY ',F5.0,' DEGREES / M')
!/T1 9045 FORMAT (' TEST W3PRNC : IX, IY, IXI(2), IYI(2), RD(4)')
!/T1 9046 FORMAT ('     ',2I4,2X,4I4,2X,4F6.2)
!
!/T1a 9050 FORMAT (' TEST W3PRNC : LAT-LONG OF INPUT FILE ')
!/T1a 9051 FORMAT ('             ',2I4,2F8.2,I4)
!
!/T2 9060 FORMAT (' TEST W3PRNC : INPUT FIELD (',I1,') :'/)
!/T2 9061 FORMAT (' TEST W3PRNC : INPUT DATA RECORDS :')
!/T2 9062 FORMAT ('      ',I6,' : ',6E11.3)
!/T2 9063 FORMAT ('               ',6E11.3)
!/T3 9065 FORMAT (' TEST W3PRNC : OUTPUT FIELD(S) :'/)
!/
!/ End of W3PRNC ----------------------------------------------------- /
!/

END PROGRAM W3PRNC

!*********************************************************************

  subroutine check_err(iret)

  USE NETCDF

  implicit none

  integer iret

  if (iret .ne. NF90_NOERR) then
     print *, nf90_strerror(iret)
     stop
  endif
  return
  END subroutine check_err


