SUBROUTINE getflt(mode,iunit,itime,ident,ldata,fdata,igtype,gparam,nfskip,ierr)

!// 16 bit input and unpack (not handling 'undefined')

!// input:
!//    mode:   0 = read field
!//            1 = read field, skip fields until time > itime
!//            2 = read field if field time = itime
!//                (otherwise next read starts at the same field)
!//          100 = read field identification
!//                (next read starts at the same field)
!//          101 = read field identification, skip fields
!//                until time > itime
!//                (next read starts at the same (last) field)
!//          102 = read field identification, skip fields
!//                until time >= itime
!//                (next read starts at the same (last) field)
!//          200 = scan rest of the file and read field with
!//                matching identification, specified identification
!//                input in ident(1:20) where -32767 means any value
!//          201 = scan the whole file and read field with
!//                matching identification, specified identification
!//                input in ident(1:20) where -32767 means any value
!//           -1 = clean up after a file is closed, and the same
!//                file unit no. is used for another file.
!//    iunit: file unit no.
!//    itime: used if mode is 1, 2, 101 or 102.
!//           itime(1) - year
!//                (2) - month
!//                (3) - day
!//                (4) - hour (utc)
!//                (5) - forecast time in hours
!//                      (added to previous date and time)
!//    ident:  used if mode is 200 or 201
!//    ldata:  length of fdata (max field size)

!// output:
!//    ident(20): field identification
!//    fdata(..): field (unscaled, according to identification)
!//    igtype:    grid type (1=polar. 2=geo. 3+rot.sph. ...)
!//    gparam(6): grid parameters
!//    nfskip:    no. of fields skipped (for mode=2,101,102)
!//    ierr = 0:  read o.k.
!//           1:  read error
!//           2:  read error, end_of_file
!//           3:  no field returned due to mode=2 and itime spec.

!// note: mode=100,101,102 and possibly mode=2 will only read
!//       a field's identification and not the field (data).
!//       to avoid (cray) backspace problems, the identification
!//       is stored until the next read.

!// warning: using file unit no. (not file name) to identify
!//          files when storing field identification.
!//          If more than one file is opened with the same
!//          unit, use 'call getflt(-1,...)' after closing
!//          a file to avoid errors.

!// no computer dependant i/o methodes
!//                 ------
!// Does also work with CRAY's from PrgEnv 3.0.2.1
!// and above, provided an assign -F f77 -N ieee u:iunit
!// is applied somewhere.
!// NTNU/ITEA 1998-06-23 Jorn Amundsen
!//                 ------

!// DNMI/FoU  12.04.1992  Anstein Foss
!// DNMI/FoU  19.08.1993  Anstein Foss
!// DNMI/FoU  16.04.1998  Anstein Foss ... cray.t3e
!// DNMI/FoU  02.07.1998  Anstein Foss ... integer*2 also for Cray

IMPLICIT NONE

INTEGER, PARAMETER :: mfsize=500000  ! Maximum field size, input and output

!// input/output

INTEGER            :: mode
INTEGER            :: iunit
INTEGER            :: ldata
INTEGER            :: igtype
INTEGER            :: nfskip         
INTEGER            :: ierr
INTEGER            :: itime(5)
INTEGER            :: ident(20)
REAL               :: fdata(ldata)
REAL               :: gparam(6)

!// local

INTEGER*2          :: ident2(20),idata2(mfsize)

INTEGER            :: itimef(5),idents(20)

INTEGER, PARAMETER :: maxsav=5, maxs20=maxsav*20
INTEGER            :: iusave(maxsav),idsave(20,maxsav)

INTEGER, PARAMETER :: mlgeom=50
INTEGER            :: lgtype
INTEGER*2          :: idgeom(20+1+mlgeom)
REAL               :: gplast(6)

INTEGER            :: i,nsave,isave,ios,iexit,iskip,n,ihours,ierr1,ierr2
INTEGER            :: nxin,nyin,nword,lgeom,nw,iscale,newgrid,ix,iy,ierror
REAL               :: scale

DATA iusave/maxsav*0/
DATA idsave/maxs20*0/

DATA lgtype/0/
DATA idgeom/71*0/
DATA gplast/6*999999./

IF (iunit < 1 .OR. iunit > 99 .OR.&
   (mode /= 0   .AND. mode /= 1   .AND. mode /= 2   .AND.&
    mode /= 100 .AND. mode /= 101 .AND. mode /= 102 .AND.&
    mode /= 200 .AND. mode /= 201 .AND. mode /= -1)) THEN
  PRINT *," **getflt** error ** mode,iunit: ",mode,iunit
  STOP
ENDIF

IF (mode == 200 .OR. mode.eq.201) THEN
  DO i=1,20
    idents(i)=ident(i)
  ENDDO
ENDIF

IF (mode == 201) REWIND(iunit)

!//      a) read field identification (one record)
!//      b) read field data           (one record)

ierr = 0
nfskip = 0

nsave = 0
isave = 0
DO i=1,maxsav
  IF (iusave(i) > 0)      nsave = nsave+1
  IF (iusave(i) == iunit) isave = i
ENDDO

IF (mode == -1) THEN
  IF (isave > 0) iusave(isave) = 0
  GOTO 990
ENDIF

IF (mode == 201 .AND. isave > 0) THEN
  iusave(isave) = 0
  isave = 0
  nsave = nsave-1
ENDIF

IF (isave > 0) THEN
  DO i=1,20
    ident(i) = idsave(i,isave)
  ENDDO
  iusave(isave) = 0
  nsave = nsave-1
  GOTO 30
