Bugfix for where some fields are not in the input increment file (#53)

* skip possibly missing hydrometeors
* add warning print
* handle case of missing increment in calc_analysis, add omga to calc_analysis
* correct erroneous calc_analysis runtime printout

src/netcdf_io/calc_analysis.fd/inc2anl.f90

@@ -17,7 +17,7 @@ module inc2anl
 
   integer, parameter :: nincv=13
   character(len=7) :: incvars_nemsio(nincv), incvars_netcdf(nincv), incvars_ncio(nincv)
-  integer, parameter :: nnciov=22
+  integer, parameter :: nnciov=23
   integer, parameter :: naero=14
   integer, parameter :: naero_copy=6
   character(len=7) :: iovars_netcdf(nnciov), iovars_aero(naero), copyvars_aero(naero_copy)
@@ -36,7 +36,7 @@ module inc2anl
                         'delz   ', 'dpres  ', 'dzdt   ', 'grle   ', 'hgtsfc ',&
                         'icmr   ', 'o3mr   ', 'pressfc', 'rwmr   ', 'snmr   ',&
                         'spfh   ', 'tmp    ', 'ugrd   ', 'vgrd   ', 'cld_amt',&
-                        'nccice ', 'nconrd '/
+                        'nccice ', 'nconrd ', 'omga   '/
   data iovars_aero / 'so4    ', 'bc1    ', 'bc2    ', 'oc1    ', 'oc2    ', &
                      'dust1  ', 'dust2  ', 'dust3  ', 'dust4  ', 'dust5  ',&
                      'seas1  ', 'seas2  ', 'seas3  ', 'seas4  '/
@@ -60,11 +60,11 @@ subroutine gen_anl
   !            increment, add the two together, and write out
   !            the analysis to a new file
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    use vars_calc_analysis, only: mype, do_aero
+    use vars_calc_analysis, only: mype, do_aero, jedi
     implicit none
     ! variables local to this subroutine
     integer :: i, j, iincvar
-    logical :: use_increment
+    logical :: use_increment,readinc
 
 
     ! loop through each variable in the background file
@@ -85,9 +85,15 @@ subroutine gen_anl
           call add_psfc_increment
         else
           ! call generic subroutine for all other fields
-          if (mype==0) print *, 'Adding Increment to ', iovars_netcdf(i), incvars_netcdf(iincvar)
-          call add_increment(iovars_netcdf(i), incvars_netcdf(iincvar))
-        end if
+          call add_increment(iovars_netcdf(i), incvars_netcdf(iincvar), readinc)
+          if (mype==0) then
+            if ((.not.jedi) .or. (jedi.and.readinc)) then
+               print *, 'Adding Increment to ', iovars_netcdf(i), incvars_netcdf(iincvar)
+            else
+               print *, 'Copying from Background ', iovars_netcdf(i)
+            endif
+          endif
+        endif
       else
         ! otherwise just write out what is in the input to the output
         if (mype==0) print *, 'Copying from Background ', iovars_netcdf(i)
@@ -241,7 +247,7 @@ subroutine add_aero_inc(fcstvar, incvar)
 
   end subroutine add_aero_inc
 
-  subroutine add_increment(fcstvar, incvar)
+  subroutine add_increment(fcstvar, incvar, readinc)
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! subroutine add_increment
   !            generic subroutine for adding increment to background
@@ -250,6 +256,7 @@ subroutine add_increment(fcstvar, incvar)
   !       fcstvar - input string of netCDF fcst/anal var name
   !       incvar  - input string of netCDF increment var name
   !                 (without _inc suffix added)
+  !       readinc - .true. if read increment, .false otherwise
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     use vars_calc_analysis, only: fcstncfile, anlncfile, incr_file,&
                                   nlat, nlon, nlev, anlfile, use_nemsio_anl, &
@@ -260,6 +267,7 @@ subroutine add_increment(fcstvar, incvar)
     implicit none
     ! input vars
     character(7), intent(in) :: fcstvar, incvar
+    logical, intent(out):: readinc
     ! local variables
     real, allocatable, dimension(:,:,:) :: work3d_bg
     real, allocatable, dimension(:,:) :: work3d_inc_gsi
@@ -268,6 +276,7 @@ subroutine add_increment(fcstvar, incvar)
     integer :: j,jj,k,krev,iret
     type(Dataset) :: incncfile
 
+    readinc=.true.
     if (has_var(fcstncfile, fcstvar)) then
       do k=1,nlev
         if (mype == levpe(k)) then
@@ -277,12 +286,16 @@ subroutine add_increment(fcstvar, incvar)
           incncfile = open_dataset(incr_file, paropen=.true.)
           ! JEDI-derived increments have a time dimension, so read to appropriate array
           if ( jedi ) then
-             call read_vardata(incncfile, trim(incvar)//"_inc", work3d_inc_jedi, nslice=k, slicedim=3)
-             ! add increment to background
-             do j=1,nlat
-                jj=nlat+1-j ! increment is S->N, history files are N->S
-                work3d_bg(:,j,1) = work3d_bg(:,j,1) + work3d_inc_jedi(:,jj,1)
-             end do
+             if (has_var(incncfile, trim(incvar)//"_inc")) then
+                call read_vardata(incncfile, trim(incvar)//"_inc", work3d_inc_jedi, nslice=k, slicedim=3)
+                ! add increment to background, otherwise do nothing
+                do j=1,nlat
+                   jj=nlat+1-j ! increment is S->N, history files are N->S
+                   work3d_bg(:,j,1) = work3d_bg(:,j,1) + work3d_inc_jedi(:,jj,1)
+                end do
+             else
+                readinc = .false.
+             endif
           else
              call read_vardata(incncfile, trim(incvar)//"_inc", work3d_inc_gsi, nslice=k, slicedim=3)
              ! add increment to background
@@ -306,7 +319,7 @@ subroutine add_increment(fcstvar, incvar)
       end do
       ! clean up and close
       if ( jedi ) then
-        deallocate(work3d_inc_jedi)
+        if (allocated(work3d_inc_jedi)) deallocate(work3d_inc_jedi)
       else
         deallocate(work3d_inc_gsi)
       end if

src/netcdf_io/interp_inc.fd/driver.F90

@@ -81,6 +81,7 @@ program interp_inc
  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', &
@@ -366,11 +367,24 @@ program interp_inc
 
    if (mype == rec-1) then
      print*,'- PROCESS RECORD: ', trim(records(rec))
+     readvar = .true.
 
      error = nf90_inq_varid(ncid_in, trim(records(rec)), id_var)
-     call netcdf_err(error, 'inquiring ' // trim(records(rec)) // ' id for file='//trim(infile) )
-     error = nf90_get_var(ncid_in, id_var, dummy_in)
-     call netcdf_err(error, 'reading ' //  trim(records(rec)) // ' for file='//trim(infile) )
+     ! 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