!
      module m_interp_NR
!
! Module for reading binary files of NR required fields 
! and performing interpolations for simulation of observations for an OSSE.
!
! This present version treats a rawindsonde sounding as though it were at 
! a single geographic location (the station location) and time. Also, there 
! is no re-determination of significant levels: whatever levels were recorded
! for the original (real) observation are used here, irrespective of their 
! "significance" in the NR profile. 
!
! See the main program sim_obs_conv.f90 for an example of use.
!
! Initial code by Ronald Errico January 2008
!   June-09-2008:    R.E.-- vertical log(p) interpolation
!   August-08-2008:  R.E.-- reconfigure to use single 3D field array
!   Sept-17-2008:    R.E.-- compute hydrostatic height if requested
!
      use m_kinds
      use m_relhum
!
      implicit none
!
      private
      public :: clean_m_interp
      public :: find_name
      public :: setup_m_interp
      public :: setup_NR_fields
      public :: interp_NR_fields 
      public :: interp_fields_for_rad
!
! Parameters for transforming t or q variables: KEEP lq2rh as FALSE 
      logical, parameter :: lq2rh=.false. ! interpolate using rh rather than q
      integer, parameter :: kq2rh=15      ! change q to rh for levs k >= this
!
! Global variables
      integer, save :: k_iter     ! numbers of factors of 2 (+1) in nlevs1
      integer, save :: nfields2d  ! number of 2d fields read from NR
      integer, save :: nfields3d  ! number of 3d fields read from NR
      integer, save :: nfields    ! nfield2sd+nfields3d
      integer, save :: nfdim      ! number of points on horz NR grid with poles
      integer, save :: nlevs1     ! number of levels on NR grid, incl. surface
      integer, save :: nlats2     ! number of lats on grid with poles included
      integer, save :: ntimes     ! number of NR times for each obs period
      integer, save :: nft, nfp, nfq, nfo, nfu, nfz, nfc ! indexes for some fields
      logical, save :: lqfield    ! true if q field is used
      logical, save :: ltfield    ! true if t field is used
      logical, save :: luvfield   ! true if u and v fields are used
!
! Global allocatable arrays
      integer, allocatable, save :: nlons2(:)
      integer, allocatable, save :: nlonsP(:)
! p=ak+bk*ps define interfaces (:,1) and data levels (:,2), including surface
      real(rkind1), allocatable, save :: ak(:,:) 
      real(rkind1), allocatable, save :: bk(:,:) 
      real(rkind2), allocatable, save :: f2d(:,:,:,:)   ! 2-D fields all times
      real(rkind2), allocatable, save :: f3d(:,:,:,:)   ! 3-D fields all times
      real(rkind1), allocatable, save :: glats2(:) ! lats (deg), incl. poles 
      real(rkind1), allocatable, save :: time_files(:)  ! times of the NR files
!
! List of field names on binary NR data files:
!
! 3D fields:  
!      u:  'uwnd'   unit: m/s
!      v:  'vwnd'   unit: m/s
!      T:  'temp'   unit: K
!      q:  'sphu'   unit: Specific humidity [kg/ kg]
!      O3: 'ozon'   unit: Ozone mass mixing ratio [kg/ kg]
!      c:  'clwc'   cloud liquid water: [kg water/ kg air]
!
! 2D ps and zs fields
!      ps: 'pres'   unit: Pa
!      zs: 'zsfc'   unit: m
!  
! 2D fields (excl. z and ps)
!          'ismk'   surface sea-ice cover (0,1)
!          'rain'   surface stratiform precip. (m)
!          'conp'   surface convective precip. (m)
!          'u10m'   uwind at 10 meter  m/s
!          'v10m'   vwind at 10 meter  m/s
!          'tp2m'   surface temperature at 2 meter    K
!          'dt2m'   dew point temperature at 2 meter   K
!          'almk'   sea (Aqua) land mask (0,1)
!          'lcld'   low cloud cover fraction (0,1)
!          'mcld'   middle  cloud cover fraction (0,1)
!          'hcld'   high cloud cover fraction (0,1)
!          'sktp'   surface skin temperature  (k)
!
      Contains
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      Subroutine clean_m_interp 
!
      implicit none
!
      deallocate (nlons2,nlonsP,ak,bk,glats2,time_files)
      deallocate (f2d,f3d)
!
      print *,' '
      print *,' Grid arrays and fields deallocated'
!
      end subroutine clean_m_interp 
!
! 
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine compute_p_20m (ltest,jtime1,h_index,obs_time, &
                                   ps,h_weights,obs_plev)
!
!  Compute p for observation specified as 20m.
!
      use m_parameters
      implicit none
! input
      logical :: ltest  ! true if extra info printed to test routine
      integer :: jtime1
      integer :: h_index(2,2)
      real(rkind1) :: h_weights(2,3)       ! weights for horiz. interp.
      real(rkind1) :: obs_time             ! time of observation
      real(rkind1) :: ps                   ! surface pressure
! output
      real(rkind1) :: obs_plev             ! pressure level of obs.
! parameter
      real(rkind1), parameter :: obs_zlev=20. ! observation at 20 m.
! local
      integer :: j, jt     ! indexes for time
      integer :: n, nf     ! indexes for field numbers
      integer :: nlevs1m1, nlevs1m2  ! indexes for 2 levels above surface
      real(rkind1) :: dz
      real(rkind1) :: f3d_interp(2,3,2) ! (levels,times,fields)
      real(rkind1) :: rfac              ! R/g
      real(rkind1) :: rt                ! dz/d(ln(p))
      real(rkind1) :: pa                ! p at 1st interface above sfc
      real(rkind1) :: za                ! z at 1st interface above sfc
!  
! Check if nft and nfq exist
      if ((nft==0) .or. (nfq==0)) then
        print *,' '
        print *,'Error in routine compute_z_obs:'
        print *,'nft=',nft,' nfq=',nfq,' must not be 0' 
        stop
      endif
!
! Check that q is specific humidity at this point and not rel. humidity
      if (lq2rh) then
        print *,' '
        print *,'Error in routine compute_z_obs:'
        print *,'lq2rh=',lq2rh,' is invalid since z needs q, not rh.' 
        stop
      endif   
!
! Interpolate T and q fields at first two atmos levels above surface
      nlevs1m2=nlevs1-2
      do n=1,2
        if (n==1) then
          nf=nft  ! index for T field
        else
          nf=nfq  ! index for q field
        endif
        do j=1,2
          jt=jtime1+j-1
          call interp_horz (nfdim,2,h_index,h_weights, &
                            f3d(1,nlevs1m2,nf,jt),f3d_interp(1,j,n))  
        enddo
        call interp_time (2,obs_time,time_files(jtime1:jtime1+1), &
                          f3d_interp(1:2,:,n))
      enddo
!
! Calculate height above surface for first interface (klev=2 in f3d_interp) 
      nlevs1m1=nlevs1-1
      rfac=Rdrygas/gravity
      pa=ak(nlevs1m1,1)+bk(nlevs1m1,1)*ps   ! value of 1st interface
      rt=rfac*f3d_interp(2,3,1)*(one_k1+ratio4R*f3d_interp(2,3,2))
      za=rt*log(ps/pa)
!
      dz=obs_zlev-za
      if (dz > zero_k1) then
        rt=rfac*f3d_interp(1,3,1)*(one_k1+ratio4R*f3d_interp(1,3,2))
      endif
      obs_plev=pa*exp(-dz/rt)
!
! Code for printing information for testing purposes follows
      if (ltest) then
        print *,' '
        print *,'TEST of compute_p_20m'
        print *,'  T=',f3d_interp(:,3,1),'  q=',f3d_interp(:,3,2)
        print *,' obs_zlev,za,rt=',obs_zlev,za,rt
        print *,' ps,pa,obs_plev=',ps,pa,obs_plev
       endif
!
      end subroutine compute_p_20m
!
! 
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine compute_z_obs (ltest,obs_nlevs,jtime1,k_above_z,    &
                     h_index,ps,zs,obs_time,special_value,h_weights, &
                     obs_plevs,obs_zlevs,c_type)
!
! Compute hydrostatic heights of observations
!
! Assume that k_above_z(i) > nlevs1 if obs i is not in range of NR 
!   atmosphere pressures in that column. 
!
      implicit none
!
! input
      logical :: ltest  ! true if extra info printed to test routine
      character(len=*) :: c_type  ! 'sing' or 'mult'
      integer :: obs_nlevs
      integer :: jtime1
      integer :: k_above_z(obs_nlevs)
      integer :: h_index(2,2)
      real(rkind1) :: h_weights(2,3)       ! weights for horiz. interp.
      real(rkind1) :: obs_plevs(obs_nlevs) ! pressure levels of observations
      real(rkind1) :: obs_time             ! time of observation
      real(rkind1) :: ps      ! surface pressure
      real(rkind1) :: special_value  ! value for missing data
      real(rkind1) :: zs      ! surface height
! output
      real(rkind1) :: obs_zlevs(obs_nlevs) ! heights of observations
! local
      integer :: n, nf     ! indexes for field numbers
      integer :: k, ka, kb ! indexes for vertical levels
      integer :: klevs, klevsm1 ! number of levels to process
      integer :: j, jt  ! indexes for time
      real(rkind1) :: f3d_interp(nlevs1,3,2) 
      real(rkind1) :: p_2(2)  ! pressures within and at bottom of a layer
      real(rkind1) :: z_grid(nlevs1)  ! z at required model interfaces
      real(rkind1) :: z_2(2)  ! heights corresponding to p_2
!  
! Check if nft and nfq exist
      if ((nft==0) .or. (nfq==0)) then
        print *,' '
        print *,'Error in routine compute_z_obs:'
        print *,'nft=',nft,' nfq=',nfq,' must not be 0' 
        stop
      endif
!
! Check that q is specific humidity at this point and not rel. humidity
      if (lq2rh) then
        print *,' '
        print *,'Error in routine compute_z_obs:'
        print *,'lq2rh=',lq2rh,' is invalid since z needs q, not rh.' 
        stop
      endif   
!
      if (c_type == 'sing') then
        ka=k_above_z(1)
        if (ka < nlevs1) then  ! obs in range of nature run fields 
          klevs=nlevs1-ka+1
!
! Interpolate T and q fields between the surface and required level
          do n=1,2
            if (n==1) then
              nf=nft  ! index for T field
            else
              nf=nfq  ! index for q field
            endif
!
            do j=1,2
              jt=jtime1+j-1
              call interp_horz (nfdim,klevs,h_index,h_weights, &
                                f3d(1,ka,nf,jt),f3d_interp(ka,j,n))  
            enddo
            call interp_time (klevs,obs_time,time_files(jtime1:jtime1+1), &
                              f3d_interp(ka:nlevs1,:,n))
          enddo