ENDIF

20 CONTINUE

READ(iunit,IOSTAT=ios,ERR=910,END=920) (ident2(i),i=1,20)
DO i=1,20
  ident(i) = ident2(i)
ENDDO
!IF (ident(9) == 1010) ident(9) = 1

30 CONTINUE

!// ..iexit set to 1 if the next field is not read
!// ..iskip set to 1 if the next field is skipped
iexit = 0
iskip = 0

IF (mode == 100) THEN
  iexit = 1
ELSEIF (mode == 200 .OR. mode == 201) THEN
  n = 0
  DO i=1,20
    IF (idents(i) == -32767 .OR. idents(i) == ident(i)) n = n+1
  ENDDO
  IF (n /= 20) iskip = 1
ELSEIF (mode /= 0) THEN
  itimef(1) = ident(12)
  itimef(2) = ident(13)/100
  itimef(3) = ident(13)-(ident(13)/100)*100
  itimef(4) = ident(14)/100
  itimef(5) = ident(4)
  CALL hrdiff(0,0,itime,itimef,ihours,ierr1,ierr2)
  IF (mode == 1   .AND. ihours <=  0) iskip = 1
  IF (mode == 2   .AND. ihours /= 0)  iexit = 1
  IF (mode == 101 .AND. ihours <=  0) iskip = 1
  IF (mode == 101 .AND. ihours > 0)   iexit = 1
  IF (mode == 102 .AND. ihours < 0)   iskip = 1
  IF (mode == 102 .AND. ihours >=  0) iexit = 1
ENDIF

IF (iexit == 1) THEN
  !// Return without reading the field (data). save identification.
  IF (nsave == maxsav) iusave(1) = 0
  isave = 0
  DO n=1,maxsav
    IF (iusave(n) > 0) THEN
      isave = isave+1
      iusave(isave) = iusave(n)
      DO i=1,20
        idsave(i,isave) = idsave(i,n)
      ENDDO
    ENDIF
  ENDDO
  isave = maxsav
  iusave(isave) = iunit
  DO i=1,20
    idsave(i,isave) = ident(i)
  ENDDO
ENDIF

IF (iexit == 1) THEN
  ierr = 0
  IF (mode == 2) ierr=3
  GOTO 990
ENDIF

nxin = ident(10)
nyin = ident(11)
nword = nxin*nyin

lgeom = 0
IF (ident(9) >=  1000) lgeom=ident(9)-(ident(9)/1000)*1000

IF (nword > ldata .AND. iskip == 0) THEN
  PRINT *," **getflt** field length too big, (input ldata too small)"
  PRINT *," **ldata = ",ldata
  PRINT *," **ident: ",(ident(i),i=1,11)
  PRINT *," **       ",(ident(i),i=12,20)
  PRINT *," **nx,ny,nx*ny: ",nxin,nyin,nword
  ierr = 1
  GOTO 990
ENDIF

IF (nword+lgeom > mfsize) THEN
  PRINT *," **getflt** field length too big, (buffer too small)"
  PRINT *," **mfsize = ",mfsize
  PRINT *," **ident: ",(ident(i),i=1,11)
  PRINT *," **       ",(ident(i),i=12,20)
  PRINT *," **nx,ny,nx*ny,lgeom: ",nxin,nyin,nword,lgeom
  ierr = 1
  GOTO 990
ENDIF

nw = nword+lgeom
READ(iunit,IOSTAT=ios,ERR=910,END=920) (idata2(i),i=1,nw)	

IF (iskip == 0) THEN
  iscale = ident(20)
  scale = 10.**iscale
  !// Following (fast) loop does not handle 'undefined' values
  DO i=1,nword
    fdata(i) = scale*idata2(i)
  ENDDO
!//.............................................................
  !// The next (slow) loop handles 'undefined' values
!  DO i=1,nword
!    IF (idata2(i) /= -32767) THEN
!      fdata(i) = scale*idata2(i)
!    ELSE
!      fdata(i)=undef
!    ENDIF
!  ENDDO
ENDIF

IF (iskip /= 0) THEN
  nfskip = nfskip+1
  GOTO 20
ENDIF

IF (ident(9) > 0) THEN
  newgrid = 0
  IF (ident(9) /= idgeom(9)) THEN
    newgrid = 1
  ELSE
    DO i=15,18
      IF (ident(i) /= idgeom(i)) newgrid = 1
    ENDDO
    DO i=1,lgeom
      IF (idata2(nword+i) /= idgeom(21+i)) newgrid = 1
    ENDDO
  ENDIF
  IF (newgrid == 1) THEN
    DO i=1,20
      idgeom(i) = ident(i)
    ENDDO
    idgeom(10) = 1
    idgeom(11) = 1
    idgeom(21) = 0
    DO i=1,lgeom
      idgeom(21+i) = idata2(nword+i)
    ENDDO
    CALL gridpars(+1,20+1+mlgeom,idgeom,lgtype,ix,iy,gplast,ierror)
    IF (ierror /= 0) THEN
      PRINT *,"GETFLT: GRIDPARS ERROR ",ierror
      idgeom(9) = -9999
      ierr = 2
      GOTO 990
    ENDIF
  ENDIF
  igtype = lgtype
  DO i=1,6
    gparam(i) = gplast(i)
  ENDDO
ELSE
  igtype    = ident(9)
  gparam(1) = ident(15)
  gparam(2) = ident(16)
  gparam(3) = ident(17)
  gparam(4) = ident(18)
  gparam(5) = 0.
  gparam(6) = 0.
ENDIF

