driver.F90

 program interp_inc

!---------------------------------------------------------------------
!
! Read a gaussian atmospheric increment file in netcdf.  Interpolate
! all fields to another gaussian resolution.  Output the result
! in another netcdf file.
!
! Namelist variables:
! -------------------
! lon_out               - 'i' dimension of output gaussian grid
! lat_out               - 'j' dimension of output gaussian grid
! lev                   - Number of vertical levels.  Must be
!                         the same for the input and output grids.
! infile                - Path/name of input gaussian increment
!                         file (netcdf)
! outfile               - Path/name of output gaussian increment
!                         file (netcdf)
!
! 2019-10-24        Initial version.
! 2020-01-27        Martin - added in some simple MPI to speed up a bit
!
!---------------------------------------------------------------------

 use netcdf
 use mpi_f08
#ifdef IP_V4
 use ip_mod, only: ipolates, ipolatev
#endif

 implicit none

 integer, parameter :: num_recs = 12

! Declare externals
 external :: w3tagb, netcdf_err, splat, w3tage
#ifndef IP_V4
 external :: ipolates, ipolatev
#endif

 character(len=128) :: outfile, infile
 character(len=11)  :: records(num_recs) 

 integer :: i, j, mi, iret, mo, rec
 integer :: lon_in, lat_in
 integer :: lon_out, lat_out
 integer :: lev, ilev, lev_in
 integer :: ncid_in, id_var
 integer :: ncid_out, error
 integer :: dim_lon_out, dim_lat_out
 integer :: dim_lev_out, dim_ilev_out
 integer :: id_u_inc_out, id_v_inc_out
 integer :: id_lon_out, id_lat_out, id_lev_out
 integer :: id_pfull_out, id_ilev_out
 integer :: id_hyai_out, id_hybi_out
 integer :: id_delp_inc_out, id_delz_inc_out
 integer :: id_t_inc_out, id_sphum_inc_out
 integer :: id_liq_wat_inc_out, id_o3mr_inc_out
 integer :: id_icmr_inc_out, id_dim
 integer :: id_rwmr_inc_out, id_snmr_inc_out, id_grle_inc_out
 integer :: header_buffer_val = 16384
 integer :: kgds_in(200), kgds_out(200)
 integer :: ip, ipopt(20), no
 integer :: klev
 integer, allocatable :: ibi(:), ibo(:), levs(:)

 integer :: mpierr, mype, npes

 type (MPI_Status) :: mpistat

 logical*1, allocatable :: li(:,:), lo(:,:)

 real, allocatable :: dummy_in(:,:,:)
 real, allocatable :: dummy_out(:,:,:)

 real(8) :: rad2deg,dlondeg
 real(8), allocatable :: latitude_in(:), longitude_in(:)
 real(8), allocatable :: latitude_out(:), longitude_out(:)
 real(8), allocatable :: slat(:), wlat(:)
 real(8), allocatable :: rlon(:), rlat(:), crot(:), srot(:)
 real(8), allocatable :: gi(:,:), gi2(:,:), go(:,:), go2(:,:), go3(:,:)
 real(8), allocatable :: send_layer(:), recv_layer(:)

 logical :: readvar

 ! NOTE: u_inc,v_inc must be consecutive
 data records /'u_inc', 'v_inc', 'delp_inc', 'delz_inc', 'T_inc', &
               'sphum_inc', 'liq_wat_inc', 'o3mr_inc', 'icmr_inc', &
               'rwmr_inc', 'snmr_inc', 'grle_inc' /

 namelist /setup/ lon_out, lat_out, outfile, infile, lev


!-----------------------------------------------------------------
! MPI initialization
call mpi_init(mpierr)
call mpi_comm_rank(mpi_comm_world, mype, mpierr)
call mpi_comm_size(mpi_comm_world, npes, mpierr)
!-----------------------------------------------------------------

!-----------------------------------------------------------------
! Open and create output file records.  These will be filled
! with data below.
!-----------------------------------------------------------------

 if (mype == npes-1)  call w3tagb('INTERP_INC', 2019, 100, 0, 'EMC')

 if (mype == npes-1) print*,'- READ SETUP NAMELIST'
 open (43, file="./fort.43")
 read (43, nml=setup, iostat=error)
 if (error /= 0) then
   print*,"- FATAL ERROR READING NAMELIST. ISTAT IS: ", error
   stop 44
 endif
 close (43)

 if (mype == npes-1) print*,"- WILL INTERPOLATE TO GAUSSIAN GRID OF DIMENSION ",lon_out, lat_out