!
! Compute z at the required observation level
          kb=ka+1 
          klevsm1=klevs-1
          z_grid(nlevs1)=zs 
          if ( kb < nlevs1 ) then
            call hydros (klevsm1,ps,ak(kb,1),bk(kb,1),f3d_interp(kb,3,1), &
                         f3d_interp(kb,3,2),z_grid(kb))
          endif
          p_2(1)=obs_plevs(1)
          p_2(2)=ak(kb,1)+bk(kb,1)*ps
          z_2(2)=z_grid(kb)
          call hydros (2,zero_k1,p_2,p_2,f3d_interp(ka,3,1), &
                     f3d_interp(ka,3,2),z_2)        
          obs_zlevs(1)=z_2(1)
!
        else ! obs below ground or above data top
           obs_zlevs(1)=special_value
        endif 
!     
      elseif (c_type == 'mult') then
!
! Find lowest p in obs dataset 
        ka=minval(k_above_z(1:obs_nlevs))
        if (ka < nlevs1) then  ! some obs in range of nature run fields 
          klevs=nlevs1-ka+1
!
! Interpolate T and q fields between the surface and top required level
          do n=1,2
            if (n==1) then
              nf=nft   ! index for T field
            else
              nf=nfq   ! index for q field
            endif
!
            do j=1,2
              jt=jtime1+j-1
              call interp_horz (nfdim,klevs,h_index,h_weights, &
                                f3d(1,ka,nf,jt),f3d_interp(ka,j,n))  
            enddo
            call interp_time (klevs,obs_time,time_files(jtime1:jtime1+1), &
                              f3d_interp(ka:nlevs1,:,n))
          enddo
!
! Compute z at all model levels up to the highest required
          kb=ka+1 
          klevsm1=klevs-1
          z_grid(nlevs1)=zs 
          if ( kb < nlevs1 ) then
            call hydros (klevsm1,ps,ak(kb,1),bk(kb,1),f3d_interp(kb,3,1), &
                       f3d_interp(kb,3,2),z_grid(kb))
          endif
!
! Compute z at all observation levels
          do k=1,obs_nlevs
            kb=k_above_z(k)+1          
            if (kb <= nlevs1) then  ! checks if this particular obs in atmos 
              p_2(1)=obs_plevs(k)
              p_2(2)=ak(kb,1)+bk(kb,1)*ps
              z_2(2)=z_grid(kb)
              call hydros (2,zero_k1,p_2,p_2,f3d_interp(ka,3,1), &
                           f3d_interp(ka,3,2),z_2)        
              obs_zlevs(k)=z_2(1)          ! obs with NR atmosphere
            else
               obs_zlevs(k)=special_value   ! obs outside range of NR data
            endif 
          enddo
!
        else ! problem with data: no valid k_z_above
           obs_zlevs(1:obs_nlevs)=special_value
        endif 
!
      endif   ! check on c_type
!
! Code for printing information for testing purposes follows
      if (ltest) then
        print *,' '
        print *,'TEST of compute_z_obs: k, k above, p obs, z obs'
        do k=1,obs_nlevs
          print ('(2i4,2e15.6)'),k,k_above_z(k),obs_plevs(k),obs_zlevs(k)
        enddo
        print *,'TEST of compute_z_obs: hydros applied to NR grid:'
        do k=nlevs1,ka-1,-1
          p_2(1)=ak(k,1)+bk(k,1)*ps
          p_2(2)=ak(k-1,1)+bk(k-1,1)*ps
          z_2(1)=alog(p_2(1)/p_2(2))
          z_2(2)=f3d_interp(k-1,3,1)*(1.+.622*f3d_interp(k-1,3,2))
          z_grid(1)=z_2(1)*z_2(2)*287.04/9.80616
          print ('(i4,1p8e14.4)'),k,p_2,z_2,z_grid(1),z_grid(k), &
               f3d_interp(k-1,3,1),f3d_interp(k-1,3,2)
        enddo
      endif    ! end of print test information block
!
      end subroutine compute_z_obs
!
! 
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine find_name (nfieldsx,f_namesx,f_namex,nfx,lstopx)
!
! Find the index (nf) corresponding to the first occurance of the name 
! of the field (f_name) within a chacter string (f_names). 
! The name is a single character.
!
! input
      logical          :: lstopx    ! stop if name not found in list
      integer          :: nfieldsx
      character(len=*) :: f_namesx(nfields)
      character(len=*) :: f_namex
!
! output
      integer :: nfx
!
! local
      integer :: n
!
      nfx=0
      do n=1,nfieldsx
        if ( (nfx == 0) .and. (f_namesx(n) == f_namex) ) then
          nfx=n       
        endif
      enddo
!
      if (lstopx .and. (nfx==0)) then
        print *,' '
        print *,' requested field name not found:'
        print *,'    requested name = ',f_namex
        print *,'    possible names = '
        print ('(10(1x,a4))'),f_namesx(1:nfields)
        stop
      endif
! 
      end subroutine find_name
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      Subroutine setup_m_interp (nfields2d_in,nfields3d_in,nlevs,ltest_all)
!
      implicit none
!
      logical :: ltest_all
      integer :: nfields2d_in
      integer :: nfields3d_in
      integer :: nlevs        ! number of data levels excluding the surface
!
      integer :: nlats        ! number of latitudes excluding the poles
!
      if (ltest_all) then     ! set low resolution for testing
         nfdim=44        
         nlevs=5         
         nlats=6         
      else                    ! set to ECMWF NR resolution 
         nfdim=348528+36      ! # gridpoints on each horiz surface, incl poles
         nlevs=91             ! # levels, excluding the surface
         nlats=512            ! number of latitudes
      endif
      nlats2=nlats+2
      nlevs1=nlevs+1
!
! Set times of NR data files that will reside in memory together 
      ntimes=3             ! number of file times in simulation period
      allocate (time_files(ntimes))
      time_files(1)=-3.
      time_files(2)= 0.    ! central time of simulation period
      time_files(3)= 3.
!
! number of iterations required to find NR grid pressures that bound
! pressure of observation 
      k_iter=int(log(real(nlevs1))/log(2.)) + 1
!
      nfields2d=nfields2d_in
      nfields3d=nfields3d_in
      nfields=nfields2d+nfields3d
!
      allocate (nlons2(nlats2))
      allocate (nlonsP(nlats2))
      allocate (ak(nlevs1,2)) 
      allocate (bk(nlevs1,2)) 
      allocate (f2d(nfdim,     1,nfields2d,ntimes)) 
      allocate (f3d(nfdim,nlevs1,nfields3d,ntimes)) 
      allocate (glats2(nlats2))
!
      print *,' '
      print ('(2a)'),' Setup_m_interp for nlevs1, nlats2, ', &
                   ' nfdim, nfields2d, nfields3d'
      print ('(18x,3i8,2i11)'),nlevs1,nlats2,nfdim,nfields2d,nfields3d
!
      end subroutine setup_m_interp 
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      Subroutine setup_NR_fields (f_names,ltest)
!
      implicit none
!
      character(len=*) :: f_names(nfields)
      logical          :: ltest
!
      integer, parameter :: nlons_pole=40 ! >= number of longitudes at pole
      real(rkind1)              :: pcoslon(nlons_pole) ! cosines and sines of 
      real(rkind1)              :: psinlon(nlons_pole) ! longitudes at 1st lat
!
! File information
      integer, parameter :: file_grid_unit=10
      integer, parameter :: file_data_unit=12
      integer, parameter :: file_data_len=10
      character(len=12), parameter :: file_grid_name='ossegrid.txt'
      character(len=file_data_len), parameter :: file_data_2D='surf_NR_01'
      character(len=file_data_len), parameter :: file_data_ps='pres_NR_01'
      character(len=file_data_len), parameter :: file_data_3D='udat_NR_01'
!
      integer :: nf
      integer :: nfx
      integer :: k
      integer :: nt

!
! Set logical variables based on field names to be used
      print *,' '
      print ('(8(1x,a))'),' f_names=',f_names
      lqfield=.false.
      ltfield=.false.
      luvfield=.false.
!
      call find_name (nfields2d,f_names(1:nfields2d),'pres',nfp,.true.)
      call find_name (nfields2d,f_names(1:nfields2d),'zsfc',nfz,.false.)
      call find_name (nfields3d,f_names(1+nfields2d:nfields),'uwnd',nfu,.false.)
      call find_name (nfields3d,f_names(1+nfields2d:nfields),'temp',nft,.false.)
      call find_name (nfields3d,f_names(1+nfields2d:nfields),'sphu',nfq,.false.)
      call find_name (nfields3d,f_names(1+nfields2d:nfields),'ozon',nfo,.false.)
      call find_name (nfields3d,f_names(1+nfields2d:nfields),'clwc',nfc,.false.)
      if (nfu > 0) then 
        luvfield=.true.  ! true if u and v fields are to be used
        call find_name (nfields3d,f_names(1+nfields2d:nfields),'vwnd',nfx,.false.)
        if (nfx /= nfu+1) then
          print *,'Problem with names:',nfu,'=nfu /= nfv=',nfx
          stop
        endif
      endif
      if (nft > 0) then 
        ltfield=.true.  ! true if t field is to be used
      endif
      if (nfq > 0) then 
        lqfield=.true.  ! true if q field is to be used
        if (nfq /= nft+1) then
          print *,'Problem with names:',nft,'=nft+1 /= nfq=',nfq
          stop
        endif
      endif
!
! Set grid information: nlons2, nlonsP, glats2, ak, bk
      if (ltest) then
        call set_gridinfo_test (nfdim,nlats2,nlevs1,nlons2,nlonsP, &
                                glats2,ak,bk)   
      else
        call read_gridinfo (nfdim,nlats2,nlevs1,file_grid_unit, &
                            file_grid_name,nlons2,nlonsP,glats2,ak,bk) 
      endif
      print *,' '
      print *,' Grid information set'
!
      if (ltest) then
        print *,'print k,nlons2(k),nlonsp(k),glats(k):'
        do k=1,nlats2
          print ('(3i8,f8.2)'),k,nlons2(k),nlonsp(k),glats2(k)
        enddo
        print *,'k,ak,bk,pk at interface levels, for computing heights'
        do k=1,nlevs1
          print ('(i3,3f15.5)'),k,ak(k,1),bk(k,1),ak(k,1)+bk(k,1)*1.e5
        enddo
        print *,'k,ak,bk,pk at data levels (includes surface)'
        do k=1,nlevs1
          print ('(i3,3f15.5)'),k,ak(k,2),bk(k,2),ak(k,2)+bk(k,2)*1.e5
        enddo
      endif  ! test of ltest
!
! Set sines and cosines for wind extrapolation to poles.
      if (luvfield) then
        if (nlons_pole < nlons2(1)) then 
          print *,' '
          print *,' nlons_pole=',nlons_pole,' < nlons2(1) =',nlons2(1)
          print *,' need larger value for nlons_pole'
          stop 
        elseif (nlons2(1) == nlons2(nlats2)) then 
          call table_sincos (nlons2(1),pcoslon(1:nlons2(1)), &
                             psinlon(1:nlons2(1)))
        else
          print *,' '
          print *,' nlons2(1)=',nlons2(1),' /= nlons2(nlats2) =', &
                    nlons2(nlats2)
          print *,' need 2 different arrays for cos and sin at poles'
          stop 
        endif
      endif  