ierr = 0
GOTO 990

910 CONTINUE
  ierr = 1
  PRINT *," **getflt** read error. file,iostat: ",iunit,ios
  IF (mode == 200 .OR. mode == 201) GOTO 940
  GOTO 990

920 CONTINUE
  ierr = 2
  PRINT *," **getflt** end_of_file. file,iostat: ",iunit,ios
  IF (mode == 200 .OR. mode == 201) GOTO 940
  GOTO 990

940 CONTINUE
  DO i=1,20
    ident(i) = idents(i)
  ENDDO

990 CONTINUE

RETURN

END


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

SUBROUTINE putflt(iunit,ident,igtype,gparam,ldata,fdata,scale,ierr)

!// 16 bit pack and output (handling 'undefined')

!// input:
!//    iunit: file unit
!//    ident: field identification, (set ident(20)=-32767 for automatic scaling)
!//    igtype: grid type
!//    gparam: grid specification parameters
!//    ldata: length of the field
!//           (field size given in identification used when writing the field)
!//    fdata: the field, unscaled
!//    scale: additional scaling of the field before output, usually 1
!//           (this is in addition to the scaling given in the identification, 
!//            for changing the basic unit)

!// output:
!//    ierr = 0:  write o.k.
!//           1:  write error

!// no computer dependant i/o methodes.
!//                 ------
!// Does also work with CRAY's from PrgEnv 3.0.2.1
!// and above, provided an assign -F f77 -N ieee u:iunit
!// is applied somewhere.
!// NTNU/ITEA 1998-06-23 Jorn Amundsen
!//                 ------

!// DNMI/FoU  12.04.1992  Anstein Foss
!// DNMI/FoU  19.08.1993  Anstein Foss
!// DNMI/FoU  10.10.1997  Anstein Foss ... added extra geometry ident.
!// DNMI/FoU  16.04.1998  Anstein Foss ... cray.t3e
!// DNMI/FoU  02.07.1998  Anstein Foss ... integer*2 also for Cray

IMPLICIT NONE

INTEGER, PARAMETER :: mfsize=500000  ! Maximum field size, input and output
REAL,    PARAMETER :: undef=+1.e+35
REAL,    PARAMETER :: udef=undef*0.9

INTEGER, PARAMETER :: mldata=mfsize+50

!// input/output
INTEGER            :: iunit
INTEGER            :: igtype
INTEGER            :: ident(20)
INTEGER            :: ldata
INTEGER            :: ierr
REAL               :: gparam(6)
REAL               :: fdata(ldata)
REAL               :: scale

!// local
INTEGER            :: idento(20)
INTEGER*2          :: ident2(20),idata2(mldata)

INTEGER, PARAMETER :: mlgeom=50
INTEGER            :: lgtype,llgeom
INTEGER*2          :: idgeom(20+1+mlgeom)
REAL               :: gplast(6)

INTEGER            :: i,newgrid,ix,iy,ierror,lgeom,nxout,nyout,nword
INTEGER            :: sc,iscale,ifmax,ios
REAL               :: fmax

DATA lgtype,llgeom/-9999,0/
DATA idgeom/71*0/
DATA gplast/6*999999./

DO i=1,20
  idento(i) = ident(i)
ENDDO

IF (igtype > 0) THEN
  newgrid = 0
  IF (igtype /= lgtype) newgrid = 1
  DO i=1,6
    IF (gparam(i) /= gplast(i)) newgrid = 1
  ENDDO
  IF (newgrid == 1) THEN
    lgtype = igtype
    DO i=1,6
      gplast(i) = gparam(i)
    ENDDO
    ix = 1
    iy = 1
    CALL gridpars(-1,20+1+mlgeom,idgeom,igtype,ix,iy,gparam,ierror)
    IF (ierror /= 0) THEN
      PRINT *,"PUTFLT: GRIDPARS ERROR ",ierror
      lgtype = -999
      ierr = 2
      GOTO 990
    ENDIF
    llgeom = 0
    IF (idgeom(9) >=  1000) llgeom = idgeom(9)-(idgeom(9)/1000)*1000
  ENDIF
  idento(9)  = idgeom( 9)
  idento(15) = idgeom(15)
  idento(16) = idgeom(16)
  idento(17) = idgeom(17)
  idento(18) = idgeom(18)
  lgeom = llgeom
ELSE
  !// in case of "something" stored as if it were a field, "geometry"
  !// identification is left unchanged (extended "geometry" not possible)
  lgeom = 0
ENDIF

nxout = ident(10)
nyout = ident(11)
nword = nxout*nyout

ierr = 0

IF (nword+lgeom > mldata) THEN
  PRINT *," **putflt** field length too big, (buffer too small)"
  PRINT *," **mfsize = ",mfsize
  PRINT *," **mldata = ",mldata
  PRINT *," **ident: ",(idento(i),i=1,11)
  PRINT *," **       ",(idento(i),i=12,20)
  PRINT *," **nx,ny,nx*ny,lgeom: ",nxout,nyout,nword,lgeom
  ierr = 1
  GOTO 990
ENDIF

sc = scale
IF (sc == 0.) sc = 1.
iscale = ident(20)

IF (iscale == -32767) THEN
  !// ..automatic scaling
  fmax = 0.
  !// Fast code, not handling undefines values
