program timeseries
implicit none

!     This program extracts a number of meteorological quantities from HIRLAM data
!     on nxn grid cell squares for selected regions of the Netherlands,
!     and CONStructs TimeSeries from them.

!     Compilation:
!
!     gfortran -g -Ofast -o timeseries timeseries.f90 port.a gcod.a vari.a

integer, parameter :: jpfclen = 54, jpfcint = 1
integer, parameter :: jpb1 = 30, jpb2 = 40, jpsup = 2, jplen = 1600000
integer, parameter :: maxregions = 100, maxregside = 5
real, parameter    :: zrad = 180.0 / (4.0 * atan(1.0))

! Region declarations.

integer      :: nregs
character*30 :: regname(maxregions)
integer      :: indxreg(maxregside,maxregside,maxregions)
integer      :: iregctr(maxregions), nregside(maxregions)
real         :: reglat(maxregions), reglon(maxregions), regdir(maxregions)

namelist /regionlist/ regname, reglat, reglon, nregside

! Observation declarations.

integer    :: jobs, nobs
type observation
   integer :: ilat, ilon
   real    :: height, wdir, wspd, temp, dewp, rh, vis, cc, rr, gust, snowdepth
end type observation
type(observation) :: obs(100000)

! ASIMOF I/O declarations.

integer :: ibpw = 0, ibpwio = 0, imdi = 0  ! ASIMOF parameters
integer :: ilen, ierr, indxwnd
integer :: ib1(jpb1), ib2(jpb2), zb2(jpb2), zsup(jpsup)
real    :: field(jplen)

! Grid declarations.

integer :: nlon = 0, nlat = 0
real    :: east, west, south, north, polat, polon, dlat, dlon

! Meteorological quantities for time series.

type modelfield
   integer   :: block1(3)
   character :: format*8
   character :: name*12
end type modelfield
type(modelfield), parameter :: fld(*) = [                  &
   modelfield([ 11, 105,   2], 'f8.2',    'TEMPERATUUR' ), &
   modelfield([ 17, 105,   2], 'f8.2',    'DAUWPUNT'    ), &
   modelfield([ 33, 105,  10], 'f8.2',    'WINDRICHTING'), &
   modelfield([ 34, 105,  10], 'f8.2',    'WINDSNELHEID'), &
   modelfield([ 62, 105,   0], 'f8.2',    'FRNTNEERSLAG'), &
   modelfield([ 63, 105,   0], 'f8.2',    'CONVNEERSLAG'), &
   modelfield([ 71, 105,   0], 'f8.2',    'BEWOLKING'   ), &
   modelfield([ 20, 105,   0], 'f8.1',    'ZICHT'       ), &
   modelfield([ 52, 105,   2], 'f8.3',    'RELVOCHT'    ), &
   modelfield([224, 105,  32], 'f8.1',    'CIN'         ), &
   modelfield([225, 105,  32], 'f8.1',    'CAPE'        ), &
   modelfield([226, 105,   0], 'f8.2',    'NEERSLAGSRT' ), &
   modelfield([228, 105,  10], 'f8.2',    'WINDSTOTEN'  ), &
   modelfield([ 66, 105,   0], 'f8.3',    'SNEEUWMASSA' ), &
   modelfield([138, 105,   0], 'f8.3',    'SNEEUWHOOGTE'), &
   modelfield([121, 105, 904], 'f8.1',    'LATHFLUX'    ), &
   modelfield([122, 105, 904], 'f8.1',    'SENSHFLUX'   ), &
   modelfield([111, 105,   0], 'es9.2e1', 'ZONLICHT'    ), &
   modelfield([112, 105,   0], 'es9.2e1', 'INFRAROOD'   )  &
]
integer, parameter :: jpval = size(fld)
real               :: values(maxregside**2,maxregions,jpval,0:jpfclen/jpfcint)
character*30       :: format

! Time series declarations.

integer :: npts, ival, itloff                ! Number of points, value index, time offset.
integer :: jtl, jreg, jval, jpts, jsq1, jsq2 ! Loop counts.
real    :: val, prv

! "Weercijfer" declarations (td=today, tm=tomorrow, da=day after; l=low, h=high).

integer     :: iwctdl, iwctdh, iwctml, iwctmh, iwcdal, iwcdah
integer     :: iwctd(maxregside**2), iwctm(maxregside**2), iwcda(maxregside**2)
character*7 :: chwctd, chwctm, chwcda