!
! Set table of saturation vapor pressures if required
      if (lqfield) then
        call relhum_setup
      endif
!      
! Read 3 times of data or fill arrays with low resolution fake data 
      if (ltest) then 
        call read_NR_test (nfdim,nlevs1,nlats2,ntimes,nfields2d,nfields3d, &
                           nfields,nlons2,nlonsP,f_names,f2d,f3d)
      else
        do nt=1,ntimes
          call read_NR (nfdim,nlevs1,nlats2,1,nfields2d,nfields3d,nfields, &
                        file_data_unit,file_data_len,nlons2(:),nlonsP(:),  &
                        f_names,file_data_2D,file_data_ps,file_data_3D,    &
                        f2d(1,1,1,nt),f3d(1,1,1,nt),nt)
        enddo
      endif
!
! Fill pole points
      do nf=1,nfields
        nfx=nf-nfields2d
        if (f_names(nf) == 'uwnd') then 
          if (nf > nfields2d) then
            do nt=1,ntimes
              call pole_uv (nfdim,nlevs1,nlats2,nlons2,nlonsP,nlons_pole, &
                            pcoslon,psinlon,f3d(1,1,nfx,nt))
            enddo
          else
            do nt=1,ntimes
              call pole_uv (nfdim,    1,nlats2,nlons2,nlonsP,nlons_pole, &
                            pcoslon,psinlon,f2d(:,:,nf:nf+1,nt))
            enddo
          endif
        elseif (f_names(nf) /= 'vwnd') then 
          if (nf > nfields2d) then
            do nt=1,ntimes
              call pole_tq (nfdim,nlevs1,nlats2,nlons2,nlonsP,f3d(1,1,nfx,nt))
            enddo
          else
            do nt=1,ntimes
              call pole_tq (nfdim,1,nlats2,nlons2,nlonsP,f2d(1,1,nf,nt))
            enddo
          endif
        endif
      enddo ! loop over fields
      print *,' Pole values set'
!
      if (ltest .and. nfdim<100) then
        do nt=1,ntimes
          call write_field_test (nfdim,1,nlats2,nlons2,nlonsP,nt, &
                                 f2d(:,:,nfp,nt),f_names(nfp:nfp)) 
          call write_field_test (nfdim,nlevs1,nlats2,nlons2,nlonsP,nt, &
                               f3d(:,:,1,nt),f_names(nfields2d+1:nfields2d+1)) 
          call write_field_test (nfdim,nlevs1,nlats2,nlons2,nlonsP,nt, &
                               f3d(:,:,2,nt),f_names(nfields2d+2:nfields2d+2)) 
        enddo  ! loop over nt
      endif    ! test on ltest
!
! If 3-d field of q is present, transform dew point temperature to q at surface
      if (lqfield) then
        do nt=1,ntimes
          call td2q (nfdim,f2d(1,1,nfp,nt),f3d(1,nlevs1,nfq,nt))
        enddo
        print *,' td converted to q at surface'
      endif 
!
! Convert q to rh: All interpolations of q done in terms of rh
      if (lqfield .and. lq2rh) then 
        do nt=1,ntimes
          call q2rh (kq2rh,nfdim,nlevs1,ak(:,2),bk(:,2),                &
             f3d(:,:,nft,nt),f2d(:,1,nfp,nt),f3d(:,:,nfq,nt))
        enddo
        print *,' q converted to rh'
      endif  
!
! For all 3d fields other than T, q, u, v, copy lowest level values to surface.
      do nf=1,nfields3d
        if ( (nf /= nft) .and. (nf /= nfq) .and. ( .not. &
             (luvfield .and. ((nf == nfu) .or. (nf == nfu+1))) )) then
          do nt=1,ntimes
            f3d(:,nlevs1,nf,nt)=f3d(:,nlevs1-1,nf,nt)
          enddo
        endif
      enddo
!
      print *,'Setup of NR fields completed'
!
      end subroutine setup_NR_fields
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine interp_NR_fields (obs_ndim1,obs_ndim2,obs_ndim3,        &
                                   obs_ndim4,obs_nlevs,obs_nfields,      &
                                   nobs,n_errors,ierrors,bmiss,l_comp_z, &
                                   obs_info8,obs_levels8,obs_values8,ltest)
!
      implicit none
!
      logical :: l_comp_z    ! flag to compute height of obs
      integer :: obs_ndim1, obs_ndim2, obs_ndim3, obs_ndim4
      integer :: obs_nlevs   ! number of pressure levels for this obs set
      integer :: obs_nfields ! number of free-atmosphere fields observed
      integer :: nobs        ! observation set counter 
      integer :: n_errors    ! number of error counters in ierrors array
      integer :: ierrors(n_errors)  ! error counter
      real(rkind3) :: bmiss         ! flag for missing observation value
      real(rkind3) :: obs_info8(obs_ndim4)  ! header info for obs. set
      real(rkind3) :: obs_levels8(obs_ndim1,obs_ndim3) ! pressure levels
      real(rkind3) :: obs_values8(obs_ndim1,obs_ndim2) ! obs fields
      logical :: ltest              ! if true, then print lots of stuff
!
      logical :: lskip0             
      logical :: lskip1
      logical :: lskip2
      logical :: lsurface ! true if observations are of surface only
!         
      integer :: h_index(2,2)
      integer :: jtime1
      integer :: k
      integer :: k_above
      integer :: k_below
      integer :: k_above_z(obs_ndim1) ! stored values for computing z 
      integer :: nf
      integer :: np
      integer :: nt
      integer :: nt12 
      real(rkind1) :: h_weights(2,3)
      real(rkind1) :: obs_lat               ! latitude of observation 
      real(rkind1) :: obs_lon               ! longitude of observation 
      real(rkind1) :: obs_height            ! height of surface report
      real(rkind1) :: obs_plev              ! pressure of 1 observation 
      real(rkind1) :: obs_plevs(obs_ndim1)  ! pressures of multiple obs.
      real(rkind1) :: obs_zlevs(obs_ndim1)  ! heights of multiple obs.
      real(rkind1) :: obs_time              ! time of observation 
      real(rkind1) :: obs_type              ! flag for type of observation  
      real(rkind1) :: obs_info(obs_ndim4)
      real(rkind1) :: obs_levels(obs_ndim1,obs_ndim3)
      real(rkind1) :: obs_values(obs_ndim1,obs_ndim2)
      real(rkind1) :: p_above               ! p above in box bounding data 
      real(rkind1) :: p_below               ! p below in box bounding data
      real(rkind1) :: special_value         ! flag for unusable observations  

      real(rkind1), allocatable :: f2d_interp(:,:,:) ! horiz interp. 2d fields
      real(rkind1), allocatable :: f3d_interp(:,:,:) ! horiz interp. 3d fields
!
      allocate (f2d_interp(1,3,nfields2d))
      allocate (f3d_interp(nlevs1,3,2))
!
      obs_info=obs_info8
      obs_levels=obs_levels8
      obs_values=obs_values8
!
! Set or initialize some variables
      special_value=bmiss
!
! Set some observation information
      obs_lon=   obs_info(2)
      obs_lat=   obs_info(3)
      obs_time=  obs_info(4)
      obs_type=  obs_info(5)
      obs_height=obs_info(6)
      obs_plevs(1:obs_nlevs)=obs_levels(1:obs_nlevs,1)
!
      ierrors(8)=ierrors(8)+obs_nfields*obs_nlevs ! tot number of obs values 
!
      if (obs_nlevs<1) then
        ierrors(10)=ierrors(10)+1
        lskip0=.true.
      else 
        ierrors(11)=ierrors(11)+1          
        lskip0=.false.
      endif
!
! Check appropriateness of data time and find time index 
      lskip1=.false.
      if (lskip0) then
        lskip1=.true. 
      else 
        call index_time (ntimes,obs_time,time_files(:),jtime1)
        if (ltest) then
          print *,' '
          print *,' '
          print *,' TEST OBSERVATION NUMBER ',nobs
          print *,'obs_time,time_files(:),jtime1=', &
                   obs_time,time_files(:),jtime1
        endif
        if (jtime1 == -1) then
          ierrors(1)=ierrors(1)+1
          lskip1=.true.
        endif
        if (jtime1 == -2) then
          ierrors(2)=ierrors(2)+1
          lskip1=.true.
        endif
!       
! Check lat and lon and find grid location indexes
        if ( (obs_lon < -180.) .or. (obs_lon > 360.) ) then
          ierrors(3)=ierrors(3)+1
          lskip1=.true.
        endif
        if ( (obs_lat < -90.) .or. (obs_lat > 90.) ) then
          ierrors(4)=ierrors(4)+1
          lskip1=.true.
        endif
!
      endif   ! test on lskip0
!
! If acceptable location then continue with interpolation
!
      if (.not.lskip1) then  
!
!
! Initialize obs. heights to missing values if they are to be computed here
        if (l_comp_z) then
          obs_zlevs(1:obs_nlevs)=special_value
        endif
!
! Check flag that indicates pressure to be set to surface
        if ( (obs_nlevs==1) .and. (nint(obs_plevs(1))<zero_k1) ) then
          lsurface=.true.
        else
          lsurface=.false.
        endif
!
! Compute grid box location that surrounds observation location
        call index_horz (obs_lat,obs_lon,nlats2,nlons2,nlonsP, &
                         glats2,h_index,h_weights,ltest)
!  
! Compute horizontally interpolated values of 2d fields at 2 times
        nt12=0
        do nt=jtime1,jtime1+1
          nt12=nt12+1
          do nf=1,nfields2d
            call interp_horz (nfdim,1,h_index,h_weights,    &
                      f2d(1,1,nf,nt),f2d_interp(1,nt12,nf))        
          enddo 
        enddo 
! 
! Compute temporally interpolated values of 2d fields 
        do nf=1,nfields2d
          call interp_time (1,obs_time,time_files(jtime1:jtime1+1), &
                            f2d_interp(:,:,nf))          
        enddo
!
! If flag on obervation indicates p corresponds to surface pressure, then
! set value of obs_plevs to the interpolated surface pressure and if
! surface height field is present, copy that to report elevation 
        do k=1,obs_nlevs
          if (obs_plevs(k) < zero_k1 ) then
            obs_plevs(k)=f2d_interp(1,3,nfp) 
            obs_levels(k,1)=obs_plevs(k)
            if (nfz > 0) then
              obs_info(6)=f2d_interp(1,3,nfz) 
            endif
          endif
        enddo
!
        if (ltest) then
          print *,' nf,f2d_interp(1,1:3,nf):'
          do nf=1,nfields2d
            print ('(i2,3e17.7)'),nf,f2d_interp(1,1:3,nf)
          enddo 
        endif  ! test on ltest