!  DO i=1,nword
!    fmax = MAX(fmax,ABS(fdata(i)))
!  ENDDO
!//.............................................................
  !//  Slow code, handling undefines values
  DO i=1,nword
    IF (fdata(i) < udef) fmax = MAX(fmax,ABS(fdata(i)))
  ENDDO
  IF (fmax > 0.) THEN
    fmax = fmax*ABS(sc)
    iscale = LOG10(fmax)-4.
    ifmax = NINT(fmax*10.**(-iscale))
    IF (ifmax < 3278) THEN
      iscale = iscale-1
      ifmax = NINT(fmax*10.**(-iscale))
    ENDIF
    IF (ifmax > 32766) iscale = iscale+1
      iscale = MAX(iscale,-30)
  ELSE
    iscale = 0
  ENDIF
  idento(20) = iscale
ENDIF

!// Scale and 'pack' data
sc = sc*10.**(-iscale)

!// Fast code not handling undefines values
!DO i=1,nword
!  idata2(i)= NINT(sc*fdata(i))
!ENDDO
!//.............................................................
!// Slow code handling undefines values
DO i=1,nword
  IF (fdata(i) < udef) THEN
    idata2(i) = NINT(sc*fdata(i))
  ELSE
    idata2(i) = -32767
  ENDIF
ENDDO
!//.............................................................
!//Fixed scaling, check integer*2 value range
!IF (ident(20) /= -32767) THEN
!  DO i=1,nword
!    idata2(i) = MAX(MIN(idata2(i),+32767),-32768)
!  ENDDO
!ENDIF
!//.............................................................
IF (lgeom > 0) THEN
  DO i=1,lgeom
    idata2(nword+i) = idgeom(21+i)
  ENDDO
  nword = nword+lgeom
ENDIF

DO i=1,20
  ident2(i) = idento(i)
ENDDO

!// Write field identification and field data

WRITE(iunit,IOSTAT=ios,ERR=900) (ident2(i),i=1,20)
WRITE(iunit,IOSTAT=ios,ERR=900) (idata2(i),i=1,nword)

GOTO 990

900 CONTINUE
  ierr = 1
  PRINT *," **putflt** write error. file,iostat: ",iunit,ios

990 CONTINUE

RETURN

END


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

SUBROUTINE hrdiff(iup1,iup2,itime1,itime2,ihours,ierr1,ierr2)

!// Calculate interval in hours between 'itime1' and 'itime2',
!// ihours = 'itime2' - 'itime1'  (positive,zero or negative)

!// itime: year,month,day,time in hours, forecast length in hours
!// (positive,zero or negative)

!// If iup1=1: 'itime1' is updated to give valid date,time
!// If iup2=1: 'itime2' is updated to give valid date,time
!//    (forecast length = 0)

!// ierr1,ierr2: 0 = o.k. itime1,itime2
!//              1 = not o.k. itime1,itime2

!// DNMI/FoU  xx.xx.1990  Anstein Foss

INTEGER :: itime1(5)
INTEGER :: itime2(5)

INTEGER :: mdays(12)
INTEGER :: it(5,2)

DATA mdays/31,28,31,30,31,30,31,31,30,31,30,31/

DO i=1,5
  it(i,1) = itime1(i)
  it(i,2) = itime2(i)
ENDDO

!// Put prog.time into year,month,day,time

CALL vtime(it(1,1),ierr1)
CALL vtime(it(1,2),ierr2)

!//  'it' now gives "veryfing" time of 'itime1' and 'itime2'

IF (iup1 == 1 .AND. ierr1 == 0) THEN
  DO i=1,5
    itime1(i) = it(i,1)
  ENDDO
ENDIF
IF (iup2 == 1 .AND. ierr2 == 0) THEN
  DO i=1,5
    itime2(i) = it(i,2)
  ENDDO
ENDIF

IF (ierr1 /= 0 .OR. ierr2 /= 0) THEN
  ihours = -32767
  RETURN
ENDIF

DO i=1,4
  IF (it(i,1) /= it(i,2)) GOTO 55
ENDDO

!// No time difference
ihours = 0
RETURN

55 CONTINUE

nt1 = 1
nt2 = 2
IF (it(i,1) > it(i,2)) THEN
  nt1 = 2
  nt2 = 1
ENDIF
nhh = 0

IF (it(1,nt1) == it(1,nt2)) GOTO 70
iy = it(1,nt1)
!// Remaining hours first year:
DO im=it(2,nt1),12
  md = mdays(im)
  IF (im == 2) THEN
    IF (iy/4*4 == iy) md = 29
    IF (iy/100*100 == iy .AND. iy/400*400 /= iy) md = 28
  ENDIF
  nhh = nhh+md*24
ENDDO
nhh = nhh-(it(3,nt1)-1)*24-it(4,nt1)
!// One year steps
DO iy=it(1,nt1)+1,it(1,nt2)-1
  nd = 365
  IF (iy/4*4 == iy) nd = 366
  IF (iy/100*100 == iy .AND. iy/400*400 /= iy) nd = 365
  nhh = nhh+nd*24
ENDDO
it(1,nt1) = it(1,nt2)
it(2,nt1) = 1
it(3,nt1) = 1
it(4,nt1) = 0

70 CONTINUE

IF (it(2,nt1) == it(2,nt2)) GOTO 80
!// Remaining hours first month
iy = it(1,nt1)
im = it(2,nt1)
md = mdays(im)
IF (im == 2) THEN
  IF (iy/4*4 == iy) md = 29
  IF (iy/100*100 == iy .AND. iy/400*400 /= iy) md = 28
ENDIF
nhh = nhh+(md-it(3,nt1)+1)*24-it(4,nt1)
!// One month steps
DO im=it(2,nt1)+1,it(2,nt2)-1
  md = mdays(im)
  IF (im == 2) THEN
    IF (iy/4*4 == iy) md = 29
    IF (iy/100*100 == iy .AND. iy/400*400 /= iy) md = 28
  ENDIF
  nhh = nhh+md*24