! Input file declarations.

integer       :: iunit, iostatus
character*256 :: file
character*10  :: dtg

!     Get sole argument of this program

if (command_argument_count() /= 1) then
   write(0,*)'TIMESERIES: SPECIFY DTG ON COMMAND LINE'
   call abort
endif
call get_command_argument(1, dtg)

!     Read namelist with region attributes

regname = ''
regdir  = 0.0
read(*, nml = regionlist)
nregs = count(regname > '')

!     Set ASIMOF parameters.

call ASIMHM(ibpw, ibpwio, imdi)

!     Loop over forecast times.

do jtl = 0, jpfclen, jpfcint

!        Read the observation file for this hour, if any.

   nobs = 0
   iunit = 10 + jtl
   open(unit = iunit, status = 'old', iostat = iostatus)
   if (iostatus == 0) then
      do
         nobs = nobs + 1
         read(unit = iunit, fmt = *, iostat = iostatus) &
            obs(nobs) % ilon, obs(nobs) % ilat,         &
            obs(nobs) % height,                         &
            obs(nobs) % wdir, obs(nobs) % wspd,         &
            obs(nobs) % temp, obs(nobs) % dewp,         &
            obs(nobs) % rh, obs(nobs) % vis,            &
            obs(nobs) % cc, obs(nobs) % rr,             &
            obs(nobs) % gust, obs(nobs) % snowdepth
         if (iostatus /= 0) then
            nobs = nobs - 1
            exit
         endif
      enddo
      close(unit = iunit)
   endif
   print*, 'uur:',jtl,' aantal waarnemingen:',nobs