!
! Interpolation for observations of the surface.
! For T and q, use 2m values from nature run. For wind, use 10m values.
! Above, p of read observation is replace by value in nature run.
! Above, station elevation is replaced by topographic height from nature run.
! 
        if (lsurface) then
!
          f3d_interp(:,:,:)=special_value
          do nf=1,2
            call interp_horz (nfdim,1,h_index,h_weights,   &
                              f3d(1,nlevs1,nf,jtime1),     &
                              f3d_interp(nlevs1,1,nf))
            call interp_horz (nfdim,1,h_index,h_weights,   &
                              f3d(1,nlevs1,nf,jtime1+1),   &
                              f3d_interp(nlevs1,2,nf))
          enddo  
!
! Temporally interpolate the horizontally interpolated fields at sfc.
          do nf=1,2            
            call interp_time (1,obs_time,time_files(jtime1:jtime1+1), &
                              f3d_interp(nlevs1:nlevs1,:,nf))
          enddo          
          obs_values(1,1:2)=f3d_interp(nlevs1,3,1:2)
!
! Set observation height, if not missing in original report, as surface 
! elevation in nature run. 
!
          if (l_comp_z .and. (nfz > 0) ) then
            obs_zlevs(1)=f2d_interp(1,3,nfz)
          endif
!
          if (ltest) then
            print *,' f3d_interp for surface observations:'
            print ('(6e14.5)'),f3d_interp(nlevs1,:,:)
          endif  ! test on ltest
!
! End of computations for surface observations 
!
        elseif (obs_nlevs == 1) then
!
! Interpolation for observations at a single level not indicated as 
! at the surface.  

! This includes some wind observations near the surface:
! Specifically, wind observations from ships and SSMI determined surface 
! winds are to be treated as observations at 20m.  For those particular 
! observations, the pressure p corresponding to the height of the observation
! is determined hydrostatically from a corresponding given z=20m. For all 
! other observation types, the p provided in the report as read is used.
! For the surface wind observations, the p written for the observation 
! is whatever was in the original file as is the station elevation. The 
! observation height, if not missing in the original report, is written 
! as the height determined from the nature run.
!
          if ((nint(obs_type) == 280) .or. (nint(obs_type) == 282) &
                .or. (nint(obs_type) == 283) ) then
            call compute_p_20m (ltest,jtime1,h_index,obs_time,  &
                                f2d_interp(1,3,nfp),h_weights,obs_plev)
          else
             obs_plev=obs_plevs(1)
          endif
!
! Check that pressure of observation is in range of NR
          lskip2=.false. 
          if (obs_plev < ak(1,2)+bk(1,2)*f2d_interp(1,3,nfp) ) then
            ierrors(5)=ierrors(5)+obs_nfields
            lskip2=.true.
          endif
          if (obs_plev > ak(nlevs1,2)+bk(nlevs1,2)*f2d_interp(1,3,nfp) ) then
            ierrors(6)=ierrors(6)+obs_nfields
            lskip2=.true.
          endif
!
          if (lskip2) then              
            obs_values(1,1:obs_nfields)=special_value 
            k_above_z(1)=nlevs1+1 ! set to > nlevs1 as special value
          else 
!
!   Determine which 2 hybrid levels need to be horizontally interpolated   
            call index_vert (k_iter,nlevs1,obs_plev,f2d_interp(1,3,nfp), &
                             ak(:,2),bk(:,2),k_above,k_below,p_above,p_below)
            k_above_z(1)=k_above  ! save index for computing z of obs
            if (ltest) then
              print *,'k_above,k_below,p_above,p_below=', &
                       k_above,k_below,p_above,p_below
            endif
!
!   Horizontally interpolate 2 levels
            f3d_interp(:,:,:)=special_value
            do nf=1,2
              call interp_horz (nfdim,2,h_index,h_weights,   &
                                f3d(1,k_above,nf,jtime1),    &
                                f3d_interp(k_above,1,nf)) 
              call interp_horz (nfdim,2,h_index,h_weights,   &
                                f3d(1,k_above,nf,jtime1+1),  &
                                f3d_interp(k_above,2,nf)) 
            enddo  
!
! Temporally and vertically fields but just at 2 required levels 
            do nf=1,2            
              call interp_time (2,obs_time,time_files(jtime1:jtime1+1), &
                                f3d_interp(k_above:k_below,:,nf))
              call interp_vert (obs_plev,p_above,p_below,          &
                                f3d_interp(k_above:k_below,3,nf),  &
                                obs_values(1,nf))  
            enddo          
!
! Compute height corresponding to observation at one level
            if (l_comp_z) then
              call compute_z_obs (ltest,1,jtime1,k_above_z,h_index,    &
                     f2d_interp(1,3,nfp),f2d_interp(1,3,nfz),obs_time, &
                     special_value,h_weights,obs_plevs,obs_zlevs,'sing')
            endif
!
            if (ltest) then
              print *,' f3d_interp for observation at a singel level:'
              print *,' k,f3d_interp(k,1:3times,1:fields3d)'
              do k=1,nlevs1
                if (f3d_interp(k,1,1) /= special_value) then 
                  print ('(i3,6e14.5)'),k,f3d_interp(k,:,:)
                endif
              enddo 
            endif  ! test on ltest
!              
          endif    ! test on lksip2
!
! End of computations for case of observations at a single level
!
        else 
!
! Begin case for interpolation of observations at mutiple levels
!   First horizontally interpolate full column of 3d fields at 2 times
          do nf=1,2
            call interp_horz (nfdim,nlevs1,h_index,h_weights,         &
                              f3d(1,1,nf,jtime1  ),f3d_interp(1,1,nf)) 
            call interp_horz (nfdim,nlevs1,h_index,h_weights,         &
                              f3d(1,1,nf,jtime1+1),f3d_interp(1,2,nf)) 
          enddo
!
!   Second, temporally interpolate columns
          do nf=1,2
            call interp_time (nlevs1,obs_time,time_files(jtime1:jtime1+1), &
                              f3d_interp(:,:,nf))
          enddo          
!
          if (ltest) then
            print *,' f3d_interp for observation at multiple levels:'
            print *,' k,f3d_interp(k,1:3times,1:fields3d)'
            do k=1,nlevs1
              print ('(i3,6e14.5)'),k,f3d_interp(k,:,:)
            enddo 
          endif  ! test on ltest
!
!   Third, perform vertical interpolations for each observation level
          do np=1,obs_nlevs
            obs_plev=obs_plevs(np)
!   Check that pressure of each observation is in range of NR
            lskip2=.false. 
            if (obs_plev < ak(1,2)+bk(1,2)*f2d_interp(1,3,nfp) ) then
              ierrors(5)=ierrors(5)+obs_nfields
              lskip2=.true.
            endif
            if (obs_plev > ak(nlevs1,2)+bk(nlevs1,2)*f2d_interp(1,3,nfp) ) then
              ierrors(6)=ierrors(6)+obs_nfields
              lskip2=.true.
            endif
!
            if (lskip2) then              
              obs_values(np,1:obs_nfields)=special_value 
              k_above_z(np)=nlevs1+1 ! set to > nlevs1 as special value 
            else         
              call index_vert (k_iter,nlevs1,obs_plev,f2d_interp(1,3,nfp), &
                               ak(:,2),bk(:,2),k_above,k_below,p_above,p_below)
              k_above_z(np)=k_above  ! save index for computing z of obs
              if (ltest) then
                print *,'np,k_above,k_below,p_above,p_below,obs_plev=', &
                         np,k_above,k_below,p_above,p_below,obs_plev
              endif
              do nf=1,2
                call interp_vert (obs_plev,p_above,p_below,          &
                                  f3d_interp(k_above:k_below,3,nf),  &
                                  obs_values(np,nf))  
              enddo
!
            endif  ! test on lskip2
!
          enddo    ! loop over observational p-levels
!
! Compute height corresponding to observations at all levels
          if (l_comp_z) then
            call compute_z_obs (ltest,obs_nlevs,jtime1,k_above_z,h_index, &
                       f2d_interp(1,3,nfp),f2d_interp(1,3,nfz),obs_time,  &
                       special_value,h_weights,obs_plevs,obs_zlevs,'mult')
          endif
!
! End of computations for case of observations at a multiple levels
!
        endif   ! check on whether surface, single level, or multiple level obs
!
! End required interpolations of NR fields
!
! Copy computed hydrostatic heights of observations 
        if (l_comp_z .and. (obs_nlevs>0) ) then
          obs_levels(1:obs_nlevs,2)=obs_zlevs(1:obs_nlevs)
        endif
!
! If q observation, and interpolation done in terms of RH, 
! then transform from RH to q
          if (lqfield .and. lq2rh .and. (obs_nlevs > 0)) then
            do np=1,obs_nlevs
              if (obs_values(np,1) /= special_value) then
                call rh2q (obs_values(np,nft),obs_plevs(np), &
                           obs_values(np,nfq),ierrors(12))
              endif
            enddo 
          endif
!
! If observations not processed since error detected, set to special values
        else     ! case for lskip1=.true. follows
          if (obs_nlevs > 0) then
            obs_values(1:obs_nlevs,1:obs_nfields)=special_value
          else
            obs_values(1,1:obs_nfields)=special_value
          endif  ! test on obs_nlevs>0 
        endif    ! test on lskip1
!
        if (ltest) then
          if (obs_nlevs > 0) then
            print *,' np, plev(np), obs_value(np,:):'
            do np=1,obs_nlevs
              print ('(i3,4f14.5)'),np,obs_plevs(np), &
                    obs_values(np,1:obs_nfields)
            enddo
          else
            print *,' obs skipped because obs_nlevs=',obs_nlevs 
          endif  ! test on obs_nlevs>0
        endif    ! test on ltest
!
! Count number of problems found
      if (obs_nlevs > 0) then       
        do np=1,obs_nlevs
          do nf=1,obs_nfields
            if (obs_values(np,nf) == special_value) then
              ierrors(7)=ierrors(7)+1 
            endif
          enddo
        enddo
      endif
!
! copy from rkind1(=4) to rkind3(=8) variables
      obs_info8=obs_info
!
      do nf=1,obs_ndim3
        do np=1,obs_nlevs
          if (obs_levels(np,nf) == special_value) then
            obs_levels8(np,nf)=bmiss
          else 
            obs_levels8(np,nf)=obs_levels(np,nf)
          endif
        enddo
      enddo
!
      do nf=1,obs_nfields
        do np=1,obs_nlevs
          if (obs_values(np,nf) == special_value) then
            obs_values8(np,nf)=bmiss
          else
            obs_values8(np,nf)=obs_values(np,nf)
          endif
        enddo
      enddo
!
      deallocate (f2d_interp,f3d_interp)
!
      end subroutine interp_NR_fields
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine interp_fields_for_rad (prof_ndim1,prof_ndim2,prof_ndim3, &
                             prof_ndim4,prof_nspot,nobs,n_errors,ierrors, &
                             prof_info,prof_2d,prof_3d,prof_plevs,        &
                             lprofile,lerror,ltest)