ENDDO
it(2,nt1) = it(2,nt2)
it(3,nt1) = 1
it(4,nt1) = 0

80 CONTINUE

IF (it(3,nt1) /= it(3,nt2)) THEN
  nhh = nhh+(it(3,nt2)-it(3,nt1))*24-it(4,nt1)
  it(3,nt1) = it(3,nt2)
  it(4,nt1) = 0
ENDIF
!// Hours last day
nhh = nhh+it(4,nt2)-it(4,nt1)
it(4,nt1) = it(4,nt2)

IF (nt1 == 1) THEN
  ihours = nhh
ELSE
  ihours = -nhh
ENDIF

RETURN

END


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

SUBROUTINE vtime(itime,ierror)

!// 'itime' is updated to give "valid" date,time (with prog.time = 0)

!// input:  itime(5) - itime(1): year
!//  itime(2): month (1-12)
!//  itime(3): day (1-28/29/30/31)
!//  itime(4): time in hours (00-23)
!//  itime(5): time in hours of prognosis
!//            (negative, zero or positive)
!// output: itime(5) -  as above, itime(5) = 0
!//   ierror   -  0 = o.k. input date/time
!//               1 = not o.k. input date/time
!// ('itime' not changed)

!// DNMI/FoU  xx.xx.1992  Anstein Foss

INTEGER :: itime(5)

INTEGER :: mdays(12)

DATA mdays/31,28,31,30,31,30,31,31,30,31,30,31/

iy = itime(1)
im = itime(2)
id = itime(3)
ih = itime(4)

!// Test input time
ierror = 0
IF (im < 1 .OR. im > 12) THEN
  ierror = 1
ELSE
  md = mdays(im)
  IF (im == 2) THEN
    IF (iy/4*4 == iy) md = 29
    IF (iy/100*100 == iy .AND. iy/400*400 /= iy) md = 28
  ENDIF
  IF (id < 1 .OR. id > md) ierror = 1
ENDIF
IF (ih < 00 .OR. ih > 23) ierror = 1

IF (ierror /= 0) RETURN

ih = ih+itime(5)
IF (ih >=  0 .AND. ih <=  23) GOTO 50
IF (ih < 0) GOTO 30

nd = ih/24
ih = ih-24*nd
DO n=1,nd
  id = id+1
  IF (id > md) THEN
    im = im+1
    IF (im > 12) THEN
      iy = iy+1
      im = 1
      md = mdays(im)
    ELSEIF (im == 2) THEN
      md = mdays(im)
      IF (iy/4*4 == iy) md = 29
      IF (iy/100*100 == iy .AND. iy/400*400 /= iy) md = 28
    ELSE
      md = mdays(im)
    ENDIF
    id = 1
  ENDIF
ENDDO

GOTO 50

30 CONTINUE

nd = (-ih+23)/24
ih = ih+24*nd
DO n=1,nd
  id = id-1
  IF (id < 1) THEN
    im = im-1
    IF (im < 1) THEN
      iy = iy-1
      im = 12
      md = mdays(im)
    ELSEIF(im == 2) THEN
      md = mdays(im)
      IF (iy/4*4 == iy) md = 29
      IF (iy/100*100 == iy .AND. iy/400*400 /= iy) md = 28
    ELSE
      md = mdays(im)
    ENDIF
    id = md
  ENDIF
ENDDO

50 CONTINUE

itime(1) = iy
itime(2) = im
itime(3) = id
itime(4) = ih
itime(5) = 0

RETURN

END


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

SUBROUTINE gridpars(icall,ldata,idata,igtype,nx,ny,grid,ierror)

!// NAME:
!//    gridpar

!// PURPOSE:
!//    Conversion between integer*2 field identification and
!//    variables used in programs. Note that some gridtypes also
!//    has extended (geometry) identification behind the field data.

!// SYNOPSIS:
!//    SUBROUTINE gridpar(icall,ldata,idata,igtype,nx,ny,grid,ierror)
!//    integer   icall,ldata,igtype,nx,ny,ierror
!//    integer*2 idata(ldata)
!//    realgrid(6)

!// INPUT:
!//    icall  - +1 : from idata to igtype,nx,ny,grid
!//             -1 : from igtype,nx,ny,grid to idata
!//    ldata  - length of the idata array

!// INPUT/OUTPUT:
!//    idata  - field identification, field data (not touched)
!//             and possibly extra geometry specifications
!//             (max 20 words in this version)
!//    igtype - grid type, 1 = polarstereographic grid (true at 60N)
!//                        2 = geographic
!//                        3 = spherical rotated grid
!//                        4 = polarstereographic grid
!//                        5 = mercator grid (unrotated)
!//                        * = unknown grid type (-32767 to +32 accepted)
!//             (note that 'codes' are used in the field identifcation
!//             when extra geometry identification is stored after
!//             the field data, the reason for the +32 limit)
!//    nx     - no. of gridpoints in x direction (1 - 32767)
!//    ny     - no. of gridpoints in y direction (1 - 32767)
!//    grid   - grid parameters

!//    igtype=1,4, polarstereographic grid:
!//      grid(1) = x position of north pole (xp)
!//      grid(2) = y position of north pole (yp)
!//      grid(3) = no. of grid units between North Pole and Equator
!//      grid(4) = grid rotation angle (degrees), positive east, negative west
!//      grid(5) = projection latitude (degrees), standard is 60 (60 deg. N)
!//      grid(6) = 0. (not used)