! Set constants
 rad2deg = 180.0_8 / (4.0_8 * atan(1.0_8))
 dlondeg = 360.0_8 / real(lon_out,8)

 ilev=lev+1

 call mpi_barrier(mpi_comm_world, mpierr)
 if (mype == npes-1) then
   print*,'- OPEN OUTPUT FILE: ', trim(outfile)
  
   error = nf90_create(outfile, cmode=IOR(NF90_CLOBBER,NF90_NETCDF4), ncid=ncid_out)
   call netcdf_err(error, 'CREATING FILE='//trim(outfile) )
  
   error = nf90_def_dim(ncid_out, 'lon', lon_out, dim_lon_out)
   call netcdf_err(error, 'defining dimension lon for file='//trim(outfile) )
  
   error = nf90_def_dim(ncid_out, 'lat', lat_out, dim_lat_out)
   call netcdf_err(error, 'defining dimension lat for file='//trim(outfile) )
  
   error = nf90_def_dim(ncid_out, 'lev', lev, dim_lev_out)
   call netcdf_err(error, 'defining dimension lev for file='//trim(outfile) )
  
   error = nf90_def_dim(ncid_out, 'ilev', ilev, dim_ilev_out)
   call netcdf_err(error, 'defining dimension ilev for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'lon', nf90_double, (/dim_lon_out/), id_lon_out)
   call netcdf_err(error, 'defining variable lon for file='//trim(outfile) )
  
   error = nf90_put_att(ncid_out, id_lon_out, "units", "degrees_east")
   call netcdf_err(error, 'define lon attribute for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'lat', nf90_double, (/dim_lat_out/), id_lat_out)
   call netcdf_err(error, 'defining varable lat for file='//trim(outfile) )
  
   error = nf90_put_att(ncid_out, id_lat_out, "units", "degrees_north")
   call netcdf_err(error, 'defining lat att for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'lev', nf90_float, (/dim_lev_out/), id_lev_out)
   call netcdf_err(error, 'defining variable lev for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'pfull', nf90_float, (/dim_lev_out/), id_pfull_out)
   call netcdf_err(error, 'defining variable pfull for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'ilev', nf90_float, (/dim_ilev_out/), id_ilev_out)
   call netcdf_err(error, 'defining variable ilev for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'hyai', nf90_float, (/dim_ilev_out/), id_hyai_out)
   call netcdf_err(error, 'defining variable hyai for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'hybi', nf90_float, (/dim_ilev_out/), id_hybi_out)
   call netcdf_err(error, 'defining variable hybi for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'u_inc', nf90_float, (/dim_lon_out,dim_lat_out,dim_lev_out/), id_u_inc_out)
   call netcdf_err(error, 'defining variable u_inc for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'v_inc', nf90_float, (/dim_lon_out,dim_lat_out,dim_lev_out/), id_v_inc_out)
   call netcdf_err(error, 'defining variable v_inc for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'delp_inc', nf90_float, (/dim_lon_out,dim_lat_out,dim_lev_out/), id_delp_inc_out)
   call netcdf_err(error, 'defining variable delp_inc for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'delz_inc', nf90_float, (/dim_lon_out,dim_lat_out,dim_lev_out/), id_delz_inc_out)
   call netcdf_err(error, 'defining variable delz_inc for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'T_inc', nf90_float, (/dim_lon_out,dim_lat_out,dim_lev_out/), id_t_inc_out)
   call netcdf_err(error, 'defining variable t_inc for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'sphum_inc', nf90_float, (/dim_lon_out,dim_lat_out,dim_lev_out/), id_sphum_inc_out)
   call netcdf_err(error, 'defining variable sphum_inc for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'liq_wat_inc', nf90_float, (/dim_lon_out,dim_lat_out,dim_lev_out/), id_liq_wat_inc_out)
   call netcdf_err(error, 'defining variable liq_wat_inc for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'o3mr_inc', nf90_float, (/dim_lon_out,dim_lat_out,dim_lev_out/), id_o3mr_inc_out)
   call netcdf_err(error, 'defining variable o3mr_inc for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'icmr_inc', nf90_float, (/dim_lon_out,dim_lat_out,dim_lev_out/), id_icmr_inc_out)
   call netcdf_err(error, 'defining variable icmr_inc for file='//trim(outfile) )

   error = nf90_def_var(ncid_out, 'rwmr_inc', nf90_float, (/dim_lon_out,dim_lat_out,dim_lev_out/), id_rwmr_inc_out)
   call netcdf_err(error, 'defining variable rwmr_inc for file='//trim(outfile) )
  
   error = nf90_def_var(ncid_out, 'snmr_inc', nf90_float, (/dim_lon_out,dim_lat_out,dim_lev_out/), id_snmr_inc_out)
   call netcdf_err(error, 'defining variable snmr_inc for file='//trim(outfile) )

   error = nf90_def_var(ncid_out, 'grle_inc', nf90_float, (/dim_lon_out,dim_lat_out,dim_lev_out/), id_grle_inc_out)
   call netcdf_err(error, 'defining variable grle_inc for file='//trim(outfile) )

   error = nf90_put_att(ncid_out, nf90_global, 'source', 'GSI')
   call netcdf_err(error, 'defining source attribute for file='//trim(outfile) )
  
   error = nf90_put_att(ncid_out, nf90_global, 'comment', 'interpolated global analysis increment')
   call netcdf_err(error, 'defining comment attribute for file='//trim(outfile) )
  
   error = nf90_enddef(ncid_out, header_buffer_val, 4,0,4)
   call netcdf_err(error, 'end meta define for file='//trim(outfile) )
 end if

!-----------------------------------------------------------------
! Compute latitude and longitude of output grid.  
!-----------------------------------------------------------------

 allocate(latitude_out(lat_out))
 allocate(slat(lat_out))
 allocate(wlat(lat_out))

 call splat(4, lat_out, slat, wlat)
 do j = 1, lat_out
   latitude_out(j) = -( 90.0_8 - (acos(slat(j))* rad2deg) )
 enddo
 deallocate(slat, wlat)

!print*,'lat out ',latitude_out(1), latitude_out(lat_out)

 allocate(longitude_out(lon_out))
 do i = 1, lon_out
   longitude_out(i) = real(i-1,8) * dlondeg
 enddo

!print*,'lon out ',longitude_out(1), longitude_out(lon_out)

!-----------------------------------------------------------------
! Compute grib 1 grid description section for output gaussian
! grid.  The GDS is required for ipolates.
!-----------------------------------------------------------------

 kgds_out = 0
 kgds_out(1) = 4            ! oct 6 - type of grid (gaussian)
 kgds_out(2) = lon_out      ! oct 7-8 - # pts on latitude circle
 kgds_out(3) = lat_out      ! oct 9-10 - # pts on longitude circle
 kgds_out(4) = nint(latitude_out(1)*1000.0_8)  ! oct 11-13 - lat of origin
 kgds_out(5) = 0            ! oct 14-16 - lon of origin
 kgds_out(6) = 128          ! oct 17 - resolution flag
 kgds_out(7) = nint(latitude_out(lat_out)*1000.0_8)     ! oct 18-20 - lat of extreme pt
 kgds_out(8) = nint(longitude_out(lon_out)*1000.0_8)    ! oct 21-23 - lon of extreme pt
 kgds_out(9) = nint((360.0_8 / real(lon_out,8))*1000.0_8) ! oct 24-25 - long increment
 kgds_out(10) = lat_out / 2 ! oct 26-27 - number of circles pole to equator
 kgds_out(11) = 64          ! oct 28 - scan mode
 kgds_out(12) = 255         ! oct 29 - reserved
 kgds_out(19) = 0           ! oct 4 - # vert coordinate parameters
 kgds_out(20) = 255         ! oct 5 - not used set to 255 (missing)

!print*,'kgds out ',kgds_out(1:20)

!----------------------------------------------------
! Open and read input file
!----------------------------------------------------

 if (mype == npes-1) print*,'- OPEN INPUT FILE: ', trim(infile)

 ! Opening for parallel access breaks global workflow on S4
 ! Also there is no parallel access, so there is no need to open for parallel
 ! access.
 !error = nf90_open(trim(infile), ior(nf90_nowrite, nf90_mpiio), &
 !                  comm=mpi_comm_world, info = mpi_info_null, ncid=ncid_in)
 error = nf90_open(trim(infile), ior(nf90_nowrite, nf90_mpiio), ncid=ncid_in)
 call netcdf_err(error, 'opening file='//trim(infile) )

 error = nf90_inq_dimid(ncid_in, 'lon', id_dim)
 call netcdf_err(error, 'inquiring lon dimension for file='//trim(infile) )
 error = nf90_inquire_dimension(ncid_in, id_dim, len=lon_in)
 call netcdf_err(error, 'reading lon dimension for file='//trim(infile) )
 allocate(longitude_in(lon_in))
 error = nf90_inq_varid(ncid_in, 'lon', id_dim)
 call netcdf_err(error, 'inquiring var lon dimension for file='//trim(infile) )
 error = nf90_get_var(ncid_in, id_dim, longitude_in)
 call netcdf_err(error, 'reading longitude_in for file='//trim(infile) )

!print*,'lon of input file is ',lon_in
!print*,'lon in ',longitude_in(1), longitude_in(lon_in)

 error = nf90_inq_dimid(ncid_in, 'lat', id_dim)
 call netcdf_err(error, 'inquiring lat dimension for file='//trim(infile) )
 error = nf90_inquire_dimension(ncid_in, id_dim, len=lat_in)
 call netcdf_err(error, 'reading lat dimension for file='//trim(infile) )
 allocate(latitude_in(lat_in))
 error = nf90_inq_varid(ncid_in, 'lat', id_dim)
 call netcdf_err(error, 'inquiring var lat dimension for file='//trim(infile) )
 error = nf90_get_var(ncid_in, id_dim, latitude_in)
 call netcdf_err(error, 'reading latitude_in for file='//trim(infile) )

!print*,'lat of input file is ',lat_in
!print*,'lat in ',latitude_in(1), latitude_in(lat_in)

 error = nf90_inq_dimid(ncid_in, 'lev', id_dim)
 call netcdf_err(error, 'inquiring lev dimension for file='//trim(infile) )
 error = nf90_inquire_dimension(ncid_in, id_dim, len=lev_in)
 call netcdf_err(error, 'reading lev dimension for file='//trim(infile) )

!print*,'lev of input file is ',lev_in

!-----------------------------------------------------------------
! Compute grib 1 grid description section for input gaussian
! grid.
!-----------------------------------------------------------------

 kgds_in = 0
 kgds_in(1) = 4            ! oct 6 - type of grid (gaussian)
 kgds_in(2) = lon_in       ! oct 7-8 - # pts on latitude circle
 kgds_in(3) = lat_in       ! oct 9-10 - # pts on longitude circle
 kgds_in(4) = nint(latitude_in(1)*1000.0_8)  ! oct 11-13 - lat of origin
 kgds_in(5) = 0            ! oct 14-16 - lon of origin
 kgds_in(6) = 128          ! oct 17 - resolution flag
 kgds_in(7) = nint(latitude_in(lat_in)*1000.0_8)  ! oct 18-20 - lat of extreme pt
 kgds_in(8) = nint(longitude_in(lon_in)*1000.0_8) ! oct 21-23 - lon of extreme pt
 kgds_in(9) = nint((360.0_8 / real(lon_in,8))*1000.0_8) ! oct 24-25 - long increment
 kgds_in(10) = lat_in / 2  ! oct 26-27 - number of circles pole to equator
 kgds_in(11) = 64          ! oct 28 - scan mode
 kgds_in(12) = 255         ! oct 29 - reserved
 kgds_in(19) = 0           ! oct 4 - # vert coordinate parameters
 kgds_in(20) = 255         ! oct 5 - not used set to 255 (missing)

!print*,'kgds in ',kgds_in(1:20)

 if (lev /= lev_in) then
   print*,'- FATAL ERROR: input and output levels dont match: ',lev_in, lev
   stop 56
 endif

!-----------------------------------------------------------------
! Loop over each record, then interpolate using ipolates.
! Interpolated data is then written to output file.
!-----------------------------------------------------------------

 mi = lon_in * lat_in
 mo = lon_out * lat_out
 allocate(dummy_out(lon_out,lat_out,lev))
 allocate(dummy_in(lon_in, lat_in, lev))
 allocate(ibi(lev))
 allocate(li(mi,lev))
 allocate(gi(mi,lev))
 allocate(gi2(mi,lev))
 allocate(rlat(mo),crot(mo))
 crot = 0; srot = 0
 allocate(rlon(mo),srot(mo))
 allocate(ibo(lev))
 allocate(lo(mo,lev))
 allocate(go(mo,lev))
 allocate(go2(mo,lev))
 allocate(go3(mo,lev))
 allocate(send_layer(mo))
 allocate(recv_layer(mo))

 call mpi_barrier(mpi_comm_world, mpierr)
 do rec = 1, num_recs

   ! skip v_inc (done with u_inc, which comes first)
   if (trim(records(rec)) .eq. 'v_inc') cycle

   if (mype == rec-1) then
     print*,'- PROCESS RECORD: ', trim(records(rec))
     readvar = .true.
  
     error = nf90_inq_varid(ncid_in, trim(records(rec)), id_var)
     ! handle missing hydrometeor increments
     if (error .ne. 0) then
       if (ANY((/ 'rwmr_inc', 'snmr_inc', 'grle_inc' /) == trim(records(rec)))) then
         print *, 'WARNING: ', trim(records(rec)), ' is missing in increment file. Skipping.'
         readvar = .false.
       else
         call netcdf_err(error, 'inquiring ' // trim(records(rec)) // ' id for file='//trim(infile) )
       end if
     end if
     if (readvar) then
       error = nf90_get_var(ncid_in, id_var, dummy_in)
       call netcdf_err(error, 'reading ' //  trim(records(rec)) // ' for file='//trim(infile) )
     else
       dummy_in(:,:,:) = 0.0
     end if
  
     ip = 0 ! bilinear
     ipopt = 0
     ibi = 0
     li = 0
     gi = 0.0_8
     rlat = 0.0_8
     rlon = 0.0_8
     ibo = 0
     lo = 0
     go = 0.0_8
     gi = reshape (dummy_in, (/mi, lev/))
  
     if (trim(records(rec)) .eq. 'u_inc') then
        ! do u_inc,v_inc at the same time
        error = nf90_inq_varid(ncid_in, 'v_inc', id_var)
        call netcdf_err(error, 'inquiring v_inc id for file='//trim(infile) )
        error = nf90_get_var(ncid_in, id_var, dummy_in)
        call netcdf_err(error, 'reading v_inc for file='//trim(infile) )
        gi2 = reshape (dummy_in, (/mi, lev/))
        call ipolatev(ip, ipopt, kgds_in, kgds_out, mi, mo,&
                      lev, ibi, li, gi, gi2,  &
                      no, rlat, rlon, crot, srot, ibo, lo, &
                      go, go3, iret)
        if (iret /= 0) then
          print*,'FATAL ERROR in ipolatev, iret: ',iret
          stop 76
        endif
        if (no /= mo) then
          print*,'FATAL ERROR: ipolatev returned wrong number of pts ',no
          stop 77
        endif

        do klev=1, lev
          send_layer=go(:,klev)
          call mpi_send(send_layer, size(send_layer), mpi_double_precision, &
                        npes-1, 1000+rec, mpi_comm_world, mpierr)
        enddo

        do klev=1, lev
          send_layer=go3(:,klev)
          call mpi_send(send_layer, size(send_layer), mpi_double_precision, &
                        npes-1, 2000+rec, mpi_comm_world, mpierr)
        enddo
     else
        call ipolates(ip, ipopt, kgds_in, kgds_out, mi, mo, &
                   lev, ibi, li, gi, no, rlat, rlon, ibo, &
                   lo, go, iret)
        if (iret /= 0) then
          print*,'FATAL ERROR in ipolates, iret: ',iret
          stop 76
        endif
        if (no /= mo) then
          print*,'FATAL ERROR: ipolates returned wrong number of pts ',no
          stop 77
        endif
        !dummy_out = reshape(go, (/lon_out,lat_out,lev/))
        !print *, lon_out, lat_out, lev, 'send'
        do klev=1, lev
          send_layer=go(:,klev)
          call mpi_send(send_layer, size(send_layer), mpi_double_precision, &
                        npes-1, 1000+rec, mpi_comm_world, mpierr)
        enddo
     endif
   else if (mype == npes-1) then
     !print *, lon_out, lat_out, lev, 'recv'
     do klev=1, lev
       call mpi_recv(recv_layer, size(recv_layer), mpi_double_precision, &
                   rec-1, 1000+rec, mpi_comm_world, mpistat, mpierr)
       go2(:,klev) = recv_layer
     enddo
     dummy_out = reshape(go2, (/lon_out,lat_out,lev/))
     error = nf90_inq_varid(ncid_out, trim(records(rec)), id_var)
     call netcdf_err(error, 'inquiring ' // trim(records(rec)) // ' id for file='//trim(outfile) )
     error = nf90_put_var(ncid_out, id_var, dummy_out)
     call netcdf_err(error, 'writing ' // trim(records(rec)) // ' for file='//trim(outfile) )
     if (trim(records(rec)) .eq. 'u_inc') then
        ! process v_inc also.
        do klev=1, lev
          call mpi_recv(recv_layer, size(recv_layer), mpi_double_precision, &
                        rec-1, 2000+rec, mpi_comm_world, mpistat, mpierr)
          go2(:,klev) = recv_layer
        enddo
        dummy_out = reshape(go2, (/lon_out,lat_out,lev/))
        error = nf90_inq_varid(ncid_out, 'v_inc', id_var)
        call netcdf_err(error, 'inquiring v_inc id for file='//trim(outfile) )
        error = nf90_put_var(ncid_out, id_var, dummy_out)
        call netcdf_err(error, 'writing v_inc for file='//trim(outfile) )
     endif
   endif 
 enddo  ! records

 error = nf90_close(ncid_in)
 call mpi_barrier(mpi_comm_world, mpierr)

 deallocate(dummy_out)
 deallocate(dummy_in)
 deallocate(ibi)
 deallocate(li)
 deallocate(gi,gi2)
 deallocate(rlat, rlon, ibo, lo, go, go2, go3, crot, srot)

!------------------------------------------------------------------
! Update remaining output file records according to Cory's sample.
!------------------------------------------------------------------

 if (mype == npes-1) then
   print*,"- WRITE OUTPUT FILE: ", trim(outfile)
  
  ! lev
  
   allocate(levs(lev))
   do j = 1, lev
     levs(j) = j
   enddo
  
   error = nf90_put_var(ncid_out, id_lev_out, levs)
   call netcdf_err(error, 'writing levs for file='//trim(outfile) )
  
  ! pfull
  
   error = nf90_put_var(ncid_out, id_pfull_out, levs)
   call netcdf_err(error, 'writing pfull for file='//trim(outfile) )
  
   deallocate (levs)
   allocate (levs(ilev))
   do j = 1, ilev
     levs(j) = j
   enddo
  
  ! ilev
  
   error = nf90_put_var(ncid_out, id_ilev_out, levs)
   call netcdf_err(error, 'writing ilev for file='//trim(outfile) )
  
  ! hyai
  
   error = nf90_put_var(ncid_out, id_hyai_out, levs)
   call netcdf_err(error, 'writing hyai for file='//trim(outfile) )
  
  ! hybi
  
   error = nf90_put_var(ncid_out, id_hybi_out, levs)
   call netcdf_err(error, 'writing hybi for file='//trim(outfile) )
  
  ! latitude
  
   error = nf90_put_var(ncid_out, id_lat_out, latitude_out)
   call netcdf_err(error, 'writing latitude for file='//trim(outfile) )
  
  ! longitude
  
   error = nf90_put_var(ncid_out, id_lon_out, longitude_out)
   call netcdf_err(error, 'writing longitude for file='//trim(outfile) )
  
   deallocate(levs)
  
   error = nf90_close(ncid_out)
 end if

 call mpi_barrier(mpi_comm_world, mpierr)
 if (mype == npes-1) print*,'- NORMAL TERMINATION'

 if (mype == npes-1) call w3tage('INTERP_INC')
 call mpi_barrier(mpi_comm_world, mpierr)
 call mpi_finalize(mpierr)

 end program interp_inc

 subroutine netcdf_err( err, string )

 use netcdf

 implicit none
 integer, intent(in) :: err
 character(len=*), intent(in) :: string
 character(len=256) :: errmsg

 if( err.EQ.NF90_NOERR )return
 errmsg = NF90_STRERROR(err)
 print*,''
 print*,'FATAL ERROR: ', trim(string), ': ', trim(errmsg)
 print*,'STOP.'
 stop 999

 return
 end subroutine netcdf_err