!
      implicit none
!
      integer :: prof_ndim1  ! number of model levels 
      integer :: prof_ndim2  ! number of 3-d fields
      integer :: prof_ndim3  ! number of 2-d fields
      integer :: prof_ndim4  ! number of scalar obs info
      integer :: prof_nspot  ! number of scalar obs info in spot header
      integer :: nobs        ! observation counter
      integer :: n_errors    ! size of array ierrors
      integer :: ierrors(n_errors)  ! counter for errors
      real(rkind1) :: prof_info(prof_ndim4)  
      real(rkind1) :: prof_2d(prof_ndim3)
      real(rkind1) :: prof_3d(prof_ndim1,prof_ndim2)
      real(rkind1) :: prof_plevs(prof_ndim1)
      logical :: lprofile   ! input: true if 3-d fields to be interpolated
                            ! otherwise, only 2-d fields interpolated
      logical :: lerror     ! output: true if error found
      logical :: ltest      ! input: true if test info to be printed
!
      logical :: lerror_ndim
      logical :: lskip1
      logical :: lskip2
!
      integer :: h_index(2,2)
      integer :: jtime1
      integer :: k
      integer :: nf
      integer :: np
      integer :: nt
      integer :: nt12 
      real(rkind1) :: h_weights(2,3)
      real(rkind1) :: obs_lat               ! latitude of observation 
      real(rkind1) :: obs_lon               ! longitude of observation 
      real(rkind1) :: obs_height            ! height of surface report
      real(rkind1) :: obs_time              ! time of observation 
      real(rkind1) :: p                     ! pressure
!
      real(rkind1), allocatable :: f2d_interp(:,:,:) ! horiz interp. 2d fields
      real(rkind1), allocatable :: f3d_interp(:,:,:) ! horiz interp. 3d fields
!
!
     
!          jlatN  --------------  jlonNW -------------- jlonNE              

!          jlatS  --------------      jlonSW ----- jlonSE
!
! Check that some input dimensions are correct
     lerror_ndim=.false.
     if (prof_ndim1 /= nlevs1-1)  lerror_ndim=.true.
     if (prof_ndim2 /= nfields3d) lerror_ndim=.true.
     if (prof_ndim3 /= nfields2d) lerror_ndim=.true.
     if (prof_ndim4 < 5         ) lerror_ndim=.true.
     if (lerror_ndim) then
       print *,' '
       print *,'Wrong dimensions input into subroutine ', &
               'interp_fields_for_rad'  
       print *,'  prof_ndim1,prof_ndim2,prof_ndim3,prof_ndim4 = ', &
                  prof_ndim1,prof_ndim2,prof_ndim3,prof_ndim4
       print *,'  nlevs,nfields2d,nfields3d = ',  &
                  nlevs1-1,nfields2d,nfields3d
       stop
     endif 
!
! Allocate some arrays
      allocate (f2d_interp(1,3,nfields2d))
      allocate (f3d_interp(nlevs1,3,nfields3d))
!
! Set some observation information
      obs_lat=   prof_info(2)
      obs_lon=   prof_info(3)
      obs_time=  prof_info(prof_ndim4)
!
! Check appropriateness of data time and find time index 
      lskip1=.false.
      call index_time (ntimes,obs_time,time_files(:),jtime1)
      if (ltest) then
        print *,' '
        print *,' '
        print *,' TEST OBSERVATION NUMBER ',nobs
        print *,'obs_time,time_files(:),jtime1=', &
                 obs_time,time_files(:),jtime1
      endif
      if (jtime1 == -1) then
        ierrors(1)=ierrors(1)+1
        lskip1=.true.
      endif
      if (jtime1 == -2) then
        ierrors(2)=ierrors(2)+1
        lskip1=.true.
      endif
!
! If acceptable location then continue with interpolation
      if (.not.lskip1) then  
        lerror=.false. 
!
! Compute grid box location that surrounds observation location
        call index_horz (obs_lat,obs_lon,nlats2,nlons2,nlonsP, &
                         glats2,h_index,h_weights,ltest)
!  
! Compute horizontally interpolated values of 2d fields at 2 times
        nt12=0
        do nt=jtime1,jtime1+1
          nt12=nt12+1
          do nf=1,nfields2d
            call interp_horz (nfdim,1,h_index,h_weights,    &
                      f2d(1,1,nf,nt),f2d_interp(1,nt12,nf))        
          enddo 
        enddo 
! 
! Compute temporally interpolated values of 2d fields 
        do nf=1,nfields2d
          call interp_time (1,obs_time,time_files(jtime1:jtime1+1), &
                            f2d_interp(:,:,nf))          
        enddo
!
! Fill output arrays for surface fields
        do nf=1,nfields2d 
          prof_2d(nf)=f2d_interp(1,3,nf)
        enddo
!
! Fill array of pressure levels
        do k=1,nlevs1-1
          prof_plevs(k)=ak(k,2)+bk(k,2)*f2d_interp(1,3,nfp)
        enddo
!
!       
        if (lprofile) then ! profile information required
!
! Begin case for interpolation of observations at mutiple levels
!   First horizontally interpolate full column of 3d fields at 2 times
          do nf=1,nfields3d
            call interp_horz (nfdim,nlevs1,h_index,h_weights,         &
                              f3d(1,1,nf,jtime1  ),f3d_interp(1,1,nf)) 
            call interp_horz (nfdim,nlevs1,h_index,h_weights,         &
                              f3d(1,1,nf,jtime1+1),f3d_interp(1,2,nf)) 
          enddo
!
!   Second, temporally interpolate columns
          do nf=1,nfields3d            
            call interp_time (nlevs1,obs_time,time_files(jtime1:jtime1+1), &
                              f3d_interp(:,:,nf))
          enddo          
!
          if (ltest) then
            print *,' f3d_interp for observation at multiple levels:'
            print *,' k,f3d_interp(k,1:3times,1:fields3d)'
            do k=1,nlevs1
              print ('(i3,6e14.5)'),k,f3d_interp(k,:,:)
            enddo 
          endif  ! test on ltest
!
! transform from rh back to q if interpolation was done using rh
          if (lqfield .and. lq2rh) then
            do np=kq2rh,nlevs1-1
                p=ak(k,2)+bk(k,2)*f2d_interp(1,3,nfp)
                call rh2q (f3d_interp(k,3,nft),p, &
                           f3d_interp(k,3,nfq),ierrors(12))
            enddo 
          endif
!
! Fill output arrays for vertical profiles 
          do nf=1,nfields3d 
            do k=1,nlevs1-1
              prof_3d(k,nf)=f3d_interp(k,3,nf)
            enddo
          enddo
!
        endif   !  test on lprofile
!
! End required interpolations of NR fields
!
!
! If observations not processed since error detected, set to special value
      else     ! case for lskip1=.false. follows
        lerror=.true.
        ierrors(5)=ierrors(5)+1 
      endif    ! test on lskip1
!
      if (ltest) then
        if (lerror) then 
          print *,' Error found in input info for profile, ', &
                  'so interpolations were skipped'
        else
          print *,'2-d interpolated fields follow:'
          do nf=1,nfields2d
            print ('(i3,1p1e15.5)'),nf,prof_2d(nf)
          enddo
          print *,'3-d interpolated fields follow:'
          do k=1,nlevs1-1
            print ('(i3,1p5e15.5)'),k,prof_3d(k,:)
          enddo
        endif    ! test on lerror
      endif      ! test on ltest
!
      deallocate (f2d_interp,f3d_interp)

      end subroutine interp_fields_for_rad
!
!
      end module m_interp_NR
!
! 
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine hydros (nlevs,ps,ak,bk,T,q,z)
!
! Compute hydrostatic height  (also see routine compute_p_qkscat) 
!
      use m_kinds
      use m_parameters
      implicit none
!      
      integer     :: nlevs
      real(rkind1) :: ps
      real(rkind1) :: ak(nlevs)  ! interface values
      real(rkind1) :: bk(nlevs)
      real(rkind1) :: T(nlevs)
      real(rkind1) :: q(nlevs)
      real(rkind1) :: z(nlevs)
!
      integer :: k
      real(rkind1) :: rfac
      real(rkind1) :: pa, pb
!
      rfac=Rdrygas/gravity
      pb=ak(nlevs)+bk(nlevs)*ps   ! value of p at interface below
      do k=nlevs-1,1,-1
        pa=ak(k)+bk(k)*ps         ! value of p at interface above
        z(k)=z(k+1)+rfac*T(k)*(one_k1+ratio4R*q(k))*log(pb/pa)
        pb=pa
      enddo
!
      end subroutine hydros 
!
! 
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine index_horz (obs_lat,obs_lon,nlats,nlons,nlons1, &
                             glats,h_index,h_weights,ltest)
!
! Determine parameters describing the horizontal grid point locations 
! surrounding an observation location. Then determine j-indexes and 
! weights for horizontal interpolation.
!
!          jlatN  --------------  jlonNW -------------- jlonNE              
!
!          jlatS  --------------      jlonSW ----- jlonSE
!
      use m_kinds
      implicit none
!
! input
      logical :: ltest
      integer :: nlats             ! nlats2 in calling program
      integer :: nlons(nlats)      ! nlons2 in calling program 
      integer :: nlons1(nlats)     ! nlonsP in calling program
      real(rkind1) :: glats(nlats) ! glats2 in calling program
      real(rkind1) :: obs_lat
      real(rkind1) :: obs_lon 
! output
      integer :: h_index(2,2)
      real(rkind1) :: h_weights(2,3)
! local
      integer :: jlatS 
      integer :: jlatN
      integer :: jlonSE
      integer :: jlonNE
      integer :: jlonSW
      integer :: jlonNW
      integer :: jlat1, jlat2
      integer :: jlat
      integer :: jlon
      integer :: j    
      real(rkind1) :: lonx
      real(rkind1) :: dlon
      real(rkind1) :: xlatS    
      real(rkind1) :: xlatN    
      real(rkind1) :: xlonSE    
      real(rkind1) :: xlonNE    
      real(rkind1) :: xlonSW    
      real(rkind1) :: xlonNW    
!
! Find <index for latitude directly south of requested lat 
      jlat1=(nlats*(obs_lat+90.)/180.)-4 
      jlat1=max(jlat1,1)
      jlat2=min(jlat1+9,nlats)
      jlat=0
      do j=jlat1,jlat2
        if (obs_lat > glats(j) ) jlat=j
      enddo
      jlatS=max(jlat,1)
      jlatS=min(jlatS,nlats-1)
      jlatN=jlatS+1 
      xlatS=glats(jlatS) 
      xlatN=glats(jlatN) 
!
! Find indexes for longitudes directly west of requested longitude
! Indexes may vary for latitudes if the grid is a reduced one
! Note that if the longitude diectly west is the last longitude, then 
! the index is 
      lonx=obs_lon
      if (obs_lon < 0. ) lonx=360.+obs_lon 