!//    igtype=2,3, geographic or spherical rotated grid:
!//      grid(1) = western boundary (degrees) (longitude for x=1)
!//      grid(2) = southern boundary (degrees) (latitude for y=1)
!//      grid(3) = longitude increment (degrees)
!//      grid(4) = latitude  increment (degrees)
!//      grid(5) = longitude position of rotated equator (degrees) 
!//                (0 if geographic grid)
!//      grid(6) = latitude position of rotated equator (degrees)
!//                (0 if geographic grid)

!//    igtype=5, mercator (unrotated) grid:
!//      grid(1) = western boundary (degrees) (longitude for x=1)
!//      grid(2) = southern boundary (degrees) (latitude for y=1)
!//      grid(3) = x (longitude) increment (km)
!//      grid(4) = y (latitude)  increment (km)
!//      grid(5) = reference (construction) latitude (degrees)
!//      grid(6) = 0. (not used)

!//    igtype=*, unknown grid type, only use grid type less than 1 if the
!//              grid parameters have no meaning:
!//      grid(1:6) : unknown grid parameters

!// OUTPUT:
!//    ierror - error status: 0 = no error
!//                           1 = bad value in input identification
!//                               or in input grid parameters etc.
!//                           2 = ldata too small
!//                           3 = unknown icall
!//                           4 = ldata too small for needed extra
!//                               geometry identification (icall=-1),
!//                               but the best possible identification
!//                               is done (see NOTES below)

!// DOCUMENTATION:
!//    Basic document on felt files:
!//   FILE STRUKTUR FOR "SANNTIDS" LAGRING AV GRID-DATA
!//  Forskningsavdeling DNMI, oktober 1982
!//    See also /usr/local/doc/felt.doc (at DNMI)

!// NOTES:
!//  - This routine maintain compability with old formats
!//    (see comments in the source code),
!//    new formats added when required for some reason.
!//  - Specyfing ldata too short to hold extra geometry identification
!//    (for icall=-1) will restrict this routine from making this even
!//    when it seems appropriate. In order to be compatible with
!//    old unextended (or less extended) formats you then may
!//    disregard the returned ierror=4 (in special cases).
!//  - Avoid calling this routine more often than necessary,
!//    i.e. only after reading the first field or only before output
!//    of the first field.
!//    In models with several calling routines you may do a 'setup' call
!//    and keep (part of) the first 20 words of identification and
!//    possibly extra geometry specification locally (see note below).
!//  - The value of nx and ny (in idata if icall=+1) does not influence
!//    conversion of grid parameters.
!//    You may then use this routine to convert between field
!//    identification and grid parameters with nx=ny=1 and possibly
!//    find extra identification in word 22,23,... in idata.
!//    (see DOCUMENTATION to find format description).

!// DNMI/FoU  05.05.1995  Anstein Foss
!// DNMI/FoU  09.06.1995  Anstein Foss
!// DNMI/FoU  14.05.1996  Anstein Foss ... mercator (unrotated)
!// DNMI/FoU  02.09.1996  Anstein Foss ... gridtype 2012 -> 2 (bad test)
!// DNMI/FoU  15.10.1996  Anstein Foss ... even better scaling when needed
!// DNMI/FoU  17.02.1997  Anstein Foss ... and again (for 'image' fields)

IMPLICIT NONE

!// Input/output
INTEGER   :: icall
INTEGER   :: ldata
INTEGER   :: igtype
INTEGER   :: nx,ny
INTEGER   :: ierror
INTEGER*2 :: idata(ldata)
REAL      :: grid(6)

!// Local
INTEGER   :: ngw,lgeom,i,igr,ig1,ig2,ld,igscale,kgeom,lgeom1,lgeom2
INTEGER*2 :: igeom1(12),igeom2(20)
REAL      :: gscale,glog,gws,dglim,dgmax,dgprv,grx
REAL      :: gw(6),gr(6)

ierror = 0