!        Determine the forecast file name for this hour (it'd better be there).

   write(file, '(13a,i2.2,a)')        &
     '/data/hirlam/hl_arc/OPR/',      &
     dtg(1:4),                        &
     '/',                             &
     dtg(5:6),                        &
     '/',                             &
     dtg(7:8),                        &
     '/',                             &
     dtg(9:10),                       &
     '/fc',                           &
     dtg(1:8),                        &
     '_',                             &
     dtg(9:10),                       &
     '+0',                            &
     jtl,                             &
     've'

!        Open that ASIMOF file.

   call ASIMHO(1, file, 0, ierr)
   if (ierr /= 0) then
      write(0,*)'TIMESERIES - ERROR ON OPENING FILE: ', trim(file), '; RETURN CODE =', ierr
      call abort
   endif

!        Read orography from file to determine grid.

   if (nlon == 0) then

      ib1 = imdi
      ib1(9:11) = [6, 105, 0]

      ilen = jplen
      ierr = 0
      call GETFD(1, ib1, jpb1, ib2, jpb2, zb2, jpb2, zsup, jpsup, field, ilen, 0, ierr)
      if (ierr /= 0) then
         write(0,*)'TIMESERIES - ERROR ON READING FIELD:', ib1(9:11), ' RETURN CODE =', ierr
         call abort
      endif

      nlon  = ib2(7)
      nlat  = ib2(9)
      south = ib2(11) / 1000.0
      west  = ib2(14) / 1000.0
      north = ib2(18) / 1000.0
      east  = ib2(21) / 1000.0
      dlon  = (east  - west ) / real(nlon)
      dlat  = (north - south) / real(nlat)
      polat = ib2(33) / 1000.0
      polon = ib2(36) / 1000.0

!           Determine index in HIRLAM field for the center of each of the regions.

      call getnnb(reglon, reglat, regdir, west, south, dlon, dlat, polon, polat, nlon, nlat, iregctr, nregs)

!           Subsequently, fill the square of each region with indices of the other points.

      do jreg = 1, nregs

         if (mod(nregside(jreg), 2) /= 1) then
            write(0,*)'TIMESERIES: SQUARE SHOULD HAVE ODD SIDES'
            call abort
         endif
         if (nregside(jreg) > maxregside) then
            write(0,*)'TIMESERIES: SQUARE SIDE EXCEEDS MAXIMUM:', maxregside
            call abort
         endif
         do jsq1 = 1, nregside(jreg)
            do jsq2 = 1, nregside(jreg)
               indxreg(jsq1,jsq2,jreg) =                      &
                  iregctr(jreg) +                             &
                     nlon * (jsq1 - (nregside(jreg)+1)/2) +   &
                            (jsq2 - (nregside(jreg)+1)/2)
            enddo
         enddo

      enddo

   endif

!        Now that we have the indices of the region, gather observations for them.

   do jreg = 1, nregs
      do jobs = 1, nobs
         if (         & ! Check if coordinates of the observation(s) ...
             any(     & ! ... are in the square determining this region.
                 (obs(jobs)%ilat-1) * nlon + obs(jobs)%ilon ==      &
                 indxreg(1:nregside(jreg), 1:nregside(jreg), jreg)  &
                )                                                   &
            ) then      ! If so, write their values to the related observation files.
            if (obs(jobs) % wdir < 1.0e+38) then
               open(7, file=trim(regname(jreg))//'.WINDRICHTING-obs', position='append')
               write(7, fmt = *) jtl, obs(jobs) % wdir, 0
               close(7)
               print*,jtl,regname(jreg),'WINDRICHTING', obs(jobs) % ilon, obs(jobs) % ilat, obs(jobs) % wdir
            endif
            if (obs(jobs) % wspd < 1.0e+38) then
               open(7, file=trim(regname(jreg))//'.WINDSNELHEID-obs', position='append')
               write(7, fmt = *) jtl, obs(jobs) % wspd, 0
               close(7)
               print*,jtl,regname(jreg),'WINDSNELHEID', obs(jobs) % ilon, obs(jobs) % ilat, obs(jobs) % wspd
            endif
            if (obs(jobs) % temp < 1.0e+38) then
               open(7, file=trim(regname(jreg))//'.TEMPERATUUR-obs', position='append')
               write(7, fmt = *) jtl, obs(jobs) % temp, 0
               close(7)
               print*,jtl,regname(jreg),'TEMPERATUUR', obs(jobs) % ilon, obs(jobs) % ilat, obs(jobs) % temp
            endif
            if (obs(jobs) % dewp < 1.0e+38) then
               open(7, file=trim(regname(jreg))//'.DAUWPUNT-obs', position='append')
               write(7, fmt = *) jtl, obs(jobs) % dewp, 0
               close(7)
               print*,jtl,regname(jreg),'DAUWPUNT', obs(jobs) % ilon, obs(jobs) % ilat, obs(jobs) % dewp
            endif
            if (obs(jobs) % rh < 1.0e+38) then
               open(7, file=trim(regname(jreg))//'.RELVOCHT-obs', position='append')
               write(7, fmt = *) jtl, obs(jobs) % rh, 0
               close(7)
               print*,jtl,regname(jreg),'RELVOCHT', obs(jobs) % ilon, obs(jobs) % ilat, obs(jobs) % rh
            endif
            if (obs(jobs) % vis < 1.0e+38) then
               open(7, file=trim(regname(jreg))//'.ZICHT-obs', position='append')
               write(7, fmt = *) jtl, obs(jobs) % vis, 0
               close(7)
               print*,jtl,regname(jreg),'ZICHT', obs(jobs) % ilon, obs(jobs) % ilat, obs(jobs) % vis
            endif
            if (obs(jobs) % cc < 1.0e+38) then
               open(7, file=trim(regname(jreg))//'.BEWOLKING-obs', position='append')
               write(7, fmt = *) jtl, obs(jobs) % cc, 0
               close(7)
               print*,jtl,regname(jreg),'BEWOLKING', obs(jobs) % ilon, obs(jobs) % ilat, obs(jobs) % cc
            endif
            if (obs(jobs) % rr < 1.0e+38) then
               open(7, file=trim(regname(jreg))//'.NEERSLAG-obs', position='append')
               write(7, fmt = *) jtl, obs(jobs) % rr, 0
               close(7)
               print*,jtl,regname(jreg),'NEERSLAG', obs(jobs) % ilon, obs(jobs) % ilat, obs(jobs) % rr
            endif
            if (obs(jobs) % gust < 1.0e+38) then
               open(7, file=trim(regname(jreg))//'.GUST-obs', position='append')
               write(7, fmt = *) jtl, obs(jobs) % gust, 0
               close(7)
               print*,jtl,regname(jreg),'GUST', obs(jobs) % ilon, obs(jobs) % ilat, obs(jobs) % gust
            endif
            if (obs(jobs) % snowdepth < 1.0e+38) then
               open(7, file=trim(regname(jreg))//'.SNOWDEPTH-obs', position='append')
               write(7, fmt = *) jtl, obs(jobs) % snowdepth, 0
               close(7)
               print*,jtl,regname(jreg),'SNOWDEPTH', obs(jobs) % ilon, obs(jobs) % ilat, obs(jobs) % snowdepth
            endif
         endif
      enddo
   enddo

!        Get fields to be plotted.

   do jval = 1, jpval

      ib1 = imdi
      ib1(9:11) = fld(jval) % block1

      ilen = jplen
      ierr = 0
      call GETFD(1, ib1, jpb1, ib2, jpb2, zb2, jpb2, zsup, jpsup, field, ilen, 0, ierr)
      if (ierr < 0) then
         write(0,*)'TIMESERIES - ERROR ON READING FIELD:', ib1(9:11), ' RETURN CODE =', ierr
         call abort
      endif

!        Transfer field values to region array.

      do jreg = 1, nregs
         do jsq1 = 1, nregside(jreg)
            do jsq2 = 1, nregside(jreg)
               values((jsq1 - 1) * nregside(jreg) + jsq2, jreg, jval, jtl/jpfcint) = field(indxreg(jsq1, jsq2, jreg))
            enddo
         enddo
      enddo

   enddo

!        Close the ASIMOF file.

   call ASIMHC(1, ierr)

enddo

!     Determine index in ib1 array for 10 meter wind (and assume that v follows u).

do jval = 1, jpval
   if (all(fld(jval) % block1 == [33, 105, 10])) indxwnd = jval
enddo

!     Convert wind components (u,v) to wind speed and direction (ff,ddd)

do jtl = 0, jpfclen, jpfcint
   do jreg = 1, nregs
      npts = nregside(jreg)**2
      block
         real u(npts), v(npts), dir(npts), ff(npts), ddd(npts)
         dir = regdir(jreg)
         u = values(1:npts,jreg,indxwnd+0,jtl/jpfcint)
         v = values(1:npts,jreg,indxwnd+1,jtl/jpfcint)
         call gb2r2v(u,v,dir,npts)
         ddd = - 90.0 - zrad * atan2(v, u)
         where (ddd < 0.0) ddd = ddd + 360.0
         ff = sqrt(u**2 + v**2)
         values(1:npts,jreg,indxwnd+0,jtl/jpfcint) = ddd
         values(1:npts,jreg,indxwnd+1,jtl/jpfcint) = ff
      end block
   enddo
enddo

!     Switch to the output phase of the program:

do jreg = 1, nregs
   npts = nregside(jreg)**2

!        Determine "Weercijfers" for 18 UTC run.

   if (dtg(9:10) == '18') then

   if (dtg(9:10) == '00') itloff =  7
   if (dtg(9:10) == '06') itloff =  1
   if (dtg(9:10) == '18') itloff = 13

   do jpts = 1, npts
      iwctd(jpts) = 10
      iwctm(jpts) = 10
      iwcda(jpts) = 10
   enddo

!        Cloud cover.

   do jpts = 1, npts
      val = 0.0
      do jtl = itloff, itloff + 11
         val = val + values(jpts,jreg,7,jtl)
      enddo
      if (val > 12.0 * 0.125) iwctd(jpts) = iwctd(jpts) - 1
      if (val > 12.0 * 0.500) iwctd(jpts) = iwctd(jpts) - 1
      if (val > 12.0 * 0.875) iwctd(jpts) = iwctd(jpts) - 1
      val = 0.0
      do jtl = itloff + 24, itloff + 24 + 11
         val = val + values(jpts,jreg,7,jtl)
      enddo
      if (val > 12.0 * 0.125) iwctm(jpts) = iwctm(jpts) - 1
      if (val > 12.0 * 0.500) iwctm(jpts) = iwctm(jpts) - 1
      if (val > 12.0 * 0.875) iwctm(jpts) = iwctm(jpts) - 1
      val = 0.0
      if (jpfclen == 78) then
      do jtl = itloff + 48, itloff + 48 + 11
         val = val + values(jpts,jreg,7,jtl)
      enddo
      if (val > 12.0 * 0.125) iwcda(jpts) = iwcda(jpts) - 1
      if (val > 12.0 * 0.500) iwcda(jpts) = iwcda(jpts) - 1
      if (val > 12.0 * 0.875) iwcda(jpts) = iwcda(jpts) - 1
      endif
   enddo

!        Fog (not implemented yet).

!        Precipitation.

   do jpts = 1, npts
      ival = 0
      prv = values(jpts,jreg,5,itloff-1) + values(jpts,jreg,6,itloff-1)
      do jtl = itloff, itloff + 11
         val = values(jpts,jreg,5,jtl) + values(jpts,jreg,6,jtl)
         if (val - prv > 0.3) ival = ival + 1
         prv = val
      enddo
      if (ival >=  1) iwctd(jpts) = iwctd(jpts) - 1
      if (ival >=  3) iwctd(jpts) = iwctd(jpts) - 1
      if (ival >=  7) iwctd(jpts) = iwctd(jpts) - 1
      if (ival >= 11) iwctd(jpts) = iwctd(jpts) - 1
      ival = 0
      prv = values(jpts,jreg,5,itloff-1+24) + values(jpts,jreg,6,itloff-1+24)
      do jtl = itloff + 24, itloff + 24 + 11
         val = values(jpts,jreg,5,jtl) + values(jpts,jreg,6,jtl)
         if (val - prv > 0.3) ival = ival + 1
         prv = val
      enddo
      if (ival >=  1) iwctm(jpts) = iwctm(jpts) - 1
      if (ival >=  3) iwctm(jpts) = iwctm(jpts) - 1
      if (ival >=  7) iwctm(jpts) = iwctm(jpts) - 1
      if (ival >= 11) iwctm(jpts) = iwctm(jpts) - 1
      if (jpfclen == 78) then
      ival = 0
      prv = values(jpts,jreg,5,itloff-1+24) + values(jpts,jreg,6,itloff-1+24)
      do jtl = itloff + 48, itloff + 48 + 11
         val = values(jpts,jreg,5,jtl) + values(jpts,jreg,6,jtl)
         if (val - prv > 0.3) ival = ival + 1
         prv = val
      enddo
      if (ival >=  1) iwcda(jpts) = iwcda(jpts) - 1
      if (ival >=  3) iwcda(jpts) = iwcda(jpts) - 1
      if (ival >=  7) iwcda(jpts) = iwcda(jpts) - 1
      if (ival >= 11) iwcda(jpts) = iwcda(jpts) - 1
      endif
   enddo

!        Wind.

   do jpts = 1, npts
      ival = 0
      do jtl = itloff, itloff + 9
         if (ival < 1 .and.                            &
             values(jpts,jreg,4,jtl+0) +               &
             values(jpts,jreg,4,jtl+1) +               &
             values(jpts,jreg,4,jtl+2) > 3.0 * 3.4)    &
                ival = 1
         if (ival < 2 .and.                            &
             values(jpts,jreg,4,jtl+0) +               &
             values(jpts,jreg,4,jtl+1) +               &
             values(jpts,jreg,4,jtl+2) > 3.0 * 5.5)    &
                ival = 2
         if (ival < 3 .and.                            &
             values(jpts,jreg,4,jtl+0) +               &
             values(jpts,jreg,4,jtl+1) +               &
             values(jpts,jreg,4,jtl+2) > 3.0 * 10.8)   &
                ival = 3
      enddo
      iwctd(jpts) = iwctd(jpts) - ival
      ival = 0
      do jtl = itloff + 24, itloff + 24 + 9
         if (ival < 1 .and.                            &
             values(jpts,jreg,4,jtl+0) +               &
             values(jpts,jreg,4,jtl+1) +               &
             values(jpts,jreg,4,jtl+2) > 3.0 * 3.4)    &
                ival = 1
         if (ival < 2 .and.                            &
             values(jpts,jreg,4,jtl+0) +               &
             values(jpts,jreg,4,jtl+1) +               &
             values(jpts,jreg,4,jtl+2) > 3.0 * 5.5)    &
                ival = 2
         if (ival < 3 .and.                            &
             values(jpts,jreg,4,jtl+0) +               &
             values(jpts,jreg,4,jtl+1) +               &
             values(jpts,jreg,4,jtl+2) > 3.0 * 10.8)   &
                ival = 3
      enddo
      iwctm(jpts) = iwctm(jpts) - ival
      if (jpfclen == 78) then
      ival = 0
      do jtl = itloff + 48, itloff + 48 + 9
         if (ival < 1 .and.                            &
             values(jpts,jreg,4,jtl+0) +               &
             values(jpts,jreg,4,jtl+1) +               &
             values(jpts,jreg,4,jtl+2) > 3.0 * 3.4)    &
                ival = 1
         if (ival < 2 .and.                            &
             values(jpts,jreg,4,jtl+0) +               &
             values(jpts,jreg,4,jtl+1) +               &
             values(jpts,jreg,4,jtl+2) > 3.0 * 5.5)    &
                ival = 2
         if (ival < 3 .and.                            &
             values(jpts,jreg,4,jtl+0) +               &
             values(jpts,jreg,4,jtl+1) +               &
             values(jpts,jreg,4,jtl+2) > 3.0 * 10.8)   &
                ival = 3
      enddo
      iwcda(jpts) = iwcda(jpts) - ival
      endif
   enddo

!        Total.

   iwctdl = 10
   iwctdh =  0
   iwctml = 10
   iwctmh =  0
   iwcdal = 10
   iwcdah =  0
   do jpts = 1, npts
      iwctd(jpts) = max(iwctd(jpts), 0)
      iwctm(jpts) = max(iwctm(jpts), 0)
      iwcda(jpts) = max(iwcda(jpts), 0)
      if (iwctd(jpts) < iwctdl) iwctdl = iwctd(jpts)
      if (iwctd(jpts) > iwctdh) iwctdh = iwctd(jpts)
      if (iwctm(jpts) < iwctml) iwctml = iwctm(jpts)
      if (iwctm(jpts) > iwctmh) iwctmh = iwctm(jpts)
      if (iwcda(jpts) < iwcdal) iwcdal = iwcda(jpts)
      if (iwcda(jpts) > iwcdah) iwcdah = iwcda(jpts)
   enddo

!        Write the "weercijfers" out.

   open(7, file=trim(regname(jreg))//'.WEERCIJFER')
   if (iwctdl == iwctdh) then
      write(chwctd, '(i2)') iwctdl
   else
      write(chwctd, '(i2,a,i2)') iwctdl,' - ',iwctdh
   endif
   if (iwctml == iwctmh) then
      write(chwctm, '(i2)') iwctml
   else
      write(chwctm, '(i2,a,i2)') iwctml,' - ',iwctmh
   endif
   if (jpfclen == 78) then
   if (iwcdal == iwcdah) then
      write(chwcda, '(i2)') iwcdal
   else
      write(chwcda, '(i2,a,i2)') iwcdal,' - ',iwcdah
   endif
   write(7, '(5a)') chwctd,' , ',chwctm,' , ',chwcda
   else
   write(7, '(3a)') chwctd,' , ',chwctm
   endif
   close(7)

   endif

!        Determine statistics of fields and write them out.

   do jval = 1, jpval

      open(7, file=trim(regname(jreg))//'.'//trim(fld(jval) % name))

      do jtl = 0, jpfclen, jpfcint
      block

         real :: smean, ssdev, smax, smin

         smean = sum(values(1:npts,jreg,jval,jtl/jpfcint)) / real(npts)
         ssdev = sqrt(max(0.0, sum(values(1:npts,jreg,jval,jtl/jpfcint)**2) / real(npts) - smean**2))
         smax  = maxval(values(1:npts,jreg,jval,jtl/jpfcint))
         smin  = minval(values(1:npts,jreg,jval,jtl/jpfcint))

         write(format, fmt=*) '(i4,', npts+4, fld(jval) % format, ')'
         write(7, fmt=format) jtl, smean, ssdev, smin, smax, values(1:npts,jreg,jval,jtl/jpfcint)

     end block
     end do

     close(7)

  enddo

enddo

end
!
!**   Subroutine getnnb - Get nearest neighbour of given lat / lon
!**                       point on the grid.
!
subroutine getnnb(  &
           plon,    & ! Array of longitudes of requested points
           plat,    & ! Array of latitudes of requested points
           pdir,    & ! Array of angles to rotate vectors by
           pwest,   & ! Western boundary of area
           psouth,  & ! Southern boundary of area
           pdlon,   & ! Grid distance in longitudinal direction
           pdlat,   & ! Grid distance in latitudinal direction
           plonp,   & ! Longitude of shifted pole
           platp,   & ! Latitude of shifted pole
           klon,    & ! Number of grid points in long. dir.
           klat,    & ! Number of grid points in lat. dir.
           kpos,    & ! Array to get indices of nearest nghbrs
           knts     & ! Number of valid points found
   )
!
implicit   none
!
! ... Declare dummy arguments
!
integer    knts
integer    klon, klat
integer    kpos(knts)
real       plon(knts), plat(knts)
real       pdir(knts)
real       pwest, psouth
real       pdlon, pdlat
real       plonp, platp
!
! ... Local constants
!
integer    jlb2
parameter  (jlb2 = 400)
!<A NAME="n10">
! ... Local variables
!
integer    i, j       ! Some counters
integer    ilon, ilat ! Longitude / latitude indices
integer    ib2r(jlb2) ! Block 2 of GRIB definition regular lat/lon
integer    ib2s(jlb2) ! Block 2 of GRIB definition shifted pole
integer    iposm      ! Position of max/min fraction of land
real       zlon, zlat ! Longitude/latitude in shifted pole coords.
real       zdir       ! Dummy for wind rotation on 2nd transform
real       zb2r(jlb2) ! Block 2 (real part), regular lat/lon
real       zb2s(jlb2) ! Block 2 (real part), shifted pole.
!
ib2r( 6) = 0
ib2r(17) = 128 + 8
ib2s( 6) = 10
ib2s(17) = 128 + 8
ib2s(33) = nint(platp * 1000.)
ib2s(36) = nint(plonp * 1000.)
ib2s(39) = 1
zb2s(ib2s(39)) = 0.
!<A NAME="n20">
! ... Loop over all given lat / lon points
!
do i = 1, knts
!
!        Check if valid (i.e. non-null)
!
   if (plon(i) == 0.0 .and. plat(i) == 0.0) then
      plat(i)=99.
      exit
   endif
!
!        Convert to shifted pole coordinates
!
   call gb2gca(          &
           plon(i),      & ! Old longitude
           plat(i),      & ! Old latitude
           zlon,         & ! New longitude
           zlat,         & ! New latitude
           zdir,         & ! Dummy variable for rotation of wind
           1,            & ! Number of points to convert
           ib2r,         & ! GRIB block 2 regular lat/lon
           zb2r,         & ! GRIB block 2 real part
           ib2s,         & ! GRIB block 2 shifted pole
           zb2s          & ! GRIB block 2 real part
   )
!<A NAME="n30">
!        Determine longitude offset in area of point to SW
!
   ilon = max(1,min(klon-1,int((zlon - pwest ) / pdlon) + 1))
   zlon = zlon - pwest  - (ilon-1) * pdlon
!
!        Determine latitude offset in area of point to SW
!
   ilat = max(1,min(klat-1,int((zlat - psouth) / pdlat) + 1))
   zlat = zlat - psouth - (ilat-1) * pdlat
!
!        Calculate index in arrays
!
   kpos(i) = ilon + (ilat - 1) * klon
!
   ilat = (kpos(i)-1)/klon
   ilon = kpos(i) - ilat*klon
!<A NAME="n70">
!        Calculate real longitude of nearest neighbour
!
   zlon = pwest + (ilon - 1) * pdlon
!
!        Calculate real latitude of nearest neighbour
!
   zlat = psouth + ilat      * pdlat
!
!        Transform back to normal lat / lon coordinates
!
   call gb2gca(          &
           zlon,         & ! Longitude in shifted pole coords.
           zlat,         & ! Latitude in shifted pole coords.
           plon(i),      & ! Longitude in normal coords.
           plat(i),      & ! Latitude in normal coords.
           pdir(i),      & ! Angle of rotation
           1,            & ! Number of points to convert
           ib2s,         & ! GRIB block 2 shifted pole
           zb2s,         & ! GRIB block 2 real part
           ib2r,         & ! GRIB block 2 regular lat/lon
           zb2r          & ! GRIB block 2 real part
   )
!
enddo
!<A NAME="n80">
! ... Remove all invalid points ...
!
j = 0
do i = 1, knts
   if (plat(i) < 90.1) then
      j = j + 1
      plon(j) = plon(i)
      plat(j) = plat(i)
      pdir(j) = pdir(i)
      kpos(j) = kpos(i)
   endif
enddo
knts = j
!
return
end