!
      dlon=360./nlons(jlatS)
      jlon=int(1+lonx/dlon)
      if (jlon == nlons(jlatS) ) then 
        jlonSE=1
        jlonSW=jlon
      else
        jlonSE=jlon+1
        jlonSW=jlon
      endif
      xlonSE=dlon*(jlonSE-1)
      xlonSW=dlon*(jlonSW-1)
!
      dlon=360./nlons(jlatN)
      jlon=int(1+lonx/dlon)
      if (jlon == nlons(jlatN) ) then
        jlonNE=1
        jlonNW=jlon
      else
        jlonNE=jlon+1
        jlonNW=jlon
      endif
      xlonNE=dlon*(jlonNE-1)
      xlonNW=dlon*(jlonNW-1)
!
! Based on location, compute weights and indexes
! Weights for 2 surrounding points to the north
      h_weights(1,1)=(lonx-xlonNW)/(xlonNE-xlonNW)
      h_weights(2,1)=one_k1-h_weights(1,1)
      h_index(1,1)=nlons1(jlatN)+jlonNE
      h_index(2,1)=nlons1(jlatN)+jlonNW
!
! Weights for 2 surrounding points to the south
      h_weights(1,2)=(lonx-xlonSW)/(xlonSE-xlonSW)
      h_weights(2,2)=one_k1-h_weights(1,2)
      h_index(1,2)=nlons1(jlatS)+jlonSE
      h_index(2,2)=nlons1(jlatS)+jlonSW
!
! Weights for 2 latitudes
      h_weights(1,3)=(obs_lat-xlatS)/(xlatN-xlatS)
      h_weights(2,3)=one_k1-h_weights(1,3)
!
      if (ltest) then
        print *,'jlatS,jlatN=',jlatS,jlatN
        print *,'jlonSE,jlonNE,jlonSW,jlonNW=',jlonSE,jlonNE,jlonSW,jlonNW
        print *,'xlatS,xlatN=',xlatS,xlatN
        print *,'xlonSE,xlonNE,xlonSW,xlonNW=',xlonSE,xlonNE,xlonSW,xlonNW
        print *,'h_index=',h_index
        print *,'h_weights=',h_weights
      endif 
!
      end subroutine index_horz
!
! 
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine index_time (ntimes,obs_time,time_files,jtime1)
!
! Determine the time of the file that preceeds the observation time
!
      use m_kinds
      implicit none
!
! input
      integer :: ntimes
      real(rkind1) :: obs_time
      real(rkind1) :: time_files(ntimes)
! output
      integer :: jtime1       
! local
      integer :: n
!  
      if (obs_time < time_files(1)) then
        jtime1=-2
      else
        do n=1,ntimes-1
          if (obs_time >= time_files(n)) then
            jtime1=n
          endif
        enddo
      endif
!   
      if (obs_time > time_files(ntimes)) then
         jtime1=-1
      endif
!      
      end subroutine index_time
!
! 
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine index_vert (k_iter,nlevs1,plev,ps,ak,bk,      &
                             k_above,k_below,p_above,p_below)
!
! Determine parameters describing the vertical grid point locations 
! surrounding an observation location.
! 
      use m_kinds
      implicit none
!
! input
      integer :: nlevs1
      integer :: k_iter  ! 
      real(rkind1) :: plev
      real(rkind1) :: ps
      real(rkind1) :: ak(nlevs1)
      real(rkind1) :: bk(nlevs1)
! output
      integer :: k_above
      integer :: k_below
      real(rkind1) :: p
      real(rkind1) :: p_above
      real(rkind1) :: p_below
! local 
      integer :: k 
      integer :: k_i
!
! Find the pair of grid points that bound the observation pressure by 
! successively dividing the range it is in by factors of 2.
      k_above=1
      k_below=nlevs1
      do k_i=1,k_iter
        k=(k_above+k_below)/2 
        p=ak(k)+bk(k)*ps
        if (plev >= p) then
          k_above=k
        else
          k_below=k
        endif
      enddo 
!
      if (k_above == nlevs1) then
         k_above=nlevs1-1
      endif
      k_below=k_above+1
      p_above=ak(k_above)+bk(k_above)*ps 
      p_below=ak(k_below)+bk(k_below)*ps
!      
      end subroutine index_vert
!
! 
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine interp_horz (nfdim,nlevs1,h_index,h_weights,f,finterp)
!
!  Perform a horizontal interpolation at mutiple p-levels based on 
!  some pre-determined parameters describing the bounding grid points
!  The algorithm is bilinear in terms of longitude and latitude.
!
      use m_kinds
      implicit none
!
! input
      integer :: nfdim
      integer :: nlevs1
      integer :: h_index(2,2)
      real(rkind1) :: h_weights(2,3)
      real(rkind2) :: f(nfdim,nlevs1)
! output
      real(rkind1) :: finterp(nlevs1) 
! local
      integer :: k    
      real(rkind1) :: fS, fN
!
      do k=1,nlevs1
        fN=h_weights(1,1)*f(h_index(1,1),k)+h_weights(2,1)*f(h_index(2,1),k)
        fS=h_weights(1,2)*f(h_index(1,2),k)+h_weights(2,2)*f(h_index(2,2),k)
        finterp(k)=h_weights(1,3)*fN+h_weights(2,3)*fS
      enddo
!        
      end subroutine interp_horz
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine interp_time (nlevs1,obs_time,times,f)
!
!  Perform linear interpolation in time
!
      use m_kinds
      implicit none
!
! input
      integer      :: nlevs1
      real(rkind1) :: obs_time
      real(rkind1) :: times(2)
      real(rkind1) :: f(nlevs1,3)
! output
      real(rkind1) :: fint
! local
      real(rkind1) :: r0, r1, r2
!
      r0=times(2)-times(1)
      r1=(obs_time-times(1))/r0
      r2=1.-r1
      f(:,3)=r2*f(:,1)+r1*f(:,2) 
!      
      end subroutine interp_time
!
! 
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine interp_vert (plev,p_above,p_below,f,f_Vinterp)
!
! Perform interpolation in the vertical that is linear in log p
! (This roughly corresponds to linear in height)
!
      use m_kinds
      implicit none
!
! input
      real(rkind1) :: plev
      real(rkind1) :: p_above
      real(rkind1) :: p_below
      real(rkind1) :: f(2)
! output
      real(rkind1) :: f_Vinterp
! local 
      real(rkind1) :: r0, r1, r2 
!
      r0=log(p_below/p_above)      
      r1=log(plev/p_above)/r0
      r2=1.-r1
      f_Vinterp=r2*f(1)+r1*f(2)
!      
      end subroutine interp_vert
!      
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      SUBROUTINE pole_tq (nfdim,nlevs1,nlats2,nlons2,nlonsP,f)
!
! Specify pole value of scalar fields (e.g., t, q) as mean along 
! first non-pole grid latitude.
!
      use m_kinds
      implicit none
!
! input      
      integer :: nfdim
      integer :: nlevs1
      integer :: nlats2
      integer :: nlons2(nlats2)
      integer :: nlonsP(nlats2)
! both input and output
      real(rkind2) :: f(nfdim,nlevs1)     
! local
      integer :: k, n
      integer :: np
      real(rkind1) :: fsum, fav
!
! Loop over levels 
      do k=1,nlevs1
!
! South pole   
        fsum=0.
        do n=1,nlons2(2)
          np=nlonsP(2)+n     
          fsum=fsum+f(np,k)
        enddo
        fav=fsum/nlons2(2)
        do n=1,nlons2(1)
          np=nlonsP(1)+n     
          f(np,k)=fav
        enddo
!
! North pole   
        fsum=0.
        do n=1,nlons2(nlats2-1)
          np=nlonsP(nlats2-1)+n     
          fsum=fsum+f(np,k)
        enddo
        fav=fsum/nlons2(nlats2-1)
        do n=1,nlons2(nlats2)
          np=nlonsP(nlats2)+n     
          f(np,k)=fav
        enddo
!          
      enddo  ! loop over levels
!
      end subroutine pole_tq
!      
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      SUBROUTINE pole_uv (nfdim,nlevs1,nlats2,nlons2,nlonsP,nlons_pole, &
                          pcoslon,psinlon,f)
!
! Specify wind field at the pole as equal to the zonal wave-number one 
! component of the field at the first non-pole grid latitude. u and v 
! components must be same amplitude with a phase shift of pi
!
      use m_kinds
      implicit none
!
! input      
      integer :: nfdim
      integer :: nlevs1
      integer :: nlats2
      integer :: nlons2(nlats2)
      integer :: nlonsP(nlats2)
      integer :: nlons_pole
      real(rkind1) :: pcoslon(nlons_pole)
      real(rkind1) :: psinlon(nlons_pole)
! both input and output
      real(rkind2) :: f(nfdim,nlevs1,2)     
! local
      integer :: k, n
      integer :: nf
      integer :: jS, jN
!
! The spectral coefficients are organized levels, 
! real or imaginary, S or N pole, u or v
      real(rkind1) :: coef(nlevs1,2,2,2) ! spectral coefficients
      real(rkind1) :: coefS1, coefS2, coefN1, coefN2
      real(rkind1) :: fac
!
! compute Fourier coef for zonal wavenumber 1      
      coef(:,:,:,:)=0.
      do n=1,nlons2(2)
        jS=nlonsP(2)+n
        jN=nlonsP(nlats2-1)+n
        do nf=1,2               ! 1=u, 2=v
          do k=1,nlevs1
            coef(k,1,1,nf)=coef(k,1,1,nf)+pcoslon(n)*f(jS,k,nf)
            coef(k,2,1,nf)=coef(k,2,1,nf)-psinlon(n)*f(jS,k,nf)
            coef(k,1,2,nf)=coef(k,1,2,nf)+pcoslon(n)*f(jN,k,nf)
            coef(k,2,2,nf)=coef(k,2,2,nf)-psinlon(n)*f(jN,k,nf)
          enddo
        enddo
      enddo
      fac=2./nlons2(2)  ! factor of 2 because end result is 2*Real Part  
      coef(:,:,:,:)=fac*coef(:,:,:,:)
!
! Average the coefficients of u estimated from the u field directly 
! and the estimate from the v field using coef(u) = i*coef(v) for S.H.
! and reversed sign for N.H.
      do k=1,nlevs1
        coefS1=0.5*(coef(k,1,1,1)-coef(k,2,1,2))
        coefS2=0.5*(coef(k,2,1,1)+coef(k,1,1,2))      
        coefN1=0.5*(coef(k,1,2,1)+coef(k,2,2,2))
        coefN2=0.5*(coef(k,2,2,1)-coef(k,1,2,2))
!
! Copy the averaged estimates for the u coefficients
        coef(k,1,1,1)=coefS1
        coef(k,2,1,1)=coefS2
        coef(k,1,2,1)=coefN1
        coef(k,2,2,1)=coefN2