IF (icall == +1) THEN   ! idata(ldata) -> igtype,nx,ny,grid(6)

  igtype = 0
  nx = 0
  ny = 0
  DO i=1,6
    grid(i) = 0.
  ENDDO
  IF (ldata < 20) THEN
    ierror = 2
    RETURN
  ENDIF
  igtype = idata(9)
  lgeom = 0
  IF (igtype > 999) THEN
    i = igtype
    igtype = igtype/1000
    lgeom = i-igtype*1000
  ENDIF
  nx = idata(10)
  ny = idata(11)
  IF (nx < 1 .OR. ny < 1) THEN
    ierror = 1
    RETURN
  ENDIF

  IF (igtype == 1 .OR. igtype == 4) THEN

    !// Polarstereographic (type 1 always true at 60 degrees North)

    IF (idata(17) > 0) THEN
      !// The standard
      grid(1) = idata(15)*0.01
      grid(2) = idata(16)*0.01
      grid(3) = idata(17)*0.1
      grid(4) = idata(18)
    ELSE
      !// An old extension
      grid(1) = idata(15)
      grid(2) = idata(16)
      grid(3) = -idata(17)*0.1
      grid(4) = idata(18)
    ENDIF
    grid(5) = 60.

  ELSEIF (igtype == 2 .OR. igtype == 3) THEN

    !// Geographic (2) or spherical rotated grid (3)

    grid(1) = idata(16)*0.01
    grid(2) = idata(15)*0.01
    grid(3) = idata(18)*0.01
    grid(4) = idata(17)*0.01

  ELSEIF (igtype == 5) THEN

    !// Mercator grid (5)

    grid(1) = idata(15)*0.01
    grid(2) = idata(16)*0.01
    grid(3) = idata(17)*0.1
    grid(4) = idata(18)*0.1

  ELSE

    !// Unknown/undefined grid type

    grid(1) = idata(15)
    grid(2) = idata(16)
    grid(3) = idata(17)
    grid(4) = idata(18)

  ENDIF

  IF (lgeom > 0) THEN
    IF (igtype == 1 .OR. igtype == 4) THEN
      gscale = 100.
      ngw = 5
    ELSEIF (igtype == 2 .OR. igtype == 3) THEN
      gscale = 10000.
      ngw = 6
    ELSEIF (igtype == 5) THEN
      gscale = 10000.
      ngw = 5
    ELSE
      gscale = 100.
      ngw = 6
    ENDIF

    ld=20+nx*ny

    IF (lgeom == ngw*2 .AND. ld+lgeom <=  ldata) THEN

      !// First extended method
      DO i=1,ngw
        ig1 = idata(ld+1)
        ig2 = idata(ld+2)
        ld = ld+2
        grid(i) = FLOAT(ig1*10000+ig2)/gscale
      ENDDO

    ELSEIF (lgeom == 2+ngw*3 .AND. ld+lgeom <=  ldata) THEN

      !// Second extended method
      IF (idata(ld+1) == ngw .AND. idata(ld+2) == 3) THEN
        ld=ld+2
        DO i=1,ngw
          igscale = idata(ld+1)
          ig1 = idata(ld+2)
          ig2 = idata(ld+3)
          ld = ld+3
          gscale = 10.**igscale
          grid(i) = FLOAT(ig1*10000+ig2)*gscale
        ENDDO
      ELSE
        ierror = 2
      ENDIF

    ELSE

      ierror = 2

    ENDIF

  ENDIF

  IF (ierror == 0) THEN
    IF (igtype == 1 .OR. igtype == 4) THEN
      !// The DNMI "standard" (150km grid => an=grid(3)=79.)
      IF (grid(3) /= 0.) grid(3) = 79.*150./grid(3)
      IF (grid(3) == 0.  .OR. grid(5) == 0. .OR.&
          grid(5) < -90. .OR. grid(5) > +90.) ierror = 1
    ELSEIF (igtype == 2 .OR. igtype == 3) THEN
      IF (grid(3) == 0. .OR. grid(4) == 0.) ierror = 1
    ELSEIF (igtype == 5) THEN
      IF (grid(3) == 0. .OR. grid(4) == 0.) ierror = 1
    ENDIF
  ENDIF