!
! Specify the coefficients for v using coef(v) = -i*coef(u) for S.H.
! and reversed sign for N.H.
        coef(k,1,1,2)= coefS2
        coef(k,2,1,2)=-coefS1
        coef(k,1,2,2)=-coefN2
        coef(k,2,2,2)= coefN1
      enddo  ! loop over k
!
! Construct wind at both poles
      do n=1,nlons2(1)
        jS=nlonsP(1)+n
        jN=nlonsP(nlats2)+n
        do nf=1,2               ! 1=u, 2=v
          do k=1,nlevs1
            f(jS,k,nf)=coef(k,1,1,nf)*pcoslon(n)-coef(k,2,1,nf)*psinlon(n)
            f(jN,k,nf)=coef(k,1,2,nf)*pcoslon(n)-coef(k,2,2,nf)*psinlon(n)
          enddo
        enddo
      enddo
!
      end subroutine pole_uv
!
! 
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      Subroutine q2rh (klow,nfdim,nlevs1,ak,bk,t,ps,q)
!
! Convert specific humidity to realtive humidity using table of 
! saturation vapor pressures for levels k=klow,nlevs1
!
      use m_kinds
      use m_relhum
      implicit none
!
! input      
      integer :: klow
      integer :: nfdim
      integer :: nlevs1
      real(rkind1) :: ak(nlevs1), bk(nlevs1)
      real(rkind2) :: ps(nfdim)
      real(rkind2) :: t(nfdim,nlevs1)
! both input and output
      real(rkind2) :: q(nfdim,nlevs1)     
! local
      integer :: errors
      integer :: k, j
      real(rkind2) :: p
!
      errors=0
      do k=klow,nlevs1
        do j=1,nfdim
          p=ak(k)+bk(k)*ps(j)
          call relhum_q2rh (t(j,k),p,q(j,k),errors)
        enddo
      enddo
!
      if (errors > 0) then
        print *,' Problem in q2rh: p < v found',errors,' times'
        stop
      endif
!
      end subroutine q2rh 
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      SUBROUTINE read_gridinfo (nfdim,nlats2,nlevs1,file_unit,        &
                                file_name,nlons2,nlonsP,glats2,ak,bk) 
!
! Read file with grid info describing NR 
!
      use m_kinds
      implicit none
!
! input
      integer :: nfdim
      integer :: nlats2
      integer :: nlevs1
      integer :: file_unit
      character(len=*) :: file_name
! 
! output
      integer :: nlons2(nlats2)
      integer :: nlonsP(nlats2)
      real(rkind1) :: glats2(nlats2)   
      real(rkind1) :: ak(nlevs1,2)  ! values for interface and data levels 
      real(rkind1) :: bk(nlevs1,2)  ! (Note: surface included in both sets) 
!
! local
      integer :: n 
      integer :: nhalf, nhalf2
      integer :: n10, n1, n2
      integer :: np
      integer :: nc
      character(len=1) :: blank
!
      nhalf2=nlats2/2
      nhalf=nhalf2-1
!
! open file
      open (file_unit, file=file_name, form='formatted')
      print *,' '
      print ('(3a,i3)'),' File=',file_name, &
             ' opened for reading grid info on unit=',file_unit 
!
! Read nlons
      read (file_unit,'(a1)') blank  ! skip record
      read (file_unit,'(a1)') blank  ! skip record
      do n=2,nhalf2
        read (file_unit,'(i6)') nlons2(n)
        nlons2(nlats2-n+1)=nlons2(n)
      enddo
!
! Read Gaussian latitudes N. to S.
      n10=nhalf/10
      if (n10*10 /= nhalf) then
        n10=n10+1
      endif
      read (file_unit,'(a1)') blank  ! skip record
      read (file_unit,'(a1)') blank  ! skip record
      n2=1
      do n=1,n10
        n1=n2+1
        n2=min(n1+9,nhalf2)
        read (file_unit,'(10f8.3)') glats2(n1:n2)
      enddo
!  Fix order
      do n=2,nhalf2
        glats2(nlats2-n+1)=glats2(n)
        glats2(n)=-glats2(n)
      enddo
!
! Read ak and bk for layer interfaces 
      read (file_unit,'(a1)') blank  ! skip record
      read (file_unit,'(a1)') blank  ! skip record
      do n=1,nlevs1
        read (file_unit,'(f24.12)') ak(n,1)
      enddo 
      do n=1,nlevs1
        read (file_unit,'(f24.12)') bk(n,1)
      enddo 
!  Determine ak and bk at data levels as simple average
!  Note that the surface is included in both as an interface and data level
      do n=1,nlevs1-1 
        ak(n,2)=0.5*(ak(n,1)+ak(n+1,1))
        bk(n,2)=0.5*(bk(n,1)+bk(n+1,1))
      enddo
      ak(nlevs1,2)=ak(nlevs1,1)  ! surface value
      bk(nlevs1,2)=bk(nlevs1,1)  ! surface value
!
! Always set pole values as shown here
      glats2(1)=-90.
      glats2(nlats2)=90.
      nlons2(1)=nlons2(2)
      nlons2(nlats2)=nlons2(nlats2-1)
!
! Set nlonsP 
      call table_nlonsP (nlats2,nlons2,nlonsP)
!
! Check nfdim
      nc=sum(nlons2)
      if (nc /= nfdim) then
        print *,' '
        print *,' sum(nlons2)=',nc,' /= ','nfdim=',nfdim
        print *,' reset nfdim to nc'
        stop
      endif
!
      close (file_unit)
!
      end subroutine read_gridinfo
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      Subroutine read_NR (nfdim,nlevs1,nlats2,ntimes,nfields2d,nfields3d,     &
                          nfields,file_data_unit,file_data_len,nlons2,nlonsP, &
                          f_names,file_data_2D,file_data_ps,file_data_3D,     &
                          f2d,f3d,ntime1)
!
!  Read Nature Run at consecutive times
!
      use m_kinds
      use m_interp_NR, only : find_name
      implicit none
!
! parameter
      logical, parameter :: lprint=.false. ! true if more info to be printed
!
! input
      integer :: nfdim     ! size of field data array for 1 level 
      integer :: nlevs1    ! number of data levels + surface
      integer :: nlats2    ! number of latitudes plus poles
      integer :: ntimes    ! number of times to read during this call 
      integer :: ntime1    ! index of first time to read
      integer :: nfields2d ! number of 2 d fields to be read
      integer :: nfields3d ! number of 3 d fields to be read
      integer :: nfields   ! sum of nfields2d+nfields3d
      integer :: file_data_unit  ! unitnumber for files to be read
      integer :: file_data_len   ! character length of file name
      integer :: nlons2(nlats2)  ! number of longitudes for each latitude    
      integer :: nlonsP(nlats2)  ! index value -1 for first point in lat.  
      character(len=*) :: f_names(nfields) ! names of fields to be read
      character(len=file_data_len) :: file_data_2D ! file name for 2d fields
      character(len=file_data_len) :: file_data_ps ! file name for ps/zs fields
      character(len=file_data_len) :: file_data_3D ! file name for 3d fields
!
! output 
      real(rkind2) :: f2d(nfdim,     1,nfields2d,ntimes)
      real(rkind2) :: f3d(nfdim,nlevs1,nfields3d,ntimes)
!
      integer :: k
      integer*4 :: kl
      integer :: nf
      integer :: nfp, nfz
      integer :: nfx, nfx3
      integer :: nfdim0
      integer, parameter :: nsurf=12  ! number of fields on surface data file
      integer :: nt, nt1
      real(rkind2), allocatable :: r4data(:)
      character(len=file_data_len) :: file_name
      character(len=4) :: varchar
      character(len=len(f_names)) :: f_name
      character(len=9) :: file_times
!
      print *,' '
      print *,' Begin read of NR data'
!
! set some arrays
      nfdim0=nfdim-nlons2(1)-nlons2(nlats2)
      allocate (r4data(1:nfdim0))
!
      if (ntimes > 9) then
        print *,' '
        print *,' ntimes=',ntimes,' > 9'
        print *,' so length and use of string file_times must be changed'
        stop       
      else
        file_times='123456789'
      endif
!
! read ps and (optionally) zsfc fields
      call find_name (nfields2d,f_names(1:nfields2d),'pres',nfp,.true.)
      call find_name (nfields2d,f_names(1:nfields2d),'zsfc',nfz,.false.)
      file_name=file_data_ps
      do nt=1,ntimes
        nt1=ntime1+nt-1
        file_name(file_data_len:file_data_len)=file_times(nt1:nt1)
        open (file_data_unit, file=file_name, form='unformatted')
        print ('(3a,i3)'),' File=',file_name, &
             ' opened for reading ps data on unit=',file_data_unit 
!
        do nf=1,2  ! there are 2 records on each ps file        
          read (file_data_unit) varchar,kl,r4data
          if (varchar == 'pres') then
            call reorder_NS (nfdim0,nfdim,nlats2,nlons2,nlonsP, &
                             r4data(:),f2d(:,1,nfp,nt))
            if (lprint) print ('(3a)'),' Field=',varchar, &
                     ' on surface pressure file is used'
          elseif (varchar == 'zsfc') then 
            if (nfz /= 0) then
              call reorder_NS (nfdim0,nfdim,nlats2,nlons2,nlonsP, &
                               r4data(:),f2d(:,1,nfz,nt))
              if (lprint) print ('(3a)'),' Field=',varchar, &
                      ' on surface pressure file is used'
            else
              if (lprint) print ('(3a)'),' Field=',varchar, &
                      ' on sfc p/z file is not used'
            endif
          else
            print *,' '
            print *,' Wrong file acceesed for surface pressure:'
            print *,' file=',file_name,' varchar=',varchar
            stop 
          endif
        enddo     ! loop over 2 records 
        close (file_data_unit)
        if (lprint) print ('(3a)'),' file=',file_name,' read'
      enddo       ! loop over ntimes
!
! read required 3-D fields      
      file_name=file_data_3d
      do nf=1,nfields3d
        f_name=f_names(nfields2d+nf)
        file_name(1:1)=f_name(1:1)
        if (file_name(1:1)=='s') file_name(1:1)='q'
        do nt=1,ntimes
          nt1=ntime1+nt-1
          file_name(file_data_len:file_data_len)=file_times(nt1:nt1)
          open (file_data_unit, file=file_name, form='unformatted')
          print ('(3a,i3)'),' File=',file_name, &
               ' opened for reading 3D data on unit=',file_data_unit 
!        
          do k=1,nlevs1-1
            read (file_data_unit) varchar,kl,r4data
            if (varchar == f_name) then
              call reorder_NS (nfdim0,nfdim,nlats2,nlons2,nlonsP, &
                               r4data(:),f3d(:,kl,nf,nt))
            else
              print *,' '
              print *,' Wrong file acceesed for 3-D data: file=',file_name, &
                      '  varchar=',varchar,' f_name=',f_name
              stop 
            endif
          enddo ! loop over levels
          close (file_data_unit)
          if (lprint) print ('(3a)'),' file=',file_name,' read'
        enddo ! loop over ntimes
      enddo ! loop over nfields3d