ELSEIF (icall == -1) THEN   ! igtype,nx,ny,grid(6) -> idata(ldata)

  IF (ldata < 20) THEN
    ierror = 2
    RETURN
  ENDIF

  IF (igtype > 32 .OR. nx < 1 .OR. nx > 32767 .OR. ny < 1 .OR. ny > 32767) THEN
    ierror = 1
    RETURN
  ENDIF

  idata(9) = igtype
  idata(10) = nx
  idata(11) = ny
  DO i=1,6
    gw(i) = grid(i)
  ENDDO

  IF (igtype == 1 .OR. igtype == 4) THEN

    !// Polarstereographic (type 1 always true at 60 degrees North)

    IF (gw(3) == 0.) THEN
      ierror = 1
      RETURN
    ENDIF

    !// The DNMI "standard" (150km grid => an=grid(3)=79.)
    gw(3) = 150.*79./gw(3)

    IF (ABS(gw(1)) < 327.66 .AND. ABS(gw(2)) < 327.66) THEN
      !// The standard
      idata(15) = NINT(gw(1)*100.)
      idata(16) = NINT(gw(2)*100.)
      idata(17) = NINT(gw(3)*10.)
      idata(18) = NINT(gw(4))
      gr(1) = FLOAT(idata(15))*0.01
      gr(2) = FLOAT(idata(16))*0.01
      gr(3) = FLOAT(idata(17))*0.1
      gr(4) = FLOAT(idata(18))
      gr(5) = 60.
    ELSEIF(ABS(gw(1)) < 32766. .AND. ABS(gw(2)) < 32766.) THEN
      !// An old extension
      idata(15) = NINT(gw(1))
      idata(16) = NINT(gw(2))
      idata(17) = -NINT(gw(3)*10.)
      idata(18) = NINT(gw(4))
      gr(1) = FLOAT(idata(15))
      gr(2) = FLOAT(idata(16))
      gr(3) = -FLOAT(idata(17))*0.1
      gr(4) = FLOAT(idata(18))
      gr(5) = 60.
    ELSE
      !// Old impossible case
      idata(15) = 32767
      idata(16) = 32767
      idata(17) = 32767
      idata(18) = 32767
      gr(1) = 327.67
      gr(2) = 327.67
      gr(3) = 3276.7
      gr(4) = 32767.
      gr(5) = 60.
    ENDIF

    gscale = 100.
    ngw = 5

  ELSEIF (igtype == 2 .OR. igtype == 3) THEN

    !// Geographic (2) or spherical rotated grid (3)

    IF (gw(3) == 0. .OR. gw(4) == 0.) THEN
      ierror = 1
      RETURN
    ENDIF

    IF (gw(1) > +180.) gw(1) = gw(1)-360.
    IF (gw(1) < -180.) gw(1) = gw(1)+360.

    idata(15) = NINT(gw(2)*100.)
    idata(16) = NINT(gw(1)*100.)
    idata(17) = NINT(gw(4)*100.)
    idata(18) = NINT(gw(3)*100.)
    gr(1) = FLOAT(idata(16))*0.01
    gr(2) = FLOAT(idata(15))*0.01
    gr(3) = FLOAT(idata(18))*0.01
    gr(4) = FLOAT(idata(17))*0.01
    gr(5) = 0.
    gr(6) = 0.

    gscale = 10000.
    ngw = 6

  ELSEIF (igtype == 5) THEN
 
    !// Mercator grid (5)

    IF (gw(3) == 0. .OR. gw(4) == 0.) THEN
      ierror = 1
      RETURN
    ENDIF

    idata(15) = NINT(gw(1)*100.)
    idata(16) = NINT(gw(2)*100.)
    idata(17) = NINT(gw(3)*10.)
    idata(18) = NINT(gw(4)*10.)
    gr(1) = FLOAT(idata(15))*0.01
    gr(2) = FLOAT(idata(16))*0.01
    gr(3) = FLOAT(idata(17))*0.1
    gr(4) = FLOAT(idata(18))*0.1
    gr(5) = 0.

    gscale = 10000.
    ngw = 5

  ELSE

    !// Unknown/undefined grid type

    idata(15) = NINT(gw(1))
    idata(16) = NINT(gw(2))
    idata(17) = NINT(gw(3))
    idata(18) = NINT(gw(4))
    gr(1) = FLOAT(idata(15))
    gr(2) = FLOAT(idata(16))
    gr(3) = FLOAT(idata(17))
    gr(4) = FLOAT(idata(18))
    gr(5) = 0.
    gr(6) = 0.

    gscale = 100.
    ngw = 6

  ENDIF

  !// Check if the standard packing above was good enough
  !// or if the first or second extended method should be used,
  !// for compability with old formats the least extended is preferred

  !// A limit to avoid high precision complications
  dglim = 1.e-8

  dgmax = 0.

  IF (igtype > 0) THEN
    DO i=1,ngw
      IF (gw(i) /= 0.) dgmax = MAX(dgmax,ABS((gr(i)-gw(i))/gw(i)))
    ENDDO
  ENDIF

  IF (dgmax > dglim) THEN

    kgeom = 0
    ld = 20+nx*ny

  !// First extended method, same scaling for all grid parameters
  !// but don't use it unless it gives better results
  !// (kept in the code for compability with old formats, possibly
  !//  in models etc. not using this routine)

    dgprv = dgmax
    dgmax = 0.
    lgeom1 = 0

    DO i=1,ngw
      gws = gw(i)*gscale
      !// Check overflow
      IF (gws /= 0. .AND. ABS(gws) < 3.e+8) THEN
        igr = NINT(gws)
        ig1 = igr/10000
        ig2 = igr-ig1*10000
        igeom1(lgeom1+1) = ig1
        igeom1(lgeom1+2) = ig2
        grx = FLOAT(ig1*10000+ig2)/gscale
        dgmax = MAX(dgmax,ABS((grx-gw(i))/gw(i)))
      ELSE
        !// Zero or value not possible for this method
        igeom1(lgeom1+1) = 0
        igeom1(lgeom1+2) = 0
        IF (gws /= 0.) dgmax = 1.
      ENDIF
      lgeom1 = lgeom1+2
    ENDDO

    IF (dgmax < dgprv) THEN
      IF (ldata >= ld+lgeom1) THEN
        kgeom = 1
      ELSE
        ierror = 4
      ENDIF
    ELSE
      dgmax = dgprv
    ENDIF

    IF (dgmax > dglim .AND. ierror == 0) THEN

      !// Second extended method, good enough for any real*4 precision value
      !// but don't use it unless it gives better results

      dgprv = dgmax
      dgmax = 0.

      igeom2(1) = ngw
      igeom2(2) = 3
      lgeom2 = 2
      DO i=1,ngw
        IF (gw(i) /= 0.) THEN
          !// 8 decimals precision (more than enough for real*4)
          glog = LOG10(ABS(gw(i)))-8.
          igscale = INT(glog)
          IF (glog > 0.) igscale = igscale+1
          !// Keep scaling within real*4 range
          igscale = MAX(MIN(igscale,+25),-25)
          gscale = 10.**(-igscale)
          igr = NINT(gw(i)*gscale)
          ig1 = igr/10000
          ig2 = igr-ig1*10000
          igeom2(lgeom2+1) = igscale
          igeom2(lgeom2+2) = ig1
          igeom2(lgeom2+3) = ig2
          grx = FLOAT(ig1*10000+ig2)/gscale
          dgmax = MAX(dgmax,ABS((grx-gw(i))/gw(i)))
        ELSE
          igeom2(lgeom2+1) = 0
          igeom2(lgeom2+2) = 0
          igeom2(lgeom2+3) = 0
        ENDIF
        lgeom2 = lgeom2+3
      ENDDO

      IF (dgmax < dgprv) THEN
        IF (ldata >= ld+lgeom2) THEN
          kgeom = 2
        ELSE
          ierror = 4
        ENDIF
!      ELSE
!        dgmax=dgprv
      ENDIF

    ENDIF

    IF (kgeom == 1) THEN
      idata(9) = igtype*1000+lgeom1
      DO i=1,lgeom1
        idata(ld+i) = igeom1(i)
      ENDDO
    ELSEIF (kgeom == 2) THEN
      idata(9) = igtype*1000+lgeom2
      DO i=1,lgeom2
        idata(ld+i) = igeom2(i)
      ENDDO
    ENDIF

  ENDIF


ELSE

  !// Wrong icall

  ierror = 3

ENDIF

RETURN

END