!
! read required 2-D and surface fields (except for press and zsfc)
!
      nfields=nfields2d+nfields3d
      file_name=file_data_2d
!
      do nt=1,ntimes
        nt1=ntime1+nt-1
        file_name(file_data_len:file_data_len)=file_times(nt1:nt1)
        open (file_data_unit, file=file_name, form='unformatted')
        print ('(3a,i3)'),' File=',file_name, &
           ' opened for reading surface data on unit=',file_data_unit 
!        
        do nf=1,nsurf
          read (file_data_unit) varchar,kl,r4data
          f_name=varchar
          call find_name (nfields2d,f_names(1:nfields2d),f_name,nfx,.false.)
          if (f_name == 'dt2m') f_name='sphu'
          if (f_name == 'u10m') f_name='uwnd'
          if (f_name == 'v10m') f_name='vwnd'
          if (f_name == 'tp2m') f_name='temp'
          call find_name (nfields3d,f_names(nfields2d+1:nfields), &
                          f_name,nfx3,.false.)
!
          if (nfx3 > 0) then
            call reorder_NS (nfdim0,nfdim,nlats2,nlons2,nlonsP,r4data(:), &
                             f3d(:,nlevs1,nfx3,nt))
            if (lprint) print ('(3a)'),' Field=',varchar, &
                    ' on surface file is used as sfc of 3-D field'
          elseif (nfx > 0) then
            call reorder_NS (nfdim0,nfdim,nlats2,nlons2,nlonsP,r4data(:), &
                           f2d(:,1,nfx,nt))
            if (lprint) print ('(3a)'),' Field=',varchar, &
                    ' on surface file is used as a 2-D field'
          else
            if (lprint) print ('(3a)'),' Field=',varchar, &
                    ' on surface file is not used'
          endif
        enddo ! loop over 2-D fields
        close (file_data_unit)
        if (lprint) print ('(3a)'),' file=',file_name,' read'
      enddo ! loop over ntimes
!
      deallocate (r4data)
      print ('(a,i2,a)'),' NR fields read for',ntimes,' times'
!
      end subroutine read_NR 
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine reorder_NS (nfdim0,nfdim,nlats2,nlons2,nlonsP,r,f)
!
!  Reorder latitudes so that final result is S pole first and  N. pole last
!
      use m_kinds
      implicit none
!
! input
      integer :: nfdim0
      integer :: nfdim
      integer :: nlats2
      integer :: nlons2(nlats2)      
      integer :: nlonsP(nlats2)   
      real(rkind2) :: r(nfdim0)
! output 
      real(rkind2) :: f(nfdim)
!
      integer :: n, n1r, n2r, n1f, n2f
!
      do n=2,nlats2-1
        n1r=nlonsP(n)-nlons2(1)+1
        n2r=n1r+nlons2(n)-1
        n1f=nlonsP(nlats2-n+1)+1
        n2f=n1f+nlons2(n)-1
        f(n1f:n2f)=r(n1r:n2r)      
      enddo
!
      end subroutine reorder_NS 
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      Subroutine read_NR_test (nfdim,nlevs1,nlats2,ntimes,nfields2d, &
              nfields3d,nfields,nlons2,nlonsP,f_names,f2d,f3d)
!
!  Fill "fake" NR fields for testing interpolation software.
!  Only fill non-pole latitudes. 
!
      use m_kinds
      implicit none
!
! input
      integer :: nfdim
      integer :: nlevs1
      integer :: nlats2
      integer :: ntimes
      integer :: nfields2d
      integer :: nfields3d
      integer :: nfields
      integer :: nlons2(nlats2)      
      integer :: nlonsP(nlats2)   
      character(len=*) :: f_names(nfields) 
!
! output 
      real(rkind2) :: f2d(nfdim,    1,nfields2d,ntimes)
      real(rkind2) :: f3d(nfdim,nlevs1,nfields3d,ntimes)
!
      integer :: nt, nf, np, k, n, i
      real(rkind1) :: ft, ff, fk, fy, fx
      real(rkind1) :: pt, py, px
!
      print *,' '
!
! fake 2D fields
      do nt=1,ntimes
        pt=1.e5+(nt-1)*1000.
        do nf=1,nfields2d
          do n=2,nlats2-1  
            py=pt+100.*n
            do i=1,nlons2(n)
              px=py+10*i
              np=nlonsP(n)+i
              f2d(np,1,nf,nt)=px
            enddo
          enddo
          print *,' NR field set for time=',nt,' and 2d field ',nf, &
                  ' with name=',f_names(nf:nf) 
        enddo  
      enddo  
!
! fake 3D fields
      do nt=1,ntimes
        ft=nt*20.
        do nf=1,nfields3d
          ff=ft+nf*10.
!        
          do k=1,nlevs1
            fk=ff+k
            do n=2,nlats2-1  
              fy=fk+0.5*n
              do i=1,nlons2(n)
                fx=fy+.2*i
                np=nlonsP(n)+i
                f3d(np,k,nf,nt)=fx
              enddo
            enddo
          enddo
!
          print *,' NR field set for time=',nt,' and 3d field ',nf, &
                  ' with name=',f_names(nfields2d+nf:nfields2d+nf) 
        enddo
      enddo  
!
      end subroutine read_NR_test     
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      SUBROUTINE rh2q (t,p,q,ierror)
!
!  Convert relative humidity to specific humidity
!  Note that real input and output here assumed to be rkind1
!
      use m_kinds
      use m_relhum    
      implicit none
!
      real(rkind1) :: p
      real(rkind1) :: t
      real(rkind1) :: q     
      integer      :: ierror
      integer      :: ier
!
      ier=0
      call relhum_rh2q (t,p,q,ier)
      if (ier > 0) then
        q=zero_k1
        ierror=ierror+1
      endif
!
      end subroutine rh2q 
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      SUBROUTINE set_gridinfo_test (nfdim,nlats2,nlevs1,nlons2,nlonsP, &
                                    glats2,ak,bk) 
!
!  Fill 'fake' NR grid info for testing interpolation software 
!
      use m_kinds
      implicit none
!  
! input
      integer :: nfdim
      integer :: nlats2
      integer :: nlevs1
! 
! output
      integer :: nlons2(nlats2)
      integer :: nlonsP(nlats2)
      real(rkind1) :: glats2(nlats2)   
      real(rkind1) :: ak(nlevs1,2) ! see routine read_gridinfo  
      real(rkind1) :: bk(nlevs1,2)  
!
! local
      integer :: n, np, nc
!
! some fake ak and bk for testing
      ak(:,1)=0.
      bk(:,1)=0.
      bk(nlevs1,1)=1.
      bk(5,1)=.80
      bk(4,1)=.60   
      bk(3,1)=.20
      ak(3,1)=100.
      ak(2,1)=100.
      ak(1,1)=10. 
      do n=1,nlevs1-1 
        ak(n,2)=0.5*(ak(n,1)+ak(n+1,1))
        bk(n,2)=0.5*(bk(n,1)+bk(n+1,1))
      enddo
      ak(nlevs1,2)=ak(nlevs1,1)
      bk(nlevs1,2)=bk(nlevs1,1)

!
! some fake latitudes for testing
      glats2(2)=-70.
      glats2(3)=-40.
      glats2(4)=-10.
      glats2(5)= 10.
      glats2(6)= 50.
      glats2(7)= 70.
!
! some fake values of nlons2
      nlons2(2)=4
      nlons2(3)=6
      nlons2(4)=8
      nlons2(5)=8
      nlons2(6)=6
      nlons2(7)=4
!
! always set pole values as shown here
      glats2(1)=-90.
      glats2(nlats2)=90.
      nlons2(1)=nlons2(2)
      nlons2(nlats2)=nlons2(nlats2-1)
!
      call table_nlonsP (nlats2,nlons2,nlonsP)
!
! check nfdim
      nc=sum(nlons2)
      if (nc /= nfdim) then
        print *,' '
        print *,' sum(nlons2)=',nc,' /= ','nfdim=',nfdim
        print *,' reset nfdim to nc'
        stop
      endif
!
      end subroutine set_gridinfo_test
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      SUBROUTINE table_nlonsP (nlats,nlons,nlonsP)
!
!  Compute nlonsP from nlons
!
      use m_kinds
      implicit none
!
! input      
      integer :: nlats
      integer :: nlons(nlats)
! output
      integer :: nlonsP(nlats)
! local
      integer :: n
! 
      nlonsP(1)=0
      do n=2,nlats
        nlonsP(n)=nlonsP(n-1)+nlons(n-1)
      enddo
!
      print *,' '
      print *,' Table for nlonsP filled'
!
      end subroutine table_nlonsP
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      SUBROUTINE table_sincos (nlons,pcoslon,psinlon)
!
!  Create table of sines and cosines of longitudes
!
      use m_kinds
      implicit none
!
! input      
      integer :: nlons

! output
      real(rkind1) :: pcoslon(nlons)
      real(rkind1) :: psinlon(nlons)
!local
      integer :: n 
      real(rkind1):: twopi, pifac, a

      twopi=8.*atan(1.)
      pifac=twopi/nlons
      do n=1,nlons
        a=pifac*(n-1)
        pcoslon(n)=cos(a)
        psinlon(n)=sin(a)
      enddo
!
      end subroutine table_sincos
!
! 
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      Subroutine td2q (nfdim,ps,q)
!
! Replace dew point temperature (input in array q) to specific humidity 
! (output array q) at the surface only.
!
      use m_kinds
      use m_relhum
      implicit none
!
! input      
      integer :: nfdim
      real(rkind2) :: ps(nfdim)
! both input and output
! on input, q is actually TD, on output is q
      real(rkind2) :: q(nfdim)     
! local
      integer :: j
!
      do j=1,nfdim
        call relhum_td2q (q(j),ps(j))
      enddo
!
      end subroutine td2q 
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      SUBROUTINE write_field_test (nfdim,nlevs1,nlats2,nlons2,nlonsP,nt, &
                                   f,f_name)
!
!  Print some fields for testing (only appropriate resolution is very low) 
!
      use m_kinds
      implicit none
!
! input      
      integer :: nfdim
      integer :: nlevs1
      integer :: nlats2
      integer :: nlons2(nlats2)
      integer :: nlonsP(nlats2)
      integer :: nt
      real(rkind2) :: f(nfdim,nlevs1)
      character(len=*) :: f_name
!
      integer :: k, n, i1, i2
!
      print *,' '
      print *,'field = ',f_name,'  nt=',nt,'  levs=',nlevs1
!
      do k=1,nlevs1
        print *,'level ',k,'  for field=',f_name, ' at time=',nt
        do n=nlats2,1,-1
          i1=nlonsP(n)+1
          i2=i1+nlons2(n)-1
          print ('(i2,12e13.5)'),n,f(i1:i2,k) 
        enddo
      enddo
!
      end subroutine write_field_test
!