From f1271f69adb3838bf320e957b5b8060caa89c321 Mon Sep 17 00:00:00 2001
From: Eric Kemp <eric.kemp@nasa.gov>
Date: Fri, 8 Sep 2023 16:52:09 -0400
Subject: [PATCH 01/16] First attempt at adding IMERG Monthly support.

This only supports V07A Final Run.

Code compiles, but has not been run yet.
---
 .../IMERG_monthly/IMERG_monthly_dataMod.F90   | 170 +++++++++
 .../IMERG_monthly/readIMERGmonthlydata.F90    | 327 ++++++++++++++++++
 lvt/make/Filepath                             |   2 +-
 lvt/plugins/LVT_datastream_pluginMod.F90      |   9 +-
 lvt/plugins/LVT_pluginIndices.F90             |   2 +
 5 files changed, 508 insertions(+), 2 deletions(-)
 create mode 100644 lvt/datastreams/IMERG_monthly/IMERG_monthly_dataMod.F90
 create mode 100644 lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90

diff --git a/lvt/datastreams/IMERG_monthly/IMERG_monthly_dataMod.F90 b/lvt/datastreams/IMERG_monthly/IMERG_monthly_dataMod.F90
new file mode 100644
index 000000000..5321bdb3e
--- /dev/null
+++ b/lvt/datastreams/IMERG_monthly/IMERG_monthly_dataMod.F90
@@ -0,0 +1,170 @@
+!-----------------------BEGIN NOTICE -- DO NOT EDIT-----------------------
+! NASA Goddard Space Flight Center
+! Land Information System Framework (LISF)
+! Version 7.4
+!
+! Copyright (c) 2022 United States Government as represented by the
+! Administrator of the National Aeronautics and Space Administration.
+! All Rights Reserved.
+!-------------------------END NOTICE -- DO NOT EDIT-----------------------
+!
+! NOTE:  Currently only V07A IMERG Final Run data are supported
+module IMERG_monthly_dataMod
+
+  ! Defaults
+  implicit none
+  private
+
+  ! Public routines
+  public :: IMERG_monthly_datainit
+
+  ! Public types
+  public :: imergmonthlydata
+
+  type, public :: imergmonthlydatadec
+     character*100 :: odir
+     character*10  :: imergver, imergprd
+     real          :: datares
+     real, allocatable           :: rlat(:)
+     real, allocatable           :: rlon(:)
+
+     ! This is only used with upscale averaging
+     integer, allocatable        :: n11(:)
+
+     ! These are only used with budget interpolation
+     integer, allocatable        :: n112(:,:)
+     integer, allocatable        :: n122(:,:)
+     integer, allocatable        :: n212(:,:)
+     integer, allocatable        :: n222(:,:)
+     real,    allocatable        :: w112(:,:)
+     real,    allocatable        :: w122(:,:)
+     real,    allocatable        :: w212(:,:)
+     real,    allocatable        :: w222(:,:)
+
+     integer                     :: nc
+     integer                     :: nr
+  end type imergmonthlydatadec
+
+  type(imergmonthlydatadec), allocatable :: imergmonthlydata(:)
+
+contains
+
+  subroutine IMERG_monthly_datainit(i)
+
+    ! Imports
+    use ESMF
+    use LVT_coreMod, only: LVT_rc, LVT_config, LVT_isAtAFinerResolution
+    use LVT_logMod, only: LVT_logunit, LVT_verify, LVT_endrun
+
+    ! Defaults
+    implicit none
+
+    ! Arguments
+    integer, intent(in) :: i
+
+    ! Local variables
+    integer              :: status
+    real                 :: gridDesci(50)
+
+    ! External routines
+    external :: conserv_interp_input
+    external :: upscaleByAveraging_input
+
+    if (.not. allocated(imergmonthlydata)) then
+       allocate(imergmonthlydata(LVT_rc%nDataStreams))
+    endif
+
+    ! Get top level IMERG data directory
+    call ESMF_ConfigGetAttribute(LVT_Config, imergmonthlydata(i)%odir, &
+         label='IMERG monthly data directory:', rc=status)
+    call LVT_verify(status, 'IMERG monthly data directory: not defined')
+
+    ! Get IMERG product
+    call ESMF_ConfigFindLabel(LVT_config," IMERG monthly product:", &
+         rc=status)
+    call LVT_verify(status, 'IMERG monthly product: not defined')
+    call ESMF_ConfigGetAttribute(LVT_config, &
+         imergmonthlydata(i)%imergprd, rc=status)
+    if (imergmonthlydata(i)%imergprd .ne. 'final') then
+       write(LVT_logunit,*)'[ERR] Invalid IMERG monthly product specified'
+       write(LVT_logunit,*)"[ERR] Currently only 'final' is supported"
+       call LVT_endrun()
+    end if
+
+    ! Get IMERG version
+    call ESMF_ConfigFindLabel(LVT_config, "IMERG monthly version:", &
+         rc=status)
+    call LVT_verify(status, 'IMERG monthly version: not defined')
+    call ESMF_ConfigGetAttribute(LVT_config, &
+         imergmonthlydata(i)%imergver, rc=status)
+    if (imergmonthlydata(i)%imergver .ne. 'V07A') then
+       write(LVT_logunit,*)'[ERR] Invalid IMERG monthly version specified'
+       write(LVT_logunit,*)"[ERR] Currently only 'V07A' is supported"
+       call LVT_endrun()
+    end if
+
+    ! Allocate arrays on LVT grid
+    allocate(imergmonthlydata(i)%rlat(LVT_rc%lnc*LVT_rc%lnr))
+    allocate(imergmonthlydata(i)%rlon(LVT_rc%lnc*LVT_rc%lnr))
+
+    ! Set IMERG grid and map projection information.
+    gridDesci(:) =    0
+    gridDesci(1) =    0     ! Lat/lon
+    gridDesci(2) = 3600     ! Points along latitude circle
+    gridDesci(3) = 1800     ! Points along longitude circle
+    gridDesci(4) =  -89.95  ! Latitude of first grid point
+    gridDesci(5) = -179.95  ! Longitude of first grid point
+    gridDesci(6) =  128
+    gridDesci(7) =   89.95  ! Latitude of last grid point
+    gridDesci(8) =  179.95  ! Longitude of last grid point
+    gridDesci(9) =    0.1   ! Longitudinal direction increment
+    gridDesci(10) =   0.1   ! Latitude direction increment
+    gridDesci(20) =  64
+
+    ! Set up interpolation data
+    imergmonthlydata(i)%datares = 0.1
+    imergmonthlydata(i)%nc = 3600
+    imergmonthlydata(i)%nr = 1800
+
+    ! Use budget-bilinear interpolation if IMERG resolution (0.1 deg)
+    ! is coarser than the analysis grid. Use upscale averaging if
+    ! IMERG is finer than the analysis grid.
+    if (LVT_isAtAFinerResolution(imergmonthlydata(i)%datares)) then
+
+       ! Used only with budget interpolation
+       allocate(imergmonthlydata(i)%n112(LVT_rc%lnc*LVT_rc%lnr,25))
+       allocate(imergmonthlydata(i)%n122(LVT_rc%lnc*LVT_rc%lnr,25))
+       allocate(imergmonthlydata(i)%n212(LVT_rc%lnc*LVT_rc%lnr,25))
+       allocate(imergmonthlydata(i)%n222(LVT_rc%lnc*LVT_rc%lnr,25))
+       allocate(imergmonthlydata(i)%w112(LVT_rc%lnc*LVT_rc%lnr,25))
+       allocate(imergmonthlydata(i)%w122(LVT_rc%lnc*LVT_rc%lnr,25))
+       allocate(imergmonthlydata(i)%w212(LVT_rc%lnc*LVT_rc%lnr,25))
+       allocate(imergmonthlydata(i)%w222(LVT_rc%lnc*LVT_rc%lnr,25))
+       imergmonthlydata(i)%n112 = 0
+       imergmonthlydata(i)%n122 = 0
+       imergmonthlydata(i)%n212 = 0
+       imergmonthlydata(i)%n222 = 0
+       imergmonthlydata(i)%w112 = 0
+       imergmonthlydata(i)%w122 = 0
+       imergmonthlydata(i)%w212 = 0
+       imergmonthlydata(i)%w222 = 0
+       call conserv_interp_input(gridDesci, LVT_rc%gridDesc, &
+            LVT_rc%lnc*LVT_rc%lnr, &
+            imergmonthlydata(i)%rlat, imergmonthlydata(i)%rlon, &
+            imergmonthlydata(i)%n112, imergmonthlydata(i)%n122, &
+            imergmonthlydata(i)%n212, imergmonthlydata(i)%n222, &
+            imergmonthlydata(i)%w112, imergmonthlydata(i)%w122, &
+            imergmonthlydata(i)%w212, imergmonthlydata(i)%w222)
+    else
+
+       ! Used only with upscale averaging
+       allocate(imergmonthlydata(i)%n11(imergmonthlydata(i)%nc * &
+            imergmonthlydata(i)%nr))
+       imergmonthlydata(i)%n11 = 0
+       call upscaleByAveraging_input(gridDesci, LVT_rc%gridDesc,&
+            imergmonthlydata(i)%nc*imergmonthlydata(i)%nr, &
+            LVT_rc%lnc*LVT_rc%lnr, &
+            imergmonthlydata(i)%n11)
+    end if
+  end subroutine IMERG_monthly_datainit
+end module IMERG_monthly_dataMod
diff --git a/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90 b/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90
new file mode 100644
index 000000000..1ee98dd7b
--- /dev/null
+++ b/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90
@@ -0,0 +1,327 @@
+!-----------------------BEGIN NOTICE -- DO NOT EDIT-----------------------
+! NASA Goddard Space Flight Center
+! Land Information System Framework (LISF)
+! Version 7.4
+!
+! Copyright (c) 2022 United States Government as represented by the
+! Administrator of the National Aeronautics and Space Administration.
+! All Rights Reserved.
+!-------------------------END NOTICE -- DO NOT EDIT-----------------------
+
+#include "LVT_misc.h"
+
+!NOTE:  Currently only V07A IMERG Final Run monthly data are supported.
+subroutine readIMERGmonthlydata(source)
+
+  ! Imports
+  use ESMF
+  use IMERG_monthly_dataMod, only: imergmonthlydata
+  use LVT_coreMod, only: LVT_rc, LVT_isAtAFinerResolution
+  use LVT_histDataMod, only: LVT_logSingleDataStreamVar, &
+       LVT_MOC_totalprecip
+  use LVT_logMod, only: LVT_verify, LVT_logunit
+  use LVT_timeMgrMod, only: LVT_calendar
+
+  ! Defaults
+  implicit none
+
+  ! Arguments
+  integer, intent(in) :: source
+
+  ! Local variables
+  character*255       :: filename
+  logical             :: file_exists
+  real                :: prcp_in(imergmonthlydata(source)%nc, &
+       imergmonthlydata(source)%nr)
+  real                :: prcp_in1(imergmonthlydata(source)%nc * &
+       imergmonthlydata(source)%nr)
+  logical*1           :: lb(imergmonthlydata(source)%nc * &
+       imergmonthlydata(source)%nr)
+  logical*1           :: lo(LVT_rc%lnc*LVT_rc%lnr)
+  real                :: prcp(LVT_rc%lnc*LVT_rc%lnr)
+  real                :: prcp_final(LVT_rc%lnc, LVT_rc%lnr)
+  integer             :: t, c, r
+  integer             :: iret, ireaderr
+  real                :: currTime
+  logical             :: alarmCheck
+  integer             :: yr1, mo1, da1, hr1, mn1, ss1
+  integer             :: yr2, mo2, da2, hr2, mn2, ss2
+  type(ESMF_Time)     :: time1
+  type(ESMF_Time)     :: time2
+  type(ESMF_TimeInterval) :: lis_ts
+  integer :: status
+  integer :: days_in_month
+
+  ! External routines
+  external :: create_IMERG_monthly_filename
+  external :: read_imergmonthly_hdf
+  external :: conserv_interp
+  external :: upscaleByAveraging
+
+  ! Initialize variables
+  prcp(:) = LVT_rc%udef
+  prcp_final(:,:) = LVT_rc%udef
+  currTime = float(LVT_rc%dhr(source))*3600 + &
+       60*LVT_rc%dmn(source) + LVT_rc%dss(source)
+
+  ! Read IMERG at beginning of the month
+  alarmCheck = .false.
+  if (LVT_rc%dda(source) == 1 .and. &
+       LVT_rc%dhr(source) == 0 .and. &
+       LVT_rc%dmn(source) == 0 .and. &
+       LVT_rc%dss(source) == 0) alarmCheck = .true.
+
+  yr1 = LVT_rc%dyr(source)
+  mo1 = LVT_rc%dmo(source)
+  da1 = LVT_rc%dda(source)
+  hr1 = LVT_rc%dhr(source)
+  mn1 = LVT_rc%dmn(source)
+  ss1 = 0
+  call ESMF_TimeSet(time1,yy=yr1, mm=mo1, dd=da1, &
+       h=hr1,m=mn1,s=ss1,calendar=LVT_calendar,rc=status)
+  call LVT_verify(status)
+
+  ! Use previous month file.
+  call ESMF_TimeIntervalSet(lis_ts, mm=1, calendar=LVT_calendar, &
+       rc=status)
+  call LVT_verify(status)
+  time2 = time1 - lis_ts
+  call ESMF_TimeGet(time2, yy=yr2, mm=mo2, dd=da2, &
+        h=hr2, m=mn2, s=ss2, calendar=LVT_calendar, &
+        rc=status)
+  call LVT_verify(status)
+
+  if (alarmCheck) then
+     call create_IMERG_monthly_filename(source, yr1, mo1, filename)
+     inquire(file=trim(filename), exist=file_exists)
+
+     if (file_exists) then
+        write(LVT_logunit,*) '[INFO] Reading IMERG data ',trim(filename)
+        call read_imergmonthly_hdf(filename, &
+             imergmonthlydata(source)%nc, &
+             imergmonthlydata(source)%nr, prcp_in, ireaderr)
+        if (ireaderr .eq. 0) then
+
+           ! Use budget-bilinear interpolation if IMERG data are at
+           ! coarser resolution than the analysis grid; otherwise, use
+           ! upscale averaging.
+           prcp_in1(:) = 0
+           lb(:) = .false.
+           t = 1
+           do r = 1, imergmonthlydata(source)%nr
+              do c = 1,imergmonthlydata(source)%nc
+                 prcp_in1(t) = prcp_in(c,r)
+                 if (prcp_in1(t) .ge. 0) then
+                    lb(t) = .true.
+                 end if
+                 t = t + 1
+              end do ! c
+           end do ! r
+
+           if (LVT_isAtAFinerResolution( &
+                imergmonthlydata(source)%datares)) then
+              call conserv_interp(LVT_rc%gridDesc, lb, prcp_in1, &
+                   lo, prcp, &
+                   (imergmonthlydata(source)%nc * &
+                   imergmonthlydata(source)%nr), &
+                   (LVT_rc%lnc*LVT_rc%lnr), &
+                   imergmonthlydata(source)%rlat, &
+                   imergmonthlydata(source)%rlon, &
+                   imergmonthlydata(source)%w112, &
+                   imergmonthlydata(source)%w122, &
+                   imergmonthlydata(source)%w212, &
+                   imergmonthlydata(source)%w222, &
+                   imergmonthlydata(source)%n112, &
+                   imergmonthlydata(source)%n122, &
+                   imergmonthlydata(source)%n212, &
+                   imergmonthlydata(source)%n222, &
+                   LVT_rc%udef, iret)
+           else
+              call upscaleByAveraging( &
+                   (imergmonthlydata(source)%nc * &
+                   imergmonthlydata(source)%nr), &
+                   (LVT_rc%lnc*LVT_rc%lnr), LVT_rc%udef, &
+                   imergmonthlydata(source)%n11, lb, &
+                   prcp_in1, lo, prcp)
+           endif
+           write(LVT_logunit,*) '[INFO] Finished processing ', &
+                trim(filename)
+        else
+           write(LVT_logunit,*) '[ERR] Read error with IMERG file ', &
+                trim(filename)
+           prcp = LVT_rc%udef
+        endif
+     else
+        write(LVT_logunit,*) '[ERR] Missing IMERG file ', trim(filename)
+        prcp = LVT_rc%udef
+     end if ! file_exists
+
+     do r = 1, LVT_rc%lnr
+        do c = 1, LVT_rc%lnc
+           prcp_final(c,r) = prcp(c+(r-1)*LVT_rc%lnc)
+        end do ! c
+     end do ! r
+  end if ! alarmCheck
+
+  ! Convert mm/hr to kg/m2s.
+  do r = 1, LVT_rc%lnr
+     do c = 1, LVT_rc%lnc
+        if (prcp_final(c,r) .ge. 0) then
+           prcp_final(c,r) = prcp_final(c,r) / 3600.
+        else
+           prcp_final(c,r) = LVT_rc%udef
+        end if
+     end do ! c
+  end do ! r
+  call LVT_logSingleDataStreamVar(LVT_MOC_totalprecip, source, &
+       prcp_final, vlevel=1, units='kg/m2s')
+
+  ! Now convert from kg/m2s to kg/m2.  Need to calculate length of
+  ! month in seconds.
+  lis_ts = time1 - time2
+  call ESMF_TimeIntervalGet(lis_ts, d=days_in_month, rc=status)
+  call LVT_verify(status)
+  write(LVT_logunit,*)'EMK: days_in_month = ', days_in_month
+  do r = 1, LVT_rc%lnr
+     do c = 1, LVT_rc%lnc
+        if (prcp_final(c,r) .ge. 0) then
+           prcp_final(c,r) = &
+                prcp_final(c,r) * days_in_month * 86400 ! kg/m2
+        else
+           prcp_final(c,r) = LVT_rc%udef
+        endif
+     enddo ! c
+  enddo ! r
+  call LVT_logSingleDataStreamVar(LVT_MOC_totalprecip, source, &
+       prcp_final, vlevel=1, units='kg/m2')
+
+end subroutine readIMERGmonthlydata
+
+subroutine read_imergmonthly_hdf(filename, col, row, precipout, ireaderr)
+
+  ! Imports
+  use LVT_logMod, only: LVT_logunit
+
+#if (defined USE_HDF5)
+  use HDF5
+#endif
+
+  ! Defaults
+  implicit none
+
+  ! Arguments
+  character(len=*), intent(in) :: filename
+  integer, intent(in)          :: col, row
+  integer, intent(out)         :: ireaderr
+  real,    intent(out)         :: precipout(col,row)
+
+  ! Local variables
+  integer :: xsize, ysize
+  character(len=40) :: dsetname='/Grid/precipitation'
+  real :: precipin(row,col)
+  integer :: istatus
+#if (defined USE_HDF5)
+  integer(HSIZE_T), dimension(2) :: dims
+  integer(HID_T) :: fileid, dsetid
+#endif
+
+#if (defined USE_HDF5)
+  xsize = col
+  ysize = row
+  dims(1) = xsize
+  dims(2) = ysize
+
+  ireaderr = 0
+
+  ! Open Fortran interface
+  call h5open_f(istatus)
+  if (istatus .ne. 0) then
+     write(LVT_logunit,*) 'Error opening HDF5 fortran interface'
+     ireaderr = istatus
+     return
+  end if
+
+  ! Open HDF5 file
+  call h5fopen_f(filename, H5F_ACC_RDONLY_F, fileid, istatus)
+  if (istatus .ne. 0) then
+     write(LVT_logunit,*) 'Error opening IMERG file', trim(filename)
+     ireaderr = istatus
+     call h5close_f(istatus) ! Close HDF5 interface
+     return
+  end if
+
+  ! Open dataset
+  call h5dopen_f(fileid, dsetname, dsetid, istatus)
+  if (istatus .ne. 0) then
+     write(LVT_logunit,*) 'Error opening IMERG dataset', trim(dsetname)
+     ireaderr = istatus
+     call h5fclose_f(fileid, istatus) ! Close HDF5 file
+     call h5close_f(istatus) ! Close HDF5 interface
+     return
+  end if
+
+  ! Read dataset
+  call h5dread_f(dsetid, H5T_NATIVE_REAL, precipin, dims, istatus)
+  if (istatus .ne. 0) then
+     write(LVT_logunit,*) 'Error reading IMERG dataset', trim(dsetname)
+     ireaderr = istatus
+     call h5dclose_f(dsetid, istatus) ! Close dataset
+     call h5fclose_f(fileid, istatus) ! Close HDF5 file
+     call h5close_f(istatus) ! Close HDF5 interface
+     return
+  end if
+
+  ! Put the real(1:,1:) on the precipout(0:,0:)
+  ! precipin is (ysize,xsize) starting at (lon=-179.9,lat=-89.9)
+  precipout(1:xsize,1:ysize) = transpose(precipin)
+
+  ! Clean up.  Since the data have already been copied, we will assume
+  ! any HDF5 errors from this point will have no impact on the data, and
+  ! therefore we will ignore the status codes.
+  call h5dclose_f(dsetid, istatus) ! Close HDF5 dataset
+  call h5fclose_f(fileid, istatus) ! Close HDF5 file
+  call h5close_f(istatus) ! Close HDF5 interface
+
+#endif
+
+end subroutine read_imergmonthly_hdf
+
+subroutine create_IMERG_monthly_filename(source, yr, mo, filename)
+
+  ! Imports
+  use IMERG_monthly_dataMod, only: imergmonthlydata
+  use LVT_logMod, only: LVT_logunit, LVT_endrun
+
+  ! Defaults
+  implicit none
+
+  ! Arguments
+  integer, intent(in) :: source, yr, mo
+  character(len=*), intent(out) :: filename
+
+  ! Local variables
+  character*4   :: cyr
+  character*2   :: cmo
+  character*100 :: fstem, fext
+  character*255 :: odir
+  character*4   :: imVer
+
+  write(cyr, '(I4.4)') yr
+  write(cmo, '(I2.2)') mo
+  odir = trim(imergmonthlydata(source)%odir)
+  imVer = imergmonthlydata(source)%imergver
+
+  ! FIXME:  Add support for Early and Late Runs
+  if (imergmonthlydata(source)%imergprd == 'final') then
+     fstem = '/3B-MO.MS.MRG.3IMERG.'
+     fext = '.HDF5'
+  else
+     write(LVT_logunit,*) "[ERR] Invalid IMERG product option was chosen."
+     write(LVT_logunit,*) "[ERR] Please choose either 'final'."
+     call LVT_endrun()
+  endif
+
+  filename = trim(odir) // "/" // cyr // cmo // trim(fstem) // &
+       cyr // cmo // "01-S000000-E235959." // cmo // trim(imVer) // fext
+
+end subroutine create_IMERG_monthly_filename
diff --git a/lvt/make/Filepath b/lvt/make/Filepath
index d0034d33b..a1e42cc42 100644
--- a/lvt/make/Filepath
+++ b/lvt/make/Filepath
@@ -1 +1 @@
-dirs := .  ../main ../core ../metrics ../plugins ../interp ../lib/bil ../domains/UTM ../domains/lambert ../domains/latlon ../datastreams/template ../datastreams/LISout ../datastreams/LISDAobs ../datastreams/ISCCP_Tskin ../datastreams/CEOP ../datastreams/SCAN ../datastreams/SCANGMAO ../datastreams/NASMD ../datastreams/GHCN ../datastreams/SNOTEL ../datastreams/SURFRAD ../datastreams/WGPBMR ../datastreams/LSWG_Tb ../datastreams/FMI_SWE ../datastreams/FMI_SNWD ../datastreams/CMC_SNWD ../datastreams/SNODAS ../datastreams/NASA_AMSRE_sm ../datastreams/LPRM_AMSRE_sm ../datastreams/AMMA ../datastreams/Ameriflux ../datastreams/ARM ../datastreams/SMOSREX ../datastreams/AGRMET ../datastreams/GlobSnow ../datastreams/SNODEPmetobs ../datastreams/SNODEP ../datastreams/MOD10A1 ../datastreams/MODSCAG ../datastreams/MOD10A1V6 ../datastreams/ANSA_SNWD ../datastreams/ANSA_SWE ../datastreams/CPC_PRCP ../datastreams/USGS_streamflow_gridded ../datastreams/USGS_streamflow ../datastreams/Natural_streamflow ../datastreams/ISMN ../datastreams/FLUXNETmte ../datastreams/FLUXNET2015 ../datastreams/FLUXNET2015_NC ../datastreams/MOD16A2 ../datastreams/UWET ../datastreams/ARSsm ../datastreams/NLDAS2 ../datastreams/ALEXI ../datastreams/ALEXIesi ../datastreams/GRACE ../datastreams/simGRACE ../datastreams/USGS_GWwell ../datastreams/PBOH2O ../datastreams/SMOS_L2sm ../datastreams/SMOS_NESDIS ../datastreams/SMOS_L1TB ../datastreams/SMOS_CATDS_L3sm ../datastreams/GCOMW_AMSR2L3sm ../datastreams/GCOMW_AMSR2L3snd ../datastreams/SMOPS ../datastreams/MODIS_LST ../datastreams/GLERL_HydroData ../datastreams/JULESdata ../datastreams/JULES2Ddata ../datastreams/ESACCI_sm ../datastreams/GIMMS_AVHRR_NDVI ../datastreams/GIMMS_MODIS_NDVI ../datastreams/GLDAS1 ../datastreams/GLDAS2 ../datastreams/MERRA2 ../datastreams/MERRA2asm ../datastreams/MERRA-Land ../datastreams/SSEBop ../datastreams/GRDC ../datastreams/GOES_LST ../datastreams/ERAinterimLand ../datastreams/SMAPsm ../datastreams/SMAPvod ../datastreams/LPRMvod ../datastreams/SMAPvwc ../datastreams/SMAP_L3TB ../datastreams/SMAPTB ../datastreams/GOME2_SIF ../datastreams/LVTbenchmarkOUT ../datastreams/LIS6out ../datastreams/LISDAdiag ../datastreams/Daymet ../datastreams/CMORPH ../datastreams/CHIRPSv2 ../datastreams/3B42V7 ../datastreams/USCRNsm ../datastreams/GLEAM ../datastreams/USDM ../datastreams/LVTpercentile ../datastreams/IMD_PRCP ../datastreams/APHRO_PRCP ../datastreams/GLASSlai ../datastreams/MODISsportLAI ../datastreams/FLUXCOM ../datastreams/HAR ../datastreams/OCO2_SIF ../datastreams/ECMWFforc ../datastreams/GDASforc ../datastreams/ASO_SWE ../runmodes/DataComp ../runmodes/557post ../runmodes/LISpost ../runmodes/DAstats ../runmodes/optUE ../runmodes/Benchmarking ../runmodes/USAFSIpost ../training/LinearRegression/ ../datastreams/GLASSalbedo ../datastreams/IMERG ../datastreams/UA_SNOW ../datastreams/OzFlux ../datastreams/jasmin ../datastreams/MCD15A2H ../datastreams/ERA5 ../datastreams/FluxSat_GPP ../datastreams/THySM/ ../datastreams/GRUNrunoff ../datastreams/UASMAP/ ../datastreams/COAMPSout/ ../datastreams/SMAP_E_OPL/
+dirs := .  ../main ../core ../metrics ../plugins ../interp ../lib/bil ../domains/UTM ../domains/lambert ../domains/latlon ../datastreams/template ../datastreams/LISout ../datastreams/LISDAobs ../datastreams/ISCCP_Tskin ../datastreams/CEOP ../datastreams/SCAN ../datastreams/SCANGMAO ../datastreams/NASMD ../datastreams/GHCN ../datastreams/SNOTEL ../datastreams/SURFRAD ../datastreams/WGPBMR ../datastreams/LSWG_Tb ../datastreams/FMI_SWE ../datastreams/FMI_SNWD ../datastreams/CMC_SNWD ../datastreams/SNODAS ../datastreams/NASA_AMSRE_sm ../datastreams/LPRM_AMSRE_sm ../datastreams/AMMA ../datastreams/Ameriflux ../datastreams/ARM ../datastreams/SMOSREX ../datastreams/AGRMET ../datastreams/GlobSnow ../datastreams/SNODEPmetobs ../datastreams/SNODEP ../datastreams/MOD10A1 ../datastreams/MODSCAG ../datastreams/MOD10A1V6 ../datastreams/ANSA_SNWD ../datastreams/ANSA_SWE ../datastreams/CPC_PRCP ../datastreams/USGS_streamflow_gridded ../datastreams/USGS_streamflow ../datastreams/Natural_streamflow ../datastreams/ISMN ../datastreams/FLUXNETmte ../datastreams/FLUXNET2015 ../datastreams/FLUXNET2015_NC ../datastreams/MOD16A2 ../datastreams/UWET ../datastreams/ARSsm ../datastreams/NLDAS2 ../datastreams/ALEXI ../datastreams/ALEXIesi ../datastreams/GRACE ../datastreams/simGRACE ../datastreams/USGS_GWwell ../datastreams/PBOH2O ../datastreams/SMOS_L2sm ../datastreams/SMOS_NESDIS ../datastreams/SMOS_L1TB ../datastreams/SMOS_CATDS_L3sm ../datastreams/GCOMW_AMSR2L3sm ../datastreams/GCOMW_AMSR2L3snd ../datastreams/SMOPS ../datastreams/MODIS_LST ../datastreams/GLERL_HydroData ../datastreams/JULESdata ../datastreams/JULES2Ddata ../datastreams/ESACCI_sm ../datastreams/GIMMS_AVHRR_NDVI ../datastreams/GIMMS_MODIS_NDVI ../datastreams/GLDAS1 ../datastreams/GLDAS2 ../datastreams/MERRA2 ../datastreams/MERRA2asm ../datastreams/MERRA-Land ../datastreams/SSEBop ../datastreams/GRDC ../datastreams/GOES_LST ../datastreams/ERAinterimLand ../datastreams/SMAPsm ../datastreams/SMAPvod ../datastreams/LPRMvod ../datastreams/SMAPvwc ../datastreams/SMAP_L3TB ../datastreams/SMAPTB ../datastreams/GOME2_SIF ../datastreams/LVTbenchmarkOUT ../datastreams/LIS6out ../datastreams/LISDAdiag ../datastreams/Daymet ../datastreams/CMORPH ../datastreams/CHIRPSv2 ../datastreams/3B42V7 ../datastreams/USCRNsm ../datastreams/GLEAM ../datastreams/USDM ../datastreams/LVTpercentile ../datastreams/IMD_PRCP ../datastreams/APHRO_PRCP ../datastreams/GLASSlai ../datastreams/MODISsportLAI ../datastreams/FLUXCOM ../datastreams/HAR ../datastreams/OCO2_SIF ../datastreams/ECMWFforc ../datastreams/GDASforc ../datastreams/ASO_SWE ../runmodes/DataComp ../runmodes/557post ../runmodes/LISpost ../runmodes/DAstats ../runmodes/optUE ../runmodes/Benchmarking ../runmodes/USAFSIpost ../training/LinearRegression/ ../datastreams/GLASSalbedo ../datastreams/IMERG ../datastreams/IMERG_monthly ../datastreams/UA_SNOW ../datastreams/OzFlux ../datastreams/jasmin ../datastreams/MCD15A2H ../datastreams/ERA5 ../datastreams/FluxSat_GPP ../datastreams/THySM/ ../datastreams/GRUNrunoff ../datastreams/UASMAP/ ../datastreams/COAMPSout/ ../datastreams/SMAP_E_OPL/
diff --git a/lvt/plugins/LVT_datastream_pluginMod.F90 b/lvt/plugins/LVT_datastream_pluginMod.F90
index 1016c4e7b..9259286b7 100644
--- a/lvt/plugins/LVT_datastream_pluginMod.F90
+++ b/lvt/plugins/LVT_datastream_pluginMod.F90
@@ -172,6 +172,7 @@ subroutine LVT_datastream_plugin
     use GDASforc_dataMod,       only : GDASforc_datainit    
     use ASOSWE_obsMod,          only : ASOSWE_obsinit
     use IMERG_dataMod,          only : IMERG_datainit
+    use IMERG_monthly_dataMod,  only : IMERG_monthly_datainit
     use UASNOW_obsMod,          only : UASNOW_obsinit
     use OzFlux_obsMod,          only : OzFlux_obsinit
     use JASMINsm_obsMod,        only : JASMINsm_obsInit
@@ -284,6 +285,7 @@ subroutine LVT_datastream_plugin
     external readGDASforcdata
     external readASOSWEObs
     external readIMERGdata
+    external readIMERGmonthlydata
     external readUASNOWObs
     external readOzFluxObs
     external readJASMINsmobs
@@ -700,7 +702,12 @@ subroutine LVT_datastream_plugin
          readASOSWEObs)
 
     call registerobssetup(trim(LVT_IMERGdataId)//char(0), IMERG_datainit)
-    call registerobsread(trim(LVT_IMERGdataId)//char(0) , readIMERGdata)
+    call registerobsread(trim(LVT_IMERGdataId)//char(0), readIMERGdata)
+
+    call registerobssetup(trim(LVT_IMERGmonthlydataId)//char(0), &
+         IMERG_monthly_datainit)
+    call registerobsread(trim(LVT_IMERGdataId)//char(0), &
+         readIMERGmonthlydata)
 
     call registerobssetup(trim(LVT_UASNOWdataId)//char(0), UASNOW_obsinit)
     call registerobsread(trim(LVT_UASNOWdataId)//char(0) , readUASNOWObs)
diff --git a/lvt/plugins/LVT_pluginIndices.F90 b/lvt/plugins/LVT_pluginIndices.F90
index e188f6197..e596e60eb 100644
--- a/lvt/plugins/LVT_pluginIndices.F90
+++ b/lvt/plugins/LVT_pluginIndices.F90
@@ -262,6 +262,8 @@ module LVT_pluginIndices
    character*50, public,  parameter :: LVT_GDASdataId = "GDAS"
    character*50, public,  parameter :: LVT_ASOSWEdataId = "ASO SWE"
    character*50, public,  parameter :: LVT_IMERGdataId = "GPM IMERG"
+   character*50, public,  parameter :: LVT_IMERGMonthlydataId = &
+        "GPM IMERG Monthly"
    character*50, public,  parameter :: LVT_UASNOWdataId = "UA SNOW"
    character*50, public,  parameter :: LVT_ozFluxdataId = "OzFlux"
    character*50, public,  parameter :: LVT_JASMINsmobsId = "JASMIN soil moisture"

From 741f6c9c50221212f3cdb9bb028d8619f847688f Mon Sep 17 00:00:00 2001
From: Eric Kemp <eric.kemp@nasa.gov>
Date: Wed, 13 Sep 2023 16:57:39 -0400
Subject: [PATCH 02/16] First cut at adding NRT bias-correction for GFS/GALWEM
 precip.

Code compiles, but not tested yet.
---
 lis/metforcing/usaf/AGRMET_forcingMod.F90 |  23 ++-
 lis/metforcing/usaf/USAF_bratsethMod.F90  | 209 +++++++++++++++++++---
 lis/metforcing/usaf/readcrd_agrmet.F90    |  53 +++++-
 3 files changed, 256 insertions(+), 29 deletions(-)

diff --git a/lis/metforcing/usaf/AGRMET_forcingMod.F90 b/lis/metforcing/usaf/AGRMET_forcingMod.F90
index 4a707ac24..ad1f1f809 100644
--- a/lis/metforcing/usaf/AGRMET_forcingMod.F90
+++ b/lis/metforcing/usaf/AGRMET_forcingMod.F90
@@ -719,6 +719,12 @@ module AGRMET_forcingMod
      ! EMK Add GFS-to-GALWEM bias correction
      integer :: back_bias_corr
      real, allocatable :: pcp_back_bias_ratio(:,:)
+     ! EMK Add NRT bias correction toward IMERG-Final Run
+     ! (back_bias_corr == 2)
+     character(255) :: gfs_nrt_bias_ratio_file
+     character(255) :: galwem_nrt_bias_ratio_file
+     real, allocatable :: gfs_nrt_bias_ratio(:,:)
+     real, allocatable :: galwem_nrt_bias_ratio(:,:)
      integer :: pcp_back_bias_ratio_month
   end type agrmet_type_dec
 
@@ -747,7 +753,7 @@ module AGRMET_forcingMod
   subroutine init_AGRMET(findex)
 
 ! !USES: 
-    use LIS_coreMod,    only : LIS_rc, LIS_domain
+    use LIS_coreMod,    only : LIS_rc
     use LIS_timeMgrMod, only : LIS_update_timestep
     use LIS_logMod,     only : LIS_logunit, LIS_endrun
     use LIS_snowMod,    only : LIS_snow_setup
@@ -827,6 +833,21 @@ subroutine init_AGRMET(findex)
     integer                       :: istat1
     integer                       :: ftn
 
+    external :: readcrd_agrmet
+    external :: AGRMET_readmask
+    external :: AGRMET_readterrain
+    external :: AGRMET_readpcpcntm
+    external :: AGRMET_read_sfcalccntm
+    external :: polarToLatLon
+    external :: bilinear_interp_input_withgrid
+    external :: neighbor_interp_input_withgrid
+    external :: AGRMET_read_pcpclimodata
+    external :: gfs_reset_interp_input
+    external :: galwem_reset_interp_input
+    external :: conserv_interp_input
+    external :: bilinear_interp_input
+    external :: upscaleByAveraging_input
+
     allocate(agrmet_struc(LIS_rc%nnest))
     call readcrd_agrmet()
 
diff --git a/lis/metforcing/usaf/USAF_bratsethMod.F90 b/lis/metforcing/usaf/USAF_bratsethMod.F90
index b6dedf826..3995698ec 100644
--- a/lis/metforcing/usaf/USAF_bratsethMod.F90
+++ b/lis/metforcing/usaf/USAF_bratsethMod.F90
@@ -628,6 +628,8 @@ subroutine USAF_addSSMIObsData(this,imax,jmax,ra_tmp,nest)
       integer :: count_good_ssmi
       real, parameter :: FILL = -9999.0
 
+      external :: polarToLatLon
+
       net = "SSMI"
       platform = "SSMI"
       sigmaOSqr = agrmet_struc(nest)%bratseth_precip_ssmi_sigma_o_sqr
@@ -774,6 +776,10 @@ subroutine USAF_getBackNWP(nest,back4,pcp_src, use_twelve, j6hr, findex)
       character(len=10) :: yyyymmddhh
       integer :: c, r
 
+      external :: AGRMET_julhr_date10
+      external :: MPI_Barrier
+      external :: sleep
+
       TRACE_ENTER("bratseth_getBackNWP")
       rc = 0
 
@@ -801,6 +807,8 @@ subroutine USAF_getBackNWP(nest,back4,pcp_src, use_twelve, j6hr, findex)
          ! EMK...Fetch bias ratio if requested
          if (agrmet_struc(nest)%back_bias_corr .eq. 1) then
             call USAF_pcpBackBiasRatio_s2s(nest, yyyymmddhh)
+         else if (agrmet_struc(nest)%back_bias_corr .eq. 2) then
+            call USAF_pcpBackBiasRatio_nrt(nest, yyyymmddhh)
          end if
 
          write(LIS_logunit,*) &
@@ -824,6 +832,20 @@ subroutine USAF_getBackNWP(nest,back4,pcp_src, use_twelve, j6hr, findex)
             ! Use GALWEM Bratseth settings if we have the data
             if (rc .eq. 0) then
                pcp_src(k) = 'GALWEM'
+
+               ! Apply bias correction to background field
+               if (agrmet_struc(nest)%back_bias_corr .eq. 2) then
+                  write(LIS_logunit,*) &
+           '[INFO] Applying IMERG-based bias correction to GALWEM precip'
+                  do r = 1, LIS_rc%gnr(nest)
+                     do c = 1, LIS_rc%gnc(nest)
+                        fg_data_glb(c,r) = &
+                             fg_data_glb(c,r) * &
+                          agrmet_struc(nest)%galwem_nrt_bias_ratio(c,r)
+                     end do
+                  end do
+               end if
+
             end if
          endif
 
@@ -851,9 +873,10 @@ subroutine USAF_getBackNWP(nest,back4,pcp_src, use_twelve, j6hr, findex)
             end if
 
             ! Apply bias correction to background field
-            if (agrmet_struc(nest)%back_bias_corr .eq. 1) then
+            if (agrmet_struc(nest)%back_bias_corr .eq. 1 .and. &
+                 rc .eq. 0) then
                write(LIS_logunit,*) &
-                   '[INFO] Applying GALWEM-based bias correction to GFS precip'
+             '[INFO] Applying GALWEM-based bias correction to GFS precip'
                do r = 1, LIS_rc%gnr(nest)
                   do c = 1, LIS_rc%gnc(nest)
                      fg_data_glb(c,r) = &
@@ -861,6 +884,19 @@ subroutine USAF_getBackNWP(nest,back4,pcp_src, use_twelve, j6hr, findex)
                           agrmet_struc(nest)%pcp_back_bias_ratio(c,r)
                   end do
                end do
+
+            else if (agrmet_struc(nest)%back_bias_corr .eq. 2 .and. &
+                 rc .eq. 0) then
+               ! New IMERG option for NRT
+               write(LIS_logunit,*) &
+             '[INFO] Applying IMERG-based bias correction to GFS precip'
+               do r = 1, LIS_rc%gnr(nest)
+                  do c = 1, LIS_rc%gnc(nest)
+                     fg_data_glb(c,r) = &
+                          fg_data_glb(c,r) * &
+                          agrmet_struc(nest)%gfs_nrt_bias_ratio(c,r)
+                  end do
+               end do
             end if
 
             ! Use GFS Bratseth settings if we have the data
@@ -946,6 +982,8 @@ subroutine USAF_getSSMIObsData(nest,j6hr,use_twelve,precip3,precip6, &
       integer :: local_diff_grid
       integer :: ii,jj,kk,count_good_ssmi ! EMK TEST
 
+      external :: agrmet_ssmiprec_filename
+
       TRACE_ENTER("bratseth_getSSMI")
       routine_name = 'USAF_getSSMIObsData'
 
@@ -1135,6 +1173,9 @@ subroutine USAF_getGeoPrecipObsData(nest,j6hr,use_twelve,precip3,precip6,&
       logical, external :: is_geo_corrupted
       real, allocatable :: gest_temp(:,:,:)
 
+      external :: agrmet_geoprec_filename
+      external :: polarToLatLon
+
       TRACE_ENTER("bratseth_getGeopPrcp")
       net = "GEOPRECIP"
       platform = "GEOPRECIP"
@@ -1546,7 +1587,6 @@ subroutine USAF_analyzePrecip(precipAll,nest,back,hourindex,mrgp,precipOBA)
       use AGRMET_forcingMod, only:  agrmet_struc
       use LIS_coreMod, only: LIS_rc, LIS_ews_halo_ind, LIS_ewe_halo_ind, &
            LIS_nss_halo_ind, LIS_nse_halo_ind, LIS_localPet
-      use LIS_logMod, only: LIS_logunit
       use USAF_OBAMod, only: OBA, createOBA,assignOBA
 
       ! Defaults
@@ -1565,8 +1605,6 @@ subroutine USAF_analyzePrecip(precipAll,nest,back,hourindex,mrgp,precipOBA)
       real, allocatable :: invDataDensities(:)
       real, allocatable :: sumObsEstimates(:)
       integer :: npasses
-      integer :: r
-      character(len=10) :: new_name,type
       real :: convergeThresh
       real :: sigmaBSqr
       integer :: good_obs
@@ -1688,6 +1726,9 @@ subroutine interpBack2ObsPts(this,nest,imax,jmax,back, &
       integer, external :: get_fieldpos
       real, parameter :: FILL = MISSING
 
+      external :: compute_grid_coord
+      external :: bilinear_interp
+
       ! See if we have any observations available
       nobs = this%nobs
       if (nobs .eq. 0) return
@@ -1831,6 +1872,8 @@ subroutine calc_invDataDensities(this,sigmaBSqr,nest,max_dist, &
       integer :: imax,jmax
       logical :: verbose
 
+      external :: MPI_Barrier, MPI_ALLREDUCE
+
       verbose = .true.
       if (present(silent)) then
          if (silent) verbose = .false.
@@ -2055,6 +2098,8 @@ subroutine calc_obsAnalysis(this,sigmaBSqr,nobs,invDataDensities,nest,&
       real :: O, A
       integer :: good_obs
 
+      external :: MPI_Barrier, MPI_ALLREDUCE
+
       verbose = .true.
       if (present(silent)) then
          if (silent) verbose = .false.
@@ -2506,6 +2551,8 @@ subroutine calc_gridAnalysis(this,nest,sigmaBSqr,nobs,invDataDensities,&
       integer :: r_local, c_local
       integer :: gindex
 
+      external :: MPI_Barrier, MPI_ALLREDUCE
+
       ! Sanity checks
       if (nobs .eq. 0) return
 
@@ -2765,6 +2812,11 @@ subroutine fldbld_precip_nwp(nest,findex,julhr,src,fc_hr, &
       integer :: rc2
       logical :: first_time
 
+      external :: getAVNfilename
+      external :: AGRMET_getGALWEMfilename
+      external :: AGRMET_fg2lis_precip
+      external :: interp_galwem_first_guess
+
       rc = 0
 
       found = .false.
@@ -2974,7 +3026,6 @@ subroutine check_grib_file(gribfile,yr1,mo1,da1,hr1,found, &
 #if (defined USE_GRIBAPI)
       use grib_api
 #endif
-      use LIS_coreMod, only: LIS_masterproc
       use LIS_logMod,  only : LIS_logunit, LIS_abort, LIS_alert, &
            LIS_verify, LIS_endrun
       use LIS_mpiMod
@@ -3203,6 +3254,8 @@ subroutine USAF_superstatQC(this,nest,new_name,network,silent_rejects)
       double precision :: t1, t2
       logical :: silent_rejects_local
 
+      external :: MPI_Barrier, MPI_Allreduce
+
       ! Sanity check
       nobs = this%nobs
       if (nobs .eq. 0) then
@@ -3617,6 +3670,8 @@ subroutine USAF_dupQC(this)
       integer :: ierr
       logical :: location_issue
 
+      external :: MPI_Barrier
+
       nobs = this%nobs
       if (nobs .eq. 0) then
          write(LIS_logunit,*)&
@@ -4006,6 +4061,8 @@ subroutine USAF_backQC(this,sigmaBSqr,silent_rejects)
       double precision :: t1, t2
       logical :: silent_rejects_local
 
+      external :: MPI_Barrier
+      
       nobs = this%nobs
       if (nobs .eq. 0) then
          write(LIS_logunit,*)&
@@ -4093,7 +4150,9 @@ end function is_gauge
    logical function is_stn(net)
       implicit none
       character(len=10), intent(in) :: net
+      character(len=10) :: net_local
       logical :: answer
+      net_local = net
       answer = .true.
       is_stn = answer
    end function is_stn
@@ -4367,8 +4426,7 @@ end subroutine destroy_hash2d
    subroutine build_hash2d(this,nest,imax,jmax,hash2d)
 
       ! Imports
-      use LIS_coreMod, only: LIS_rc, LIS_domain, LIS_localPet
-      use LIS_logMod, only: LIS_logunit
+      use LIS_coreMod, only: LIS_rc
       use LIS_mpiMod
 
       ! Defaults
@@ -4382,26 +4440,21 @@ subroutine build_hash2d(this,nest,imax,jmax,hash2d)
       type(hash_list), allocatable, intent(out) :: hash2d(:,:)
 
       ! Local variables
-      real :: dlat, dlon, ctrlat, ctrlon
       integer :: nobs
       integer, allocatable :: cols(:), rows(:)
-      integer :: pet, pet_incr
-      logical :: found
-      integer :: j,r,c,gindex
-      integer :: ierr
 
       nobs = this%nobs
       if (nobs .eq. 0) return
 
       ! Get the columns and rows of the observations in the LIS grid
-      call find_LIS_cols_rows(this,nest,cols,rows)
+      call find_LIS_cols_rows(this, nest, cols, rows)
 
       ! Loop through the collected col, row information and add to the
       ! hash2d table.
-      call init_hash2d(LIS_rc%gnc(nest),LIS_rc%gnr(nest),hash2d)
+      call init_hash2d(LIS_rc%gnc(nest), LIS_rc%gnr(nest), hash2d)
       imax = LIS_rc%gnc(nest) ! Output
       jmax = LIS_rc%gnr(nest) ! Output
-      call insert_hash2d_array(imax,jmax,hash2d,nobs,cols,rows)
+      call insert_hash2d_array(imax, jmax, hash2d, nobs, cols, rows)
 !       do j = 1,nobs
 !          c = cols(j)
 !          if (c .eq. 0) cycle
@@ -4786,6 +4839,8 @@ subroutine find_LIS_cols_rows(this,nest,cols,rows)
       integer :: j,r,c,gindex
       integer :: ierr
 
+      external :: MPI_Barrier, MPI_ALLREDUCE
+
       nobs = this%nobs
       if (nobs .eq. 0) return
 
@@ -4910,12 +4965,10 @@ end subroutine find_LIS_cols_rows
    subroutine USAF_waterQC(this,nest,silent_rejects)
 
       ! Imports
-      use LIS_coreMod, only: LIS_rc
       use LIS_LMLCMod, only: LIS_LMLC
       use LIS_logMod, only: LIS_logunit
       use LIS_mpiMod
 
-
       ! Defaults
       implicit none
 
@@ -4934,6 +4987,8 @@ subroutine USAF_waterQC(this,nest,silent_rejects)
       integer :: reject_count
       logical :: silent_rejects_local
 
+      external :: MPI_Barrier
+
       ! Sanity check
       nobs = this%nobs
       if (nobs .eq. 0) then
@@ -5036,6 +5091,8 @@ subroutine USAF_snowQC(this,nest,hourindex,threshold,silent_rejects)
       real :: threshold_local
       logical :: silent_rejects_local
 
+      external :: MPI_Barrier, MPI_ALLREDUCE
+
       ! Sanity check
       nobs = this%nobs
       if (nobs .eq. 0) then
@@ -5219,6 +5276,8 @@ subroutine USAF_snowDepthQC(this,nest,silent_rejects)
       integer :: rglb,cglb
       integer :: gid
 
+      external :: MPI_Barrier, MPI_ALLREDUCE
+
       ! Sanity check
       nobs = this%nobs
       if (nobs .eq. 0) then
@@ -5440,9 +5499,9 @@ subroutine USAF_analyzeScreen(screenObs,nest,back,sigmaBSqr, &
       real, allocatable :: invDataDensities(:)
       real, allocatable :: sumObsEstimates(:)
       integer :: npasses
-      character(len=10) :: new_name,type
+      character(len=10) :: new_name, type
       real :: convergeThresh
-      integer :: i,j
+      integer :: j
 
       TRACE_ENTER("bratseth_analyzeScrn")
       ! Initialize analysis field with the background first guess.  This will
@@ -5614,6 +5673,8 @@ subroutine USAF_getCMORPHObsData(nest,j6hr,use_twelve, &
       integer :: ftn, ios
       integer :: i,j,k
 
+      external :: cmorfile_agrmet
+
       TRACE_ENTER("bratseth_getCMORPH")
       net = "CMORPH"
       platform = "CMORPH"
@@ -6397,14 +6458,15 @@ subroutine USAF_setBratsethScreenStats(src,n)
 
    end subroutine USAF_setBratsethScreenStats
 
-   ! Read bias ratio for background field
+   ! Read bias ratio for background field (S2S version, which uses
+   ! data in the LDT parameter file)
    subroutine USAF_pcpBackBiasRatio_s2s(n, yyyymmddhh)
 
      ! Imports
      use AGRMET_forcingMod, only: agrmet_struc
-     use LIS_coreMod, only: LIS_rc, LIS_localPet, &
-          LIS_ews_halo_ind, LIS_ewe_halo_ind, &
-          LIS_nss_halo_ind, LIS_nse_halo_ind
+     use LIS_coreMod, only: LIS_rc!, LIS_localPet, &
+          !LIS_ews_halo_ind, LIS_ewe_halo_ind, &
+          !LIS_nss_halo_ind, LIS_nse_halo_ind
      use LIS_logMod, only: LIS_logunit, LIS_endrun, LIS_verify
 #if ( defined USE_NETCDF3 || defined USE_NETCDF4 )
      use netcdf
@@ -6468,4 +6530,101 @@ subroutine USAF_pcpBackBiasRatio_s2s(n, yyyymmddhh)
 #endif
 
    end subroutine USAF_pcpBackBiasRatio_s2s
-end module USAF_bratsethMod
+
+   ! Read bias ratio for background field (NRT version, which uses
+   ! data in standalone netCDF file)
+   subroutine USAF_pcpBackBiasRatio_nrt(n, yyyymmddhh)
+
+     ! Imports
+     use AGRMET_forcingMod, only: agrmet_struc
+     use LIS_coreMod, only: LIS_rc!, LIS_localPet, &
+          !LIS_ews_halo_ind, LIS_ewe_halo_ind, &
+          !LIS_nss_halo_ind, LIS_nse_halo_ind
+
+     use LIS_logMod, only: LIS_logunit, LIS_endrun, LIS_verify
+#if ( defined USE_NETCDF3 || defined USE_NETCDF4 )
+     use netcdf
+#endif
+
+     ! Defaults
+     implicit none
+
+     ! Arguments
+     integer, intent(in) :: n
+     character(len=10), intent(in) :: yyyymmddhh
+
+     ! Locals
+     integer :: imonth
+     logical :: file_exists
+     integer :: ncid, ppt_ratio_Id
+
+#if ( defined USE_NETCDF3 || defined USE_NETCDF4 )
+
+     read(yyyymmddhh(5:6),'(i2.2)') imonth
+
+     ! Accumulations ending at 00Z first of month are part of the
+     ! *previous* month.  So decrement imonth by one in that situation,
+     ! and reset to 12 (December) at the year change.
+     if (agrmet_struc(n)%pcp_back_bias_ratio_month .ne. 0) then
+        if (yyyymmddhh(9:10) .eq. '00' .and. &
+             yyyymmddhh(7:8) .eq. '01') then
+           imonth = imonth - 1
+           if (imonth .eq. 0) then
+              imonth = 12
+           end if
+        end if
+     end if
+     if (imonth .eq. agrmet_struc(n)%pcp_back_bias_ratio_month .and. &
+          agrmet_struc(n)%pcp_back_bias_ratio_month .ne. 0) return
+
+     agrmet_struc(n)%pcp_back_bias_ratio_month = imonth
+
+     ! First, update the GFS ratios
+     inquire(file=agrmet_struc(n)%gfs_nrt_bias_ratio_file, &
+          exist=file_exists)
+     if (.not. file_exists) then
+        write(LIS_logunit,*) '[ERR] Cannot find ', &
+             trim(agrmet_struc(n)%gfs_nrt_bias_ratio_file)
+        call LIS_endrun()
+     end if
+     call LIS_verify( &
+          nf90_open(path=agrmet_struc(n)%gfs_nrt_bias_ratio_file, &
+          mode=NF90_NOWRITE, ncid=ncid), &
+          'Error in nf90_open in USAF_pcpBackBiasRatio_nrt')
+     call LIS_verify(nf90_inq_varid(ncid, "biasRatio", ppt_ratio_Id), &
+          'biasRatio field not found in bias file')
+     call LIS_verify(nf90_get_var(ncid, ppt_ratio_Id, &
+          agrmet_struc(n)%gfs_nrt_bias_ratio, &
+          start=(/1,1,imonth/), &
+          count=(/LIS_rc%gnc(n),LIS_rc%gnr(n),1/)), &
+          'n90_get_var failed for biasRatio in LIS param file')
+     call LIS_verify(nf90_close(ncid),&
+          'nf90_close failed in USAF_pcpBackBiasRatio_nrt')
+
+     ! Repeat for GALWEM ratios
+     inquire(file=agrmet_struc(n)%galwem_nrt_bias_ratio_file, &
+          exist=file_exists)
+     if (.not. file_exists) then
+        write(LIS_logunit,*) '[ERR] Cannot find ', &
+             trim(agrmet_struc(n)%galwem_nrt_bias_ratio_file)
+        call LIS_endrun()
+     end if
+     call LIS_verify( &
+          nf90_open(path=agrmet_struc(n)%galwem_nrt_bias_ratio_file, &
+          mode=NF90_NOWRITE, ncid=ncid), &
+          'Error in nf90_open in USAF_pcpBackBiasRatio_nrt')
+     call LIS_verify(nf90_inq_varid(ncid, "biasRatio", ppt_ratio_Id), &
+          'biasRatio field not found in bias file')
+     call LIS_verify(nf90_get_var(ncid, ppt_ratio_Id, &
+          agrmet_struc(n)%galwem_nrt_bias_ratio, &
+          start=(/1,1,imonth/), &
+          count=(/LIS_rc%gnc(n),LIS_rc%gnr(n),1/)),&
+          'n90_get_var failed for biasRatio in LIS param file')
+     call LIS_verify(nf90_close(ncid),&
+          'nf90_close failed in USAF_pcpBackBiasRatio_nrt')
+
+#endif
+
+   end subroutine USAF_pcpBackBiasRatio_nrt
+
+ end module USAF_bratsethMod
diff --git a/lis/metforcing/usaf/readcrd_agrmet.F90 b/lis/metforcing/usaf/readcrd_agrmet.F90
index 05452a3a4..819905484 100644
--- a/lis/metforcing/usaf/readcrd_agrmet.F90
+++ b/lis/metforcing/usaf/readcrd_agrmet.F90
@@ -62,6 +62,10 @@ subroutine readcrd_agrmet()
   integer, external :: LIS_create_subdirs
   integer :: tmp_imerg_plp_thresh
   integer :: ierr
+  logical :: use_nrt_bias_files ! EMK
+
+  external :: MPI_Barrier
+  external :: sleep
 
   call ESMF_ConfigFindLabel(LIS_config,"AGRMET forcing directory:",rc=rc)
   do n=1,LIS_rc%nnest
@@ -1091,16 +1095,59 @@ subroutine readcrd_agrmet()
      call LIS_verify(rc, &
           "[ERR] AGRMET PPT Background bias correction option: not specified in config file")
      if (agrmet_struc(n)%back_bias_corr .lt. 0 .or. &
-          agrmet_struc(n)%back_bias_corr .gt. 1) then
+          agrmet_struc(n)%back_bias_corr .gt. 2) then
         call LIS_verify(rc, &
-             "[ERR] AGRMET PPT Background bias correction option: bad value in config file, set 0 or 1")
+             "[ERR] AGRMET PPT Background bias correction option: bad value in config file, set 0, 1, or 2")
      end if
      if (agrmet_struc(n)%back_bias_corr .eq. 1) then
         allocate(agrmet_struc(n)%pcp_back_bias_ratio(LIS_rc%gnc(n),LIS_rc%gnr(n)))
         agrmet_struc(n)%pcp_back_bias_ratio = 1.
         agrmet_struc(n)%pcp_back_bias_ratio_month = 0
+     else if (agrmet_struc(n)%back_bias_corr .eq. 2) then
+
+        allocate(agrmet_struc(n)%gfs_nrt_bias_ratio( &
+             LIS_rc%gnc(n),LIS_rc%gnr(n)))
+        allocate(agrmet_struc(n)%galwem_nrt_bias_ratio( &
+             LIS_rc%gnc(n),LIS_rc%gnr(n)))
+        agrmet_struc(n)%gfs_nrt_bias_ratio = 1.
+        agrmet_struc(n)%galwem_nrt_bias_ratio = 1.
+        agrmet_struc(n)%pcp_back_bias_ratio_month = 0
      end if
-  enddo ! n
+  end do
+
+  ! EMK Add support for NRT bias files for GFS and GALWEM
+  use_nrt_bias_files = .false.
+  do n = 1, LIS_rc%nnest
+     if (agrmet_struc(n)%back_bias_corr .eq. 2) then
+        use_nrt_bias_files = .true.
+        exit
+     end if
+  end do
+
+  if (use_nrt_bias_files) then
+     call ESMF_ConfigFindLabel(LIS_config, &
+          "AGRMET PPT GFS NRT bias file:", rc=rc)
+     call LIS_verify(rc, &
+          "[ERR] AGRMET PPT GFS NRT bias file: not specified in config file")
+
+     do n = 1, LIS_rc%nnest
+        call ESMF_ConfigGetAttribute(LIS_config, &
+             agrmet_struc(n)%gfs_nrt_bias_ratio_file, rc=rc)
+        call LIS_verify(rc, &
+             "[ERR] AGRMET PPT GFS NRT bias file: not specified in config file")
+     end do
+
+     call ESMF_ConfigFindLabel(LIS_config, &
+          "AGRMET PPT GALWEM NRT bias file:", rc=rc)
+     call LIS_verify(rc, &
+          "[ERR] AGRMET PPT GALWEM NRT bias file: not specified in config file")
+     do n = 1, LIS_rc%nnest
+        call ESMF_ConfigGetAttribute(LIS_config, &
+             agrmet_struc(n)%galwem_nrt_bias_ratio_file, rc=rc)
+        call LIS_verify(rc, &
+             "[ERR] AGRMET PPT GALWEM NRT bias file: not specified in config file")
+     end do
+  end if
 
   do n=1,LIS_rc%nnest
      agrmet_struc(n)%radProcessInterval = 1

From 6e4b714ae36aa5c7850b1b531d9c36d4bf75888a Mon Sep 17 00:00:00 2001
From: Eric Kemp <eric.kemp@nasa.gov>
Date: Wed, 20 Sep 2023 10:39:31 -0400
Subject: [PATCH 03/16] Attempt at remapping IMERG-FR monthly data.

This version uses the precipitationQualityIndex field to screen out
apparent data voids (i.e., areas w/o gauge corrections).

HOWEVER, a better method may be to use the gaugeRelativeWeighting field,
which indicates the Kalman weight given to the GPCC gauge data in the
IMERG-Final Run.  More work and testing is anticipated.
---
 .../IMERG_monthly/IMERG_monthly_dataMod.F90   |  2 +-
 .../IMERG_monthly/readIMERGmonthlydata.F90    | 91 ++++++++++++++++---
 lvt/plugins/LVT_datastream_pluginMod.F90      |  2 +-
 3 files changed, 80 insertions(+), 15 deletions(-)

diff --git a/lvt/datastreams/IMERG_monthly/IMERG_monthly_dataMod.F90 b/lvt/datastreams/IMERG_monthly/IMERG_monthly_dataMod.F90
index 5321bdb3e..9a98a1c91 100644
--- a/lvt/datastreams/IMERG_monthly/IMERG_monthly_dataMod.F90
+++ b/lvt/datastreams/IMERG_monthly/IMERG_monthly_dataMod.F90
@@ -80,7 +80,7 @@ subroutine IMERG_monthly_datainit(i)
     call LVT_verify(status, 'IMERG monthly data directory: not defined')
 
     ! Get IMERG product
-    call ESMF_ConfigFindLabel(LVT_config," IMERG monthly product:", &
+    call ESMF_ConfigFindLabel(LVT_config,"IMERG monthly product:", &
          rc=status)
     call LVT_verify(status, 'IMERG monthly product: not defined')
     call ESMF_ConfigGetAttribute(LVT_config, &
diff --git a/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90 b/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90
index 1ee98dd7b..ddadcb844 100644
--- a/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90
+++ b/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90
@@ -59,8 +59,8 @@ subroutine readIMERGmonthlydata(source)
   external :: upscaleByAveraging
 
   ! Initialize variables
-  prcp(:) = LVT_rc%udef
-  prcp_final(:,:) = LVT_rc%udef
+  prcp = LVT_rc%udef
+  prcp_final = LVT_rc%udef
   currTime = float(LVT_rc%dhr(source))*3600 + &
        60*LVT_rc%dmn(source) + LVT_rc%dss(source)
 
@@ -77,8 +77,8 @@ subroutine readIMERGmonthlydata(source)
   hr1 = LVT_rc%dhr(source)
   mn1 = LVT_rc%dmn(source)
   ss1 = 0
-  call ESMF_TimeSet(time1,yy=yr1, mm=mo1, dd=da1, &
-       h=hr1,m=mn1,s=ss1,calendar=LVT_calendar,rc=status)
+  call ESMF_TimeSet(time1, yy=yr1, mm=mo1, dd=da1, &
+       h=hr1, m=mn1, s=ss1, calendar=LVT_calendar, rc=status)
   call LVT_verify(status)
 
   ! Use previous month file.
@@ -92,7 +92,7 @@ subroutine readIMERGmonthlydata(source)
   call LVT_verify(status)
 
   if (alarmCheck) then
-     call create_IMERG_monthly_filename(source, yr1, mo1, filename)
+     call create_IMERG_monthly_filename(source, yr2, mo2, filename)
      inquire(file=trim(filename), exist=file_exists)
 
      if (file_exists) then
@@ -105,8 +105,8 @@ subroutine readIMERGmonthlydata(source)
            ! Use budget-bilinear interpolation if IMERG data are at
            ! coarser resolution than the analysis grid; otherwise, use
            ! upscale averaging.
-           prcp_in1(:) = 0
-           lb(:) = .false.
+           prcp_in1 = LVT_rc%udef
+           lb = .false.
            t = 1
            do r = 1, imergmonthlydata(source)%nr
               do c = 1,imergmonthlydata(source)%nc
@@ -118,6 +118,12 @@ subroutine readIMERGmonthlydata(source)
               end do ! c
            end do ! r
 
+           write(LVT_logunit,*)'EMK: maxval(prcp_in1) = ', &
+                maxval(prcp_in1)
+           write(LVT_logunit,*)'EMK: minval(prcp_in1) = ', &
+                minval(prcp_in1)
+           flush(LVT_logunit)
+
            if (LVT_isAtAFinerResolution( &
                 imergmonthlydata(source)%datares)) then
               call conserv_interp(LVT_rc%gridDesc, lb, prcp_in1, &
@@ -181,10 +187,9 @@ subroutine readIMERGmonthlydata(source)
   lis_ts = time1 - time2
   call ESMF_TimeIntervalGet(lis_ts, d=days_in_month, rc=status)
   call LVT_verify(status)
-  write(LVT_logunit,*)'EMK: days_in_month = ', days_in_month
   do r = 1, LVT_rc%lnr
      do c = 1, LVT_rc%lnc
-        if (prcp_final(c,r) .ge. 0) then
+        if (prcp_final(c,r) >= 0) then
            prcp_final(c,r) = &
                 prcp_final(c,r) * days_in_month * 86400 ! kg/m2
         else
@@ -195,11 +200,17 @@ subroutine readIMERGmonthlydata(source)
   call LVT_logSingleDataStreamVar(LVT_MOC_totalprecip, source, &
        prcp_final, vlevel=1, units='kg/m2')
 
+  write(LVT_logunit,*)'EMK: LVT_rc%udef = ', LVT_rc%udef
+  write(LVT_logunit,*)'EMK: maxval(precip) = ', maxval(prcp_final)
+  write(LVT_logunit,*)'EMK: minval(precip) = ', minval(prcp_final)
+  flush(LVT_logunit)
+
 end subroutine readIMERGmonthlydata
 
 subroutine read_imergmonthly_hdf(filename, col, row, precipout, ireaderr)
 
   ! Imports
+  use LVT_coreMod, only: LVT_rc
   use LVT_logMod, only: LVT_logunit
 
 #if (defined USE_HDF5)
@@ -217,14 +228,21 @@ subroutine read_imergmonthly_hdf(filename, col, row, precipout, ireaderr)
 
   ! Local variables
   integer :: xsize, ysize
-  character(len=40) :: dsetname='/Grid/precipitation'
+  character(len=40) :: dsetname = '/Grid/precipitation'
+  character(len=40) :: dsetname_qi = '/Grid/precipitationQualityIndex'
+  real :: qiin(row,col)
   real :: precipin(row,col)
   integer :: istatus
+  integer :: i, j
 #if (defined USE_HDF5)
   integer(HSIZE_T), dimension(2) :: dims
   integer(HID_T) :: fileid, dsetid
 #endif
 
+  qiin = LVT_rc%udef
+  precipin = LVT_rc%udef
+  precipout = LVT_rc%udef
+
 #if (defined USE_HDF5)
   xsize = col
   ysize = row
@@ -250,7 +268,42 @@ subroutine read_imergmonthly_hdf(filename, col, row, precipout, ireaderr)
      return
   end if
 
-  ! Open dataset
+  ! Open quality index dataset
+  call h5dopen_f(fileid, dsetname_qi, dsetid, istatus)
+  if (istatus .ne. 0) then
+     write(LVT_logunit,*) 'Error opening IMERG dataset', &
+          trim(dsetname_qi)
+     ireaderr = istatus
+     call h5fclose_f(fileid, istatus) ! Close HDF5 file
+     call h5close_f(istatus) ! Close HDF5 interface
+     return
+  end if
+
+  ! Read quality index dataset
+  call h5dread_f(dsetid, H5T_NATIVE_REAL, &
+       qiin, dims, istatus)
+  if (istatus .ne. 0) then
+     write(LVT_logunit,*) 'Error reading IMERG dataset', &
+          trim(dsetname_qi)
+     ireaderr = istatus
+     call h5dclose_f(dsetid, istatus) ! Close dataset
+     call h5fclose_f(fileid, istatus) ! Close HDF5 file
+     call h5close_f(istatus) ! Close HDF5 interface
+     return
+  end if
+
+  ! Close quality index dataset
+  call h5dclose_f(dsetid, istatus)
+  if (istatus .ne. 0) then
+     write(LVT_logunit,*) 'Error closing IMERG dataset', &
+          trim(dsetname_qi)
+     ireaderr = istatus
+     call h5fclose_f(fileid, istatus) ! Close HDF5 file
+     call h5close_f(istatus) ! Close HDF5 interface
+     return
+  end if
+
+  ! Open precip dataset
   call h5dopen_f(fileid, dsetname, dsetid, istatus)
   if (istatus .ne. 0) then
      write(LVT_logunit,*) 'Error opening IMERG dataset', trim(dsetname)
@@ -271,6 +324,17 @@ subroutine read_imergmonthly_hdf(filename, col, row, precipout, ireaderr)
      return
   end if
 
+  ! Flag precipin values where quality index is less then 2.5 (interpreted
+  ! as less than 2.5 gauges used in analysis).  Note that polar regions
+  ! have values ~2.1.
+  do j = 1, col
+    do i = 1, row
+       if (qiin(i,j) < 2.5) then
+          precipin(i,j) = LVT_rc%udef
+       end if
+    end do
+  end do
+
   ! Put the real(1:,1:) on the precipout(0:,0:)
   ! precipin is (ysize,xsize) starting at (lon=-179.9,lat=-89.9)
   precipout(1:xsize,1:ysize) = transpose(precipin)
@@ -321,7 +385,8 @@ subroutine create_IMERG_monthly_filename(source, yr, mo, filename)
      call LVT_endrun()
   endif
 
-  filename = trim(odir) // "/" // cyr // cmo // trim(fstem) // &
-       cyr // cmo // "01-S000000-E235959." // cmo // trim(imVer) // fext
+  filename = trim(odir) // "/" // cyr // trim(fstem) // &
+       cyr // cmo // "01-S000000-E235959." // cmo // "." // &
+       trim(imVer) // fext
 
 end subroutine create_IMERG_monthly_filename
diff --git a/lvt/plugins/LVT_datastream_pluginMod.F90 b/lvt/plugins/LVT_datastream_pluginMod.F90
index 9259286b7..8a5302abc 100644
--- a/lvt/plugins/LVT_datastream_pluginMod.F90
+++ b/lvt/plugins/LVT_datastream_pluginMod.F90
@@ -706,7 +706,7 @@ subroutine LVT_datastream_plugin
 
     call registerobssetup(trim(LVT_IMERGmonthlydataId)//char(0), &
          IMERG_monthly_datainit)
-    call registerobsread(trim(LVT_IMERGdataId)//char(0), &
+    call registerobsread(trim(LVT_IMERGmonthlydataId)//char(0), &
          readIMERGmonthlydata)
 
     call registerobssetup(trim(LVT_UASNOWdataId)//char(0), UASNOW_obsinit)

From 95afa5276ea78b6da035486af25340e06190a0a0 Mon Sep 17 00:00:00 2001
From: Eric Kemp <eric.kemp@nasa.gov>
Date: Wed, 4 Oct 2023 15:23:16 -0400
Subject: [PATCH 04/16] Alternative QC check.

---
 .../IMERG_monthly/readIMERGmonthlydata.F90    | 93 ++++++++++++++-----
 1 file changed, 72 insertions(+), 21 deletions(-)

diff --git a/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90 b/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90
index ddadcb844..596b275d1 100644
--- a/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90
+++ b/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90
@@ -229,8 +229,10 @@ subroutine read_imergmonthly_hdf(filename, col, row, precipout, ireaderr)
   ! Local variables
   integer :: xsize, ysize
   character(len=40) :: dsetname = '/Grid/precipitation'
-  character(len=40) :: dsetname_qi = '/Grid/precipitationQualityIndex'
-  real :: qiin(row,col)
+  !character(len=40) :: dsetname_qi = '/Grid/precipitationQualityIndex'
+  character(len=40) :: dsetname_wgt = '/Grid/gaugeRelativeWeighting'
+  !real :: qiin(row,col)
+  real :: wgtin(row,col)
   real :: precipin(row,col)
   integer :: istatus
   integer :: i, j
@@ -239,7 +241,8 @@ subroutine read_imergmonthly_hdf(filename, col, row, precipout, ireaderr)
   integer(HID_T) :: fileid, dsetid
 #endif
 
-  qiin = LVT_rc%udef
+  !qiin = LVT_rc%udef
+  wgtin = LVT_rc%udef
   precipin = LVT_rc%udef
   precipout = LVT_rc%udef
 
@@ -268,23 +271,58 @@ subroutine read_imergmonthly_hdf(filename, col, row, precipout, ireaderr)
      return
   end if
 
-  ! Open quality index dataset
-  call h5dopen_f(fileid, dsetname_qi, dsetid, istatus)
+  ! ! Open quality index dataset
+  ! call h5dopen_f(fileid, dsetname_qi, dsetid, istatus)
+  ! if (istatus .ne. 0) then
+  !    write(LVT_logunit,*) 'Error opening IMERG dataset', &
+  !         trim(dsetname_qi)
+  !    ireaderr = istatus
+  !    call h5fclose_f(fileid, istatus) ! Close HDF5 file
+  !    call h5close_f(istatus) ! Close HDF5 interface
+  !    return
+  ! end if
+
+  ! ! Read quality index dataset
+  ! call h5dread_f(dsetid, H5T_NATIVE_REAL, &
+  !      qiin, dims, istatus)
+  ! if (istatus .ne. 0) then
+  !    write(LVT_logunit,*) 'Error reading IMERG dataset', &
+  !         trim(dsetname_qi)
+  !    ireaderr = istatus
+  !    call h5dclose_f(dsetid, istatus) ! Close dataset
+  !    call h5fclose_f(fileid, istatus) ! Close HDF5 file
+  !    call h5close_f(istatus) ! Close HDF5 interface
+  !    return
+  ! end if
+
+  ! ! Close quality index dataset
+  ! call h5dclose_f(dsetid, istatus)
+  ! if (istatus .ne. 0) then
+  !    write(LVT_logunit,*) 'Error closing IMERG dataset', &
+  !         trim(dsetname_qi)
+  !    ireaderr = istatus
+  !    call h5fclose_f(fileid, istatus) ! Close HDF5 file
+  !    call h5close_f(istatus) ! Close HDF5 interface
+  !    return
+  ! end if
+
+  ! Open gauge relative weight dataset
+  call h5dopen_f(fileid, dsetname_wgt, dsetid, istatus)
   if (istatus .ne. 0) then
      write(LVT_logunit,*) 'Error opening IMERG dataset', &
-          trim(dsetname_qi)
+          trim(dsetname_wgt)
      ireaderr = istatus
      call h5fclose_f(fileid, istatus) ! Close HDF5 file
      call h5close_f(istatus) ! Close HDF5 interface
      return
   end if
 
-  ! Read quality index dataset
+  ! Read gauge relative weight dataset
   call h5dread_f(dsetid, H5T_NATIVE_REAL, &
-       qiin, dims, istatus)
+       wgtin, dims, istatus)
   if (istatus .ne. 0) then
      write(LVT_logunit,*) 'Error reading IMERG dataset', &
-          trim(dsetname_qi)
+          trim(dsetname_wgt)
      ireaderr = istatus
      call h5dclose_f(dsetid, istatus) ! Close dataset
      call h5fclose_f(fileid, istatus) ! Close HDF5 file
@@ -292,11 +330,11 @@ subroutine read_imergmonthly_hdf(filename, col, row, precipout, ireaderr)
      return
   end if
 
-  ! Close quality index dataset
+  ! Close gauge relative weight dataset
   call h5dclose_f(dsetid, istatus)
   if (istatus .ne. 0) then
      write(LVT_logunit,*) 'Error closing IMERG dataset', &
-          trim(dsetname_qi)
+          trim(dsetname_wgt)
      ireaderr = istatus
      call h5fclose_f(fileid, istatus) ! Close HDF5 file
      call h5close_f(istatus) ! Close HDF5 interface
@@ -324,16 +362,29 @@ subroutine read_imergmonthly_hdf(filename, col, row, precipout, ireaderr)
      return
   end if
 
-  ! Flag precipin values where quality index is less then 2.5 (interpreted
-  ! as less than 2.5 gauges used in analysis).  Note that polar regions
-  ! have values ~2.1.
-  do j = 1, col
-    do i = 1, row
-       if (qiin(i,j) < 2.5) then
-          precipin(i,j) = LVT_rc%udef
-       end if
-    end do
-  end do
+  ! ! Flag precipin values where quality index is less then 2.5 (interpreted
+  ! ! as less than 2.5 gauges used in analysis).  Note that polar regions
+  ! ! have values ~2.1.
+  ! do j = 1, col
+  !   do i = 1, row
+  !      if (qiin(i,j) < 2.5) then
+  !         precipin(i,j) = LVT_rc%udef
+  !      end if
+  !   end do
+  ! end do
+
+  ! Flag precipin values where gauge relative weighting is less then
+  ! 0.5 (meaning raw IMERG was weighted at least 50% of final product)
+  ! do j = 1, col
+  !    do i = 1, row
+  !      if (wgtin(i,j) < 0.50) then
+  !      !if (wgtin(i,j) < 0.25) then
+  !      !if (wgtin(i,j) < 0.1) then
+  !      !if (wgtin(i,j) < 0.01) then
+  !         precipin(i,j) = LVT_rc%udef
+  !      end if
+  !   end do
+  ! end do
 
   ! Put the real(1:,1:) on the precipout(0:,0:)
   ! precipin is (ysize,xsize) starting at (lon=-179.9,lat=-89.9)

From 96d78c26001d0f2124e3f6d933e84c1c7d0859b4 Mon Sep 17 00:00:00 2001
From: Eric Kemp <eric.kemp@nasa.gov>
Date: Wed, 4 Oct 2023 15:33:21 -0400
Subject: [PATCH 05/16] Removed QC of IMERG data.

All monthly data will be interpolated as-is.
---
 .../IMERG_monthly/readIMERGmonthlydata.F90    | 111 ------------------
 1 file changed, 111 deletions(-)

diff --git a/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90 b/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90
index 596b275d1..b0da4236e 100644
--- a/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90
+++ b/lvt/datastreams/IMERG_monthly/readIMERGmonthlydata.F90
@@ -118,12 +118,6 @@ subroutine readIMERGmonthlydata(source)
               end do ! c
            end do ! r
 
-           write(LVT_logunit,*)'EMK: maxval(prcp_in1) = ', &
-                maxval(prcp_in1)
-           write(LVT_logunit,*)'EMK: minval(prcp_in1) = ', &
-                minval(prcp_in1)
-           flush(LVT_logunit)
-
            if (LVT_isAtAFinerResolution( &
                 imergmonthlydata(source)%datares)) then
               call conserv_interp(LVT_rc%gridDesc, lb, prcp_in1, &
@@ -200,11 +194,6 @@ subroutine readIMERGmonthlydata(source)
   call LVT_logSingleDataStreamVar(LVT_MOC_totalprecip, source, &
        prcp_final, vlevel=1, units='kg/m2')
 
-  write(LVT_logunit,*)'EMK: LVT_rc%udef = ', LVT_rc%udef
-  write(LVT_logunit,*)'EMK: maxval(precip) = ', maxval(prcp_final)
-  write(LVT_logunit,*)'EMK: minval(precip) = ', minval(prcp_final)
-  flush(LVT_logunit)
-
 end subroutine readIMERGmonthlydata
 
 subroutine read_imergmonthly_hdf(filename, col, row, precipout, ireaderr)
@@ -229,10 +218,6 @@ subroutine read_imergmonthly_hdf(filename, col, row, precipout, ireaderr)
   ! Local variables
   integer :: xsize, ysize
   character(len=40) :: dsetname = '/Grid/precipitation'
-  !character(len=40) :: dsetname_qi = '/Grid/precipitationQualityIndex'
-  character(len=40) :: dsetname_wgt = '/Grid/gaugeRelativeWeighting'
-  !real :: qiin(row,col)
-  real :: wgtin(row,col)
   real :: precipin(row,col)
   integer :: istatus
   integer :: i, j
@@ -241,8 +226,6 @@ subroutine read_imergmonthly_hdf(filename, col, row, precipout, ireaderr)
   integer(HID_T) :: fileid, dsetid
 #endif
 
-  !qiin = LVT_rc%udef
-  wgtin = LVT_rc%udef
   precipin = LVT_rc%udef
   precipout = LVT_rc%udef
 
@@ -271,76 +254,6 @@ subroutine read_imergmonthly_hdf(filename, col, row, precipout, ireaderr)
      return
   end if
 
-  ! ! Open quality index dataset
-  ! call h5dopen_f(fileid, dsetname_qi, dsetid, istatus)
-  ! if (istatus .ne. 0) then
-  !    write(LVT_logunit,*) 'Error opening IMERG dataset', &
-  !         trim(dsetname_qi)
-  !    ireaderr = istatus
-  !    call h5fclose_f(fileid, istatus) ! Close HDF5 file
-  !    call h5close_f(istatus) ! Close HDF5 interface
-  !    return
-  ! end if
-
-  ! ! Read quality index dataset
-  ! call h5dread_f(dsetid, H5T_NATIVE_REAL, &
-  !      qiin, dims, istatus)
-  ! if (istatus .ne. 0) then
-  !    write(LVT_logunit,*) 'Error reading IMERG dataset', &
-  !         trim(dsetname_qi)
-  !    ireaderr = istatus
-  !    call h5dclose_f(dsetid, istatus) ! Close dataset
-  !    call h5fclose_f(fileid, istatus) ! Close HDF5 file
-  !    call h5close_f(istatus) ! Close HDF5 interface
-  !    return
-  ! end if
-
-  ! ! Close quality index dataset
-  ! call h5dclose_f(dsetid, istatus)
-  ! if (istatus .ne. 0) then
-  !    write(LVT_logunit,*) 'Error closing IMERG dataset', &
-  !         trim(dsetname_qi)
-  !    ireaderr = istatus
-  !    call h5fclose_f(fileid, istatus) ! Close HDF5 file
-  !    call h5close_f(istatus) ! Close HDF5 interface
-  !    return
-  ! end if
-
-  ! Open gauge relative weight dataset
-  call h5dopen_f(fileid, dsetname_wgt, dsetid, istatus)
-  if (istatus .ne. 0) then
-     write(LVT_logunit,*) 'Error opening IMERG dataset', &
-          trim(dsetname_wgt)
-     ireaderr = istatus
-     call h5fclose_f(fileid, istatus) ! Close HDF5 file
-     call h5close_f(istatus) ! Close HDF5 interface
-     return
-  end if
-
-  ! Read gauge relative weight dataset
-  call h5dread_f(dsetid, H5T_NATIVE_REAL, &
-       wgtin, dims, istatus)
-  if (istatus .ne. 0) then
-     write(LVT_logunit,*) 'Error reading IMERG dataset', &
-          trim(dsetname_wgt)
-     ireaderr = istatus
-     call h5dclose_f(dsetid, istatus) ! Close dataset
-     call h5fclose_f(fileid, istatus) ! Close HDF5 file
-     call h5close_f(istatus) ! Close HDF5 interface
-     return
-  end if
-
-  ! Close gauge relative weight dataset
-  call h5dclose_f(dsetid, istatus)
-  if (istatus .ne. 0) then
-     write(LVT_logunit,*) 'Error closing IMERG dataset', &
-          trim(dsetname_wgt)
-     ireaderr = istatus
-     call h5fclose_f(fileid, istatus) ! Close HDF5 file
-     call h5close_f(istatus) ! Close HDF5 interface
-     return
-  end if
-
   ! Open precip dataset
   call h5dopen_f(fileid, dsetname, dsetid, istatus)
   if (istatus .ne. 0) then
@@ -362,30 +275,6 @@ subroutine read_imergmonthly_hdf(filename, col, row, precipout, ireaderr)
      return
   end if
 
-  ! ! Flag precipin values where quality index is less then 2.5 (interpreted
-  ! ! as less than 2.5 gauges used in analysis).  Note that polar regions
-  ! ! have values ~2.1.
-  ! do j = 1, col
-  !   do i = 1, row
-  !      if (qiin(i,j) < 2.5) then
-  !         precipin(i,j) = LVT_rc%udef
-  !      end if
-  !   end do
-  ! end do
-
-  ! Flag precipin values where gauge relative weighting is less then
-  ! 0.5 (meaning raw IMERG was weighted at least 50% of final product)
-  ! do j = 1, col
-  !    do i = 1, row
-  !      if (wgtin(i,j) < 0.50) then
-  !      !if (wgtin(i,j) < 0.25) then
-  !      !if (wgtin(i,j) < 0.1) then
-  !      !if (wgtin(i,j) < 0.01) then
-  !         precipin(i,j) = LVT_rc%udef
-  !      end if
-  !   end do
-  ! end do
-
   ! Put the real(1:,1:) on the precipout(0:,0:)
   ! precipin is (ysize,xsize) starting at (lon=-179.9,lat=-89.9)
   precipout(1:xsize,1:ysize) = transpose(precipin)

From 5ea34c7a7e28f9cb039da6d8b20fa85fc5face61 Mon Sep 17 00:00:00 2001
From: Eric Kemp <eric.kemp@nasa.gov>
Date: Wed, 4 Oct 2023 15:34:37 -0400
Subject: [PATCH 06/16] First cut of script to calculate bias ratios for NAFPA.

Targets GFS and GALWEM background fields across multiple years, producing
12 monthly fields.  IMERG-FR V07A monthly data are used as "truth".

TODO:  Add code to move most currently hardwired settings to a config file.
---
 .../calc_gfsgalwem_biasratios_multiyear.py    | 189 ++++++++++++++++++
 1 file changed, 189 insertions(+)
 create mode 100755 lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py

diff --git a/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py b/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
new file mode 100755
index 000000000..8c0c7ed21
--- /dev/null
+++ b/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
@@ -0,0 +1,189 @@
+#!/usr/bin/env python3
+
+import datetime
+import os
+import sys
+
+from netCDF4 import Dataset as nc4_dataset
+import numpy as np
+
+# Start and end date of GFS period for NAFPA
+back_source = "GFS"
+#startdt = datetime.datetime(2008, 1, 1, 0, 0, 0) # GFS-Spectral
+#enddt   = datetime.datetime(2019, 6, 1, 0, 0, 0) # GFS-Spectral
+#startdt = datetime.datetime(2019, 7, 1, 0, 0, 0) # GFS-FV3
+#enddt   = datetime.datetime(2023, 4, 1, 0, 0, 0) # GFS-FV3
+#startdt = datetime.datetime(2008, 1, 1, 0, 0, 0) # All GFS
+#enddt   = datetime.datetime(2023, 4, 1, 0, 0, 0) # All GFS
+#startdt = datetime.datetime(2008, 1, 1, 0, 0, 0) # Pre-GALWEM
+#enddt   = datetime.datetime(2017, 9, 1, 0, 0, 0) # Pre-GALWEM
+startdt = datetime.datetime(2017, 10, 1, 0, 0, 0) # GALWEM Era
+enddt   = datetime.datetime(2023,  4, 1, 0, 0, 0) # GALWEM Era
+
+topdir_back = "/discover/nobackup/projects/usaf_lis/emkemp/AFWA/lis76_imergf_biascorr/data/GFS_NAFPA_Monthly_all"
+
+# Start and end date of GALWEM period for NAFPA
+#back_source = "GALWEM"
+#startdt = datetime.datetime(2017, 10, 1, 0, 0, 0)
+#enddt = datetime.datetime(2023, 4, 1, 0, 0, 0)
+#topdir_back = "/discover/nobackup/projects/usaf_lis/emkemp/AFWA/lis76_imergf_biascorr/data/GALWEM_NAFPA_Monthly_all"
+
+timedelta = datetime.timedelta(days=1)
+
+# Other data
+topdir_imergf = "/discover/nobackup/projects/usaf_lis/emkemp/AFWA/lis76_imergf_biascorr/data/IMERGF_V07A_NAFPA_Monthly"
+topdir_biasratio = f"testdir_{back_source}"
+
+def _main():
+    """Main driver"""
+
+    ny = 1920
+    nx = 2560
+    nmon = 12
+
+    sum_blended = np.zeros([nmon,ny,nx])
+    sum_back = np.zeros([nmon,ny,nx])
+    precip_ratio = np.zeros([nmon,ny,nx])
+
+    # Loop through each month
+    curdt = startdt
+    while curdt <= enddt:
+
+        # LVT output files have data from the prior month, so we must
+        # advance one month to find the appropriate file.  Python's
+        # timedelta object isn't smart enough to jump a whole month, so
+        # we loop through each day instead.
+        nextdt = curdt + timedelta
+        while nextdt.day != 1:
+            nextdt += timedelta
+
+        # First, the IMERG file
+        filename = f"{topdir_imergf}/SUM_TS."
+        filename += f"{nextdt.year:04d}{nextdt.month:02d}{nextdt.day:02d}"
+        filename += "0000.d01.nc"
+        ncid_imerg = nc4_dataset(filename, mode='r', \
+                                 format="NETCDF4_CLASSIC")
+        precip_imerg = ncid_imerg.variables['TotalPrecip'][:,:]
+        lats_imerg = ncid_imerg.variables["latitude"][:,:]
+        lons_imerg = ncid_imerg.variables["longitude"][:,:]
+        north_south = ncid_imerg.dimensions['north_south'].size
+        east_west = ncid_imerg.dimensions['east_west'].size
+        del ncid_imerg
+
+        # Next, the GFS or GALWEM file
+        filename = f"{topdir_back}/SUM_TS."
+        filename += f"{nextdt.year:04d}{nextdt.month:02d}{nextdt.day:02d}"
+        filename += "0000.d01.nc"
+        ncid_back = nc4_dataset(filename, mode='r', \
+                               format="NETCDF4_CLASSIC")
+        precip_back = ncid_back.variables['TotalPrecip'][:,:]
+        del ncid_back
+
+        # Use IMERG from 40S to 40N, use linear tapers from 40S to 60S,
+        # and 40N to 60N, and don't use IMERG poleward of 60N and 60S.
+        # RATIONALE:  No IR data is available poleward of 60N and 60S,
+        # and PMW data has diminished performance over frozen surfaces.
+        # Also, GPCC gage coverage tends to thin out in poleward regions,
+        # (no GPCC gages at all in Antarctica), so we screen out these
+        # areas for simplicity.  Finally, the 20 degree linear taper
+        # mimicks MERRA-2 usage of CPCU gauge analyses.
+        precip_imerg_weights = np.ones(np.shape(lats_imerg))
+        precip_imerg_weights = np.where(lats_imerg < -60,
+                                        0,
+                                        precip_imerg_weights)
+        precip_imerg_weights = np.where( (lats_imerg >= -60) &
+                                         (lats_imerg <  -40),
+                                         ( (lats_imerg + 60.) / 20.),
+                                         precip_imerg_weights)
+        #precip_imerg_weights = np.where( (lats_imerg >  40) &
+        #                                (lats_imerg <= 60),
+        #                                ( (60. - lats_imerg) / 20.),
+        #                                precip_imerg_weights)
+        #precip_imerg_weights = np.where(lats_imerg > 60,
+        #                              0,
+        #                              precip_imerg_weights)
+        # EMK Alternative...Move the northern linear taper region to
+        # 51N to 71N.  Rationale is to leverage the relatively high
+        # GPCC gauge density in the Scandinavian Peninsula.
+        precip_imerg_weights = np.where( (lats_imerg >  51) &
+                                          (lats_imerg <= 71),
+                                         ( (71. - lats_imerg) / 20.),
+                                          precip_imerg_weights)
+        precip_imerg_weights = np.where(lats_imerg > 71,
+                                        0,
+                                        precip_imerg_weights)
+
+        precip_blended = precip_imerg_weights[:,:]*precip_imerg[:,:] + \
+            (1. - precip_imerg_weights[:,:])*precip_back[:,:]
+
+        sum_blended[(curdt.month-1),:,:] += precip_blended[:,:] + 0.05
+        sum_back[(curdt.month-1),:,:] += precip_back[:,:] + 0.05
+
+        # Move on to next month
+        curdt = nextdt
+
+    # Output to file
+    filename = f"{topdir_biasratio}/{back_source}_pcp_biasratio.nc"
+
+    rootgrp = nc4_dataset(filename, "w", format="NETCDF4")
+    rootgrp.missing_value = np.float32("-9999.")
+    rootgrp.title = f"Monthly bias ratio for IMERG / {back_source}"
+    rootgrp.institution = "NASA GSFC"
+    now = datetime.datetime.utcnow()
+    history = "created on date: "
+    history += f"{now.year:04d}-{now.month:02d}-{now.day:02d}"
+    history += f"T{now.hour:02d}:{now.minute:02d}:{now.second:02d}"
+    rootgrp.history = f"created on date: {history}"
+    rootgrp.comment = "website: http://lis.gsfc.nasa.gov/"
+    rootgrp.MAP_PROJECTION = "EQUIDISTANT CYLINDRICAL"
+    rootgrp.SOUTH_WEST_CORNER_LAT = np.float32("-89.95312")
+    rootgrp.SOUTH_WEST_CORNER_LON = np.float32("-179.9297")
+    rootgrp.DX = np.float32("0.140625")
+    rootgrp.DY = np.float32("0.09375")
+
+    # Define dimensions
+    months = rootgrp.createDimension("months", 12)
+    north_south = rootgrp.createDimension("north_south", north_south)
+    east_west = rootgrp.createDimension("east_west", east_west)
+
+    # Define latitude
+    latitude = rootgrp.createVariable("latitude", "f4", \
+                                      ("north_south", "east_west",))
+    latitude.units = "degree_north"
+    latitude.standard_name = "latitude"
+    latitude.long_name = "latitude"
+    latitude.scale_factor = np.float32("1.")
+    latitude.add_offset = np.float32("0.")
+    latitude.missing_value = np.float32("-9999.")
+
+    # Define longitude
+    longitude = rootgrp.createVariable("longitude", "f4", \
+                                       ("north_south", "east_west",))
+    longitude.units = "degree_east"
+    longitude.standard_name = "longitude"
+    longitude.long_name = "longitude"
+    longitude.scale_factor = np.float32("1.")
+    longitude.add_offset = np.float32("0.")
+    longitude.missing_value = np.float32("-9999.")
+
+    # Define biasRatio
+    biasRatio = rootgrp.createVariable("biasRatio", "f4", \
+                              ("months", "north_south", "east_west",))
+    biasRatio.units = "-"
+    biasRatio.long_name = f"bias_ratio_for_{back_source}_precipitation"
+    biasRatio.scale_factor = np.float32("1.")
+    biasRatio.add_offset = np.float32("0.")
+    biasRatio.missing_value = np.float32("-9999.")
+
+    latitude[:,:] = lats_imerg[:,:]
+    longitude[:,:] = lons_imerg[:,:]
+
+    precip_ratio[:,:,:] = sum_blended[:,:,:] / sum_back[:,:,:]
+    precip_ratio[:,:,:] = np.where(precip_ratio == 0, 1, precip_ratio)
+    biasRatio[:,:,:] = precip_ratio[:,:,:]
+
+    rootgrp.close()
+
+
+if __name__ == "__main__":
+    _main()

From 138a2620bcc37abd271ac4752689aefab57a946f Mon Sep 17 00:00:00 2001
From: Eric Kemp <eric.kemp@nasa.gov>
Date: Wed, 4 Oct 2023 17:35:26 -0400
Subject: [PATCH 07/16] Revisions to use config file and command line
 arguments.

Some testing and further tweaks are still pending.
---
 .../calc_gfsgalwem_biasratios_multiyear.py    | 218 ++++++++++++++----
 1 file changed, 167 insertions(+), 51 deletions(-)

diff --git a/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py b/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
index 8c0c7ed21..9b98dd40f 100755
--- a/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
+++ b/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
@@ -1,5 +1,18 @@
 #!/usr/bin/env python3
 
+"""
+SCRIPT: calc_gfsgalwem_biasratios_multiyear.py
+
+Calculates bias ratio fields using monthly IMERG-FR V07A, GFS, and GALWEM
+data, all already interpolated to the NAFPA grid and summed to monthly
+totals via LVT.  Outputs 12 fields total, each covering a month based
+on multiple years of data.
+
+REVISION HISTORY:
+04 Oct 2023:  Eric Kemp:  Initial specification.
+"""
+
+import configparser
 import datetime
 import os
 import sys
@@ -34,8 +47,64 @@
 topdir_imergf = "/discover/nobackup/projects/usaf_lis/emkemp/AFWA/lis76_imergf_biascorr/data/IMERGF_V07A_NAFPA_Monthly"
 topdir_biasratio = f"testdir_{back_source}"
 
-def _main():
-    """Main driver"""
+def _usage():
+    """Print usage message for this script."""
+    print(f"Usage: {sys.argv[0]} CFGFILE BACKSOURCE STARTDATE ENDDATE")
+    print("   CFGFILE is path to config file")
+    print("   BACKSOURCE is GFS or GALWEM")
+    print("   STARTDATE is start date of calculation (YYYYMM)")
+    print("   ENDDATE is end date of calculateion (YYYYMM)")
+
+def _process_cmd_line():
+    """Process command line arguments."""
+    if len(sys.argv) != 5:
+        _usage()
+        sys.exit(1)
+
+    cfgfile = sys.argv[1]
+    if not os.path.exists(cfgfile):
+        print(f"[ERR] {cfgfile} does not exist!")
+        sys.exit(1)
+
+    backsource = sys.argv[2]
+    if backsource not in ["GFS", "GALWEM"]:
+        print(f"[ERR] {backsource} is not a valid background source")
+        print("[ERR] Only GFS and GALWEM recognized")
+        sys.exit(1)
+
+    yyyymm = sys.argv[3]
+    year = int(yyyymm[0:4])
+    month = int(yyyymm[4:6])
+    startdate = datetime.datetime(year, month, 1)
+
+    yyyymm = sys.argv[4]
+    year = int(yyyymm[0:4])
+    month = int(yyyymm[4:6])
+    enddate = datetime.datetime(year, month, 1)
+
+    if startdate > enddate:
+        print("[ERR] STARTDATE is beyond ENDDATE!")
+        sys.exit(1)
+
+    return cfgfile, backsource, startdate, enddate
+
+def _process_cfg_file(cfgfile, backsource):
+    """Processes config file for this script."""
+    config = configparser.ConfigParser()
+    config.read(cfgfile)
+
+    if backsource == "GFS":
+        backindir = config.get('Input', 'gfsdir')
+        backoutdir = config.get('Output', 'gfsdir')
+    elif backsource == "GALWEM":
+        backindir = config.get('Input', 'galwemdir')
+        backoutdir = config.get('Output', 'galwemdir')
+    imergdir = config.get('Input', 'imergdir')
+
+    return backindir, backoutdir, imergdir
+
+def _calc_biasratios(imergdir, backindir, startdate, enddate):
+    """Calculate the bias ratios from the input data files."""
 
     ny = 1920
     nx = 2560
@@ -45,22 +114,25 @@ def _main():
     sum_back = np.zeros([nmon,ny,nx])
     precip_ratio = np.zeros([nmon,ny,nx])
 
+    timedelta = datetime.timedelta(days=1)
+
     # Loop through each month
-    curdt = startdt
-    while curdt <= enddt:
+    curdate = startdate
+    while curdate <= enddate:
 
         # LVT output files have data from the prior month, so we must
-        # advance one month to find the appropriate file.  Python's
+        # advance one month to find the appropriate file, e.g., data
+        # for May 2020 will be in SUM_TS.202006010000.d01.nc.  Python's
         # timedelta object isn't smart enough to jump a whole month, so
         # we loop through each day instead.
-        nextdt = curdt + timedelta
-        while nextdt.day != 1:
-            nextdt += timedelta
+        nextdate = curdate + timedelta
+        while nextdate.day != 1:
+            nextdate += timedelta
 
         # First, the IMERG file
-        filename = f"{topdir_imergf}/SUM_TS."
-        filename += f"{nextdt.year:04d}{nextdt.month:02d}{nextdt.day:02d}"
-        filename += "0000.d01.nc"
+        filename = f"{imergdir}/SUM_TS."
+        filename += f"{nextdate.year:04d}{nextdate.month:02d}"
+        filename += f"{nextdate.day:02d}0000.d01.nc"
         ncid_imerg = nc4_dataset(filename, mode='r', \
                                  format="NETCDF4_CLASSIC")
         precip_imerg = ncid_imerg.variables['TotalPrecip'][:,:]
@@ -70,23 +142,23 @@ def _main():
         east_west = ncid_imerg.dimensions['east_west'].size
         del ncid_imerg
 
-        # Next, the GFS or GALWEM file
-        filename = f"{topdir_back}/SUM_TS."
-        filename += f"{nextdt.year:04d}{nextdt.month:02d}{nextdt.day:02d}"
-        filename += "0000.d01.nc"
+        # Next, the background data (GFS or GALWEM)
+        filename = f"{backindir}/SUM_TS."
+        filename += f"{nextdate.year:04d}{nextdate.month:02d}"
+        filename += f"{nextdate.day:02d}0000.d01.nc"
         ncid_back = nc4_dataset(filename, mode='r', \
                                format="NETCDF4_CLASSIC")
         precip_back = ncid_back.variables['TotalPrecip'][:,:]
         del ncid_back
 
-        # Use IMERG from 40S to 40N, use linear tapers from 40S to 60S,
-        # and 40N to 60N, and don't use IMERG poleward of 60N and 60S.
-        # RATIONALE:  No IR data is available poleward of 60N and 60S,
-        # and PMW data has diminished performance over frozen surfaces.
-        # Also, GPCC gage coverage tends to thin out in poleward regions,
-        # (no GPCC gages at all in Antarctica), so we screen out these
-        # areas for simplicity.  Finally, the 20 degree linear taper
-        # mimicks MERRA-2 usage of CPCU gauge analyses.
+        # Use IMERG from 40S to 71N; use linear tapers from 40S to 60S,
+        # and 51N to 71N; and don't use IMERG south of 60S and north of
+        # 71N.  Rationale:  No IR and little gauge data is available
+        # south of 60S, so we don't expect IMERG to be useful here.
+        # In the northern hemisphere, there is no IR data north of 60N,
+        # but there is good GPCC gauge density in Scandinavia up to 71N.
+        # Finally, the 20 degree latitude linear taper mimicks MERRA-2
+        # usage of CPCU gauge analyses.
         precip_imerg_weights = np.ones(np.shape(lats_imerg))
         precip_imerg_weights = np.where(lats_imerg < -60,
                                         0,
@@ -95,16 +167,6 @@ def _main():
                                          (lats_imerg <  -40),
                                          ( (lats_imerg + 60.) / 20.),
                                          precip_imerg_weights)
-        #precip_imerg_weights = np.where( (lats_imerg >  40) &
-        #                                (lats_imerg <= 60),
-        #                                ( (60. - lats_imerg) / 20.),
-        #                                precip_imerg_weights)
-        #precip_imerg_weights = np.where(lats_imerg > 60,
-        #                              0,
-        #                              precip_imerg_weights)
-        # EMK Alternative...Move the northern linear taper region to
-        # 51N to 71N.  Rationale is to leverage the relatively high
-        # GPCC gauge density in the Scandinavian Peninsula.
         precip_imerg_weights = np.where( (lats_imerg >  51) &
                                           (lats_imerg <= 71),
                                          ( (71. - lats_imerg) / 20.),
@@ -113,26 +175,58 @@ def _main():
                                         0,
                                         precip_imerg_weights)
 
+        # Conservative alternative:  Linear taper from 40N to 60N, and
+        # screen out everything north of 60N.  We'll do this as a backup
+        # if we find unphysical IMERG patterns north of 60N
+        #precip_imerg_weights = np.where( (lats_imerg >  40) &
+        #                                (lats_imerg <= 60),
+        #                                ( (60. - lats_imerg) / 20.),
+        #                                precip_imerg_weights)
+        #precip_imerg_weights = np.where(lats_imerg > 60,
+        #                              0,
+        #                              precip_imerg_weights)
+
         precip_blended = precip_imerg_weights[:,:]*precip_imerg[:,:] + \
             (1. - precip_imerg_weights[:,:])*precip_back[:,:]
 
-        sum_blended[(curdt.month-1),:,:] += precip_blended[:,:] + 0.05
-        sum_back[(curdt.month-1),:,:] += precip_back[:,:] + 0.05
+        # Add trace precipitation every month to prevent undefined
+        # ratios in deserts.
+        sum_blended[(curdate.month-1),:,:] += precip_blended[:,:] + 0.05
+        sum_back[(curdate.month-1),:,:] += precip_back[:,:] + 0.05
 
         # Move on to next month
-        curdt = nextdt
+        curdate = nextdate
 
-    # Output to file
-    filename = f"{topdir_biasratio}/{back_source}_pcp_biasratio.nc"
+    # Finish calculation
+    precip_ratio[:,:,:] = sum_blended[:,:,:] / sum_back[:,:,:]
+    precip_ratio[:,:,:] = np.where(precip_ratio == 0, 1, precip_ratio)
+
+    return lats_imerg, lons_imerg, precip_ratio, east_west, north_south
+
+def _create_output_filename(backoutdir, backsource, startdate, enddate):
+    """Create output netCDF file name"""
+    filename = f"{backoutdir}"
+    filename += f"/{backsource}_pcp_biasratios_"
+    filename += f"{startdate.year:04d}{startdate.month:02d}_"
+    filename += f"{enddate.year:04d}{enddate.month:02d}.nc"
+    return filename
+
+def _write_biasratios(args):
+    """Write out bias ratios to netCDF file"""
+
+    os.makedirs(args['backoutdir'], exist_ok=True)
 
-    rootgrp = nc4_dataset(filename, "w", format="NETCDF4")
-    rootgrp.missing_value = np.float32("-9999.")
-    rootgrp.title = f"Monthly bias ratio for IMERG / {back_source}"
-    rootgrp.institution = "NASA GSFC"
     now = datetime.datetime.utcnow()
+
     history = "created on date: "
     history += f"{now.year:04d}-{now.month:02d}-{now.day:02d}"
     history += f"T{now.hour:02d}:{now.minute:02d}:{now.second:02d}"
+
+    rootgrp = nc4_dataset(args['outfile'], "w", format="NETCDF4")
+    rootgrp.missing_value = np.float32("-9999.")
+    rootgrp.title = f"Monthly bias ratio for IMERG / {args['backsource']}"
+    rootgrp.institution = "NASA GSFC"
+
     rootgrp.history = f"created on date: {history}"
     rootgrp.comment = "website: http://lis.gsfc.nasa.gov/"
     rootgrp.MAP_PROJECTION = "EQUIDISTANT CYLINDRICAL"
@@ -143,8 +237,10 @@ def _main():
 
     # Define dimensions
     months = rootgrp.createDimension("months", 12)
-    north_south = rootgrp.createDimension("north_south", north_south)
-    east_west = rootgrp.createDimension("east_west", east_west)
+    north_south = rootgrp.createDimension("north_south", \
+                                          args['north_south'])
+    east_west = rootgrp.createDimension("east_west", \
+                                        args['east_west'])
 
     # Define latitude
     latitude = rootgrp.createVariable("latitude", "f4", \
@@ -170,20 +266,40 @@ def _main():
     biasRatio = rootgrp.createVariable("biasRatio", "f4", \
                               ("months", "north_south", "east_west",))
     biasRatio.units = "-"
-    biasRatio.long_name = f"bias_ratio_for_{back_source}_precipitation"
+    biasRatio.long_name = \
+        f"bias_ratio_for_{args['backsource']}_precipitation"
     biasRatio.scale_factor = np.float32("1.")
     biasRatio.add_offset = np.float32("0.")
     biasRatio.missing_value = np.float32("-9999.")
 
-    latitude[:,:] = lats_imerg[:,:]
-    longitude[:,:] = lons_imerg[:,:]
-
-    precip_ratio[:,:,:] = sum_blended[:,:,:] / sum_back[:,:,:]
-    precip_ratio[:,:,:] = np.where(precip_ratio == 0, 1, precip_ratio)
-    biasRatio[:,:,:] = precip_ratio[:,:,:]
+    latitude[:,:] = args['lats_imerg'][:,:]
+    longitude[:,:] = args['lons_imerg'][:,:]
+    biasRatio[:,:,:] = args['precip_ratio'][:,:,:]
 
     rootgrp.close()
 
+def _main():
+    """Main driver"""
+
+    cfgfile, backsource, startdate, enddate = _process_cmd_line()
+    backindir, backoutdir, imergdir = \
+        _process_cfg_file(cfgfile, backsource)
+    lats_imerg, lons_imerg, precip_ratio, east_west, north_south = \
+        _calc_biasratios(imergdir, backindir, startdate, enddate)
+    outfile = _create_output_filename(backoutdir, backsource,
+                                      startdate, enddate)
+    # To satisfy pylint....
+    args = {
+        "outfile" : outfile,
+        "backoutdir" : backoutdir,
+        "backsource" : backsource,
+        "north_south" : north_south,
+        "east_west" : east_west,
+        "lats_imerg" : lats_imerg,
+        "lons_imerg" : lons_imerg,
+        "precip_ratio" : precip_ratio,
+    }
+    _write_biasratios(args)
 
 if __name__ == "__main__":
     _main()

From 6927a6bfec0f12d79a4d388e2b5c10fd4b0e420e Mon Sep 17 00:00:00 2001
From: Eric Kemp <eric.kemp@nasa.gov>
Date: Thu, 5 Oct 2023 13:39:42 -0400
Subject: [PATCH 08/16] Revised code to only use single output directory.

---
 lvt/utils/usaf/nafpa/biasratio.cfg              | 17 +++++++++++++++++
 .../calc_gfsgalwem_biasratios_multiyear.py      | 17 ++++++++---------
 2 files changed, 25 insertions(+), 9 deletions(-)
 create mode 100644 lvt/utils/usaf/nafpa/biasratio.cfg

diff --git a/lvt/utils/usaf/nafpa/biasratio.cfg b/lvt/utils/usaf/nafpa/biasratio.cfg
new file mode 100644
index 000000000..c8ea41fa9
--- /dev/null
+++ b/lvt/utils/usaf/nafpa/biasratio.cfg
@@ -0,0 +1,17 @@
+#
+# Entries for calc_gfsgalwem_biasratios_multiyear.py
+#
+
+[Input]
+gfsdir: /discover/nobackup/projects/usaf_lis/emkemp/AFWA/lis76_imergf_biascorr/data/GFS_NAFPA_Monthly_all
+
+galwemdir: /discover/nobackup/projects/usaf_lis/emkemp/AFWA/lis76_imergf_biascorr/data/GALWEM_NAFPA_Monthly_all
+
+imergdir: /discover/nobackup/projects/usaf_lis/emkemp/AFWA/lis76_imergf_biascorr/data/IMERGF_V07A_NAFPA_Monthly
+
+[Output]
+outdir: biasratios
+
+
+
+
diff --git a/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py b/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
index 9b98dd40f..067b6a4f1 100755
--- a/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
+++ b/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
@@ -95,13 +95,12 @@ def _process_cfg_file(cfgfile, backsource):
 
     if backsource == "GFS":
         backindir = config.get('Input', 'gfsdir')
-        backoutdir = config.get('Output', 'gfsdir')
     elif backsource == "GALWEM":
         backindir = config.get('Input', 'galwemdir')
-        backoutdir = config.get('Output', 'galwemdir')
     imergdir = config.get('Input', 'imergdir')
+    outdir = config.get('Output', 'outdir')
 
-    return backindir, backoutdir, imergdir
+    return backindir, imergdir, outdir
 
 def _calc_biasratios(imergdir, backindir, startdate, enddate):
     """Calculate the bias ratios from the input data files."""
@@ -203,9 +202,9 @@ def _calc_biasratios(imergdir, backindir, startdate, enddate):
 
     return lats_imerg, lons_imerg, precip_ratio, east_west, north_south
 
-def _create_output_filename(backoutdir, backsource, startdate, enddate):
+def _create_output_filename(outdir, backsource, startdate, enddate):
     """Create output netCDF file name"""
-    filename = f"{backoutdir}"
+    filename = f"{outdir}"
     filename += f"/{backsource}_pcp_biasratios_"
     filename += f"{startdate.year:04d}{startdate.month:02d}_"
     filename += f"{enddate.year:04d}{enddate.month:02d}.nc"
@@ -214,7 +213,7 @@ def _create_output_filename(backoutdir, backsource, startdate, enddate):
 def _write_biasratios(args):
     """Write out bias ratios to netCDF file"""
 
-    os.makedirs(args['backoutdir'], exist_ok=True)
+    os.makedirs(args['outdir'], exist_ok=True)
 
     now = datetime.datetime.utcnow()
 
@@ -282,16 +281,16 @@ def _main():
     """Main driver"""
 
     cfgfile, backsource, startdate, enddate = _process_cmd_line()
-    backindir, backoutdir, imergdir = \
+    backindir, imergdir, outdir = \
         _process_cfg_file(cfgfile, backsource)
     lats_imerg, lons_imerg, precip_ratio, east_west, north_south = \
         _calc_biasratios(imergdir, backindir, startdate, enddate)
-    outfile = _create_output_filename(backoutdir, backsource,
+    outfile = _create_output_filename(outdir, backsource,
                                       startdate, enddate)
     # To satisfy pylint....
     args = {
         "outfile" : outfile,
-        "backoutdir" : backoutdir,
+        "outdir" : outdir,
         "backsource" : backsource,
         "north_south" : north_south,
         "east_west" : east_west,

From bd690dbd9b2ba715e059650987d6ecc50e19b671 Mon Sep 17 00:00:00 2001
From: Eric Kemp <eric.kemp@nasa.gov>
Date: Fri, 6 Oct 2023 09:14:15 -0400
Subject: [PATCH 09/16] Removed dead code.

---
 .../calc_gfsgalwem_biasratios_multiyear.py    | 27 -------------------
 1 file changed, 27 deletions(-)

diff --git a/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py b/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
index 067b6a4f1..be534d002 100755
--- a/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
+++ b/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
@@ -20,33 +20,6 @@
 from netCDF4 import Dataset as nc4_dataset
 import numpy as np
 
-# Start and end date of GFS period for NAFPA
-back_source = "GFS"
-#startdt = datetime.datetime(2008, 1, 1, 0, 0, 0) # GFS-Spectral
-#enddt   = datetime.datetime(2019, 6, 1, 0, 0, 0) # GFS-Spectral
-#startdt = datetime.datetime(2019, 7, 1, 0, 0, 0) # GFS-FV3
-#enddt   = datetime.datetime(2023, 4, 1, 0, 0, 0) # GFS-FV3
-#startdt = datetime.datetime(2008, 1, 1, 0, 0, 0) # All GFS
-#enddt   = datetime.datetime(2023, 4, 1, 0, 0, 0) # All GFS
-#startdt = datetime.datetime(2008, 1, 1, 0, 0, 0) # Pre-GALWEM
-#enddt   = datetime.datetime(2017, 9, 1, 0, 0, 0) # Pre-GALWEM
-startdt = datetime.datetime(2017, 10, 1, 0, 0, 0) # GALWEM Era
-enddt   = datetime.datetime(2023,  4, 1, 0, 0, 0) # GALWEM Era
-
-topdir_back = "/discover/nobackup/projects/usaf_lis/emkemp/AFWA/lis76_imergf_biascorr/data/GFS_NAFPA_Monthly_all"
-
-# Start and end date of GALWEM period for NAFPA
-#back_source = "GALWEM"
-#startdt = datetime.datetime(2017, 10, 1, 0, 0, 0)
-#enddt = datetime.datetime(2023, 4, 1, 0, 0, 0)
-#topdir_back = "/discover/nobackup/projects/usaf_lis/emkemp/AFWA/lis76_imergf_biascorr/data/GALWEM_NAFPA_Monthly_all"
-
-timedelta = datetime.timedelta(days=1)
-
-# Other data
-topdir_imergf = "/discover/nobackup/projects/usaf_lis/emkemp/AFWA/lis76_imergf_biascorr/data/IMERGF_V07A_NAFPA_Monthly"
-topdir_biasratio = f"testdir_{back_source}"
-
 def _usage():
     """Print usage message for this script."""
     print(f"Usage: {sys.argv[0]} CFGFILE BACKSOURCE STARTDATE ENDDATE")

From 88dd51755330494f2cdcc540d31e4be94f74ceaf Mon Sep 17 00:00:00 2001
From: Eric Kemp <eric.kemp@nasa.gov>
Date: Fri, 6 Oct 2023 09:43:13 -0400
Subject: [PATCH 10/16] Added log10 transform to output file for plotting.

---
 .../nafpa/calc_gfsgalwem_biasratios_multiyear.py     | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py b/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
index be534d002..1be185a43 100755
--- a/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
+++ b/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
@@ -10,6 +10,7 @@
 
 REVISION HISTORY:
 04 Oct 2023:  Eric Kemp:  Initial specification.
+06 Oct 2023:  Eric Kemp:  Added log10 transform to output for plotting.
 """
 
 import configparser
@@ -244,9 +245,20 @@ def _write_biasratios(args):
     biasRatio.add_offset = np.float32("0.")
     biasRatio.missing_value = np.float32("-9999.")
 
+    # Define log10BiasRatio
+    log10BiasRatio = rootgrp.createVariable("log10BiasRatio", "f4", \
+                                ("months", "north_south", "east_west",))
+    log10BiasRatio.units = "-"
+    log10BiasRatio.long_name = \
+        f"log10_bias_ratio_for_{args['backsource']}_precipitation"
+    log10BiasRatio.scale_factor = np.float32("1.")
+    log10BiasRatio.add_offset = np.float32("0.")
+    log10BiasRatio.missing_value = np.float32("-9999.")
+
     latitude[:,:] = args['lats_imerg'][:,:]
     longitude[:,:] = args['lons_imerg'][:,:]
     biasRatio[:,:,:] = args['precip_ratio'][:,:,:]
+    log10BiasRatio[:,:,:] = np.log10(args['precip_ratio'][:,:,:])
 
     rootgrp.close()
 

From 886df022763ec1611405d5ca56f49ba3c4a59dac Mon Sep 17 00:00:00 2001
From: Eric Kemp <eric.kemp@nasa.gov>
Date: Fri, 6 Oct 2023 10:15:13 -0400
Subject: [PATCH 11/16] Some refactoring to appease pylint.

---
 .../calc_gfsgalwem_biasratios_multiyear.py    | 131 ++++++++++--------
 1 file changed, 77 insertions(+), 54 deletions(-)

diff --git a/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py b/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
index 1be185a43..d30c92080 100755
--- a/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
+++ b/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
@@ -18,7 +18,11 @@
 import os
 import sys
 
+# Disable false alarm in pylint (Dataset is not visible in netCDF4 module
+# because it is not implemented in python)
+# pylint: disable=no-name-in-module
 from netCDF4 import Dataset as nc4_dataset
+# pylint: enable=no-name-in-module
 import numpy as np
 
 def _usage():
@@ -76,16 +80,20 @@ def _process_cfg_file(cfgfile, backsource):
 
     return backindir, imergdir, outdir
 
+def _create_arrays():
+    """Allocate and zero out numpy arrays"""
+    nlat = 1920
+    nlon = 2560
+    nmon = 12
+    sum_blended = np.zeros([nmon,nlat,nlon])
+    sum_back = np.zeros([nmon,nlat,nlon])
+    precip_ratio = np.zeros([nmon,nlat,nlon])
+    return sum_blended, sum_back, precip_ratio
+
 def _calc_biasratios(imergdir, backindir, startdate, enddate):
     """Calculate the bias ratios from the input data files."""
 
-    ny = 1920
-    nx = 2560
-    nmon = 12
-
-    sum_blended = np.zeros([nmon,ny,nx])
-    sum_back = np.zeros([nmon,ny,nx])
-    precip_ratio = np.zeros([nmon,ny,nx])
+    sum_blended, sum_back, precip_ratio = _create_arrays()
 
     timedelta = datetime.timedelta(days=1)
 
@@ -184,6 +192,54 @@ def _create_output_filename(outdir, backsource, startdate, enddate):
     filename += f"{enddate.year:04d}{enddate.month:02d}.nc"
     return filename
 
+def _create_latitude(rootgrp):
+    """Create latitude dataset in output file"""
+    latitude = rootgrp.createVariable("latitude", "f4", \
+                                      ("north_south", "east_west",))
+    latitude.units = "degree_north"
+    latitude.standard_name = "latitude"
+    latitude.long_name = "latitude"
+    latitude.scale_factor = np.float32("1.")
+    latitude.add_offset = np.float32("0.")
+    latitude.missing_value = np.float32("-9999.")
+    return latitude
+
+def _create_longitude(rootgrp):
+    """Create longitude dataset in output file."""
+    longitude = rootgrp.createVariable("longitude", "f4", \
+                                       ("north_south", "east_west",))
+    longitude.units = "degree_east"
+    longitude.standard_name = "longitude"
+    longitude.long_name = "longitude"
+    longitude.scale_factor = np.float32("1.")
+    longitude.add_offset = np.float32("0.")
+    longitude.missing_value = np.float32("-9999.")
+    return longitude
+
+def _create_bias_ratio(rootgrp, backsource):
+    """Create the biasRatio dataset in the output file."""
+    bias_ratio = rootgrp.createVariable("biasRatio", "f4", \
+                              ("months", "north_south", "east_west",))
+    bias_ratio.units = "-"
+    bias_ratio.long_name = \
+        f"bias_ratio_for_{backsource}_precipitation"
+    bias_ratio.scale_factor = np.float32("1.")
+    bias_ratio.add_offset = np.float32("0.")
+    bias_ratio.missing_value = np.float32("-9999.")
+    return bias_ratio
+
+def _create_log10_bias_ratio(rootgrp, backsource):
+    """Create the log10BiasRatio dataset in the output file."""
+    log10_bias_ratio = rootgrp.createVariable("log10BiasRatio", "f4", \
+                                ("months", "north_south", "east_west",))
+    log10_bias_ratio.units = "-"
+    log10_bias_ratio.long_name = \
+        f"log10_bias_ratio_for_{backsource}_precipitation"
+    log10_bias_ratio.scale_factor = np.float32("1.")
+    log10_bias_ratio.add_offset = np.float32("0.")
+    log10_bias_ratio.missing_value = np.float32("-9999.")
+    return log10_bias_ratio
+
 def _write_biasratios(args):
     """Write out bias ratios to netCDF file"""
 
@@ -209,56 +265,23 @@ def _write_biasratios(args):
     rootgrp.DY = np.float32("0.09375")
 
     # Define dimensions
-    months = rootgrp.createDimension("months", 12)
-    north_south = rootgrp.createDimension("north_south", \
-                                          args['north_south'])
-    east_west = rootgrp.createDimension("east_west", \
-                                        args['east_west'])
-
-    # Define latitude
-    latitude = rootgrp.createVariable("latitude", "f4", \
-                                      ("north_south", "east_west",))
-    latitude.units = "degree_north"
-    latitude.standard_name = "latitude"
-    latitude.long_name = "latitude"
-    latitude.scale_factor = np.float32("1.")
-    latitude.add_offset = np.float32("0.")
-    latitude.missing_value = np.float32("-9999.")
-
-    # Define longitude
-    longitude = rootgrp.createVariable("longitude", "f4", \
-                                       ("north_south", "east_west",))
-    longitude.units = "degree_east"
-    longitude.standard_name = "longitude"
-    longitude.long_name = "longitude"
-    longitude.scale_factor = np.float32("1.")
-    longitude.add_offset = np.float32("0.")
-    longitude.missing_value = np.float32("-9999.")
-
-    # Define biasRatio
-    biasRatio = rootgrp.createVariable("biasRatio", "f4", \
-                              ("months", "north_south", "east_west",))
-    biasRatio.units = "-"
-    biasRatio.long_name = \
-        f"bias_ratio_for_{args['backsource']}_precipitation"
-    biasRatio.scale_factor = np.float32("1.")
-    biasRatio.add_offset = np.float32("0.")
-    biasRatio.missing_value = np.float32("-9999.")
-
-    # Define log10BiasRatio
-    log10BiasRatio = rootgrp.createVariable("log10BiasRatio", "f4", \
-                                ("months", "north_south", "east_west",))
-    log10BiasRatio.units = "-"
-    log10BiasRatio.long_name = \
-        f"log10_bias_ratio_for_{args['backsource']}_precipitation"
-    log10BiasRatio.scale_factor = np.float32("1.")
-    log10BiasRatio.add_offset = np.float32("0.")
-    log10BiasRatio.missing_value = np.float32("-9999.")
+    rootgrp.createDimension("months", 12)
+    rootgrp.createDimension("north_south", \
+                            args['north_south'])
+    rootgrp.createDimension("east_west", \
+                            args['east_west'])
+
+    # Define output variables
+    latitude = _create_latitude(rootgrp)
+    longitude = _create_longitude(rootgrp)
+    bias_ratio = _create_bias_ratio(rootgrp, args['backsource'])
+    log10_bias_ratio = _create_log10_bias_ratio(rootgrp,
+                                                args['backsource'])
 
     latitude[:,:] = args['lats_imerg'][:,:]
     longitude[:,:] = args['lons_imerg'][:,:]
-    biasRatio[:,:,:] = args['precip_ratio'][:,:,:]
-    log10BiasRatio[:,:,:] = np.log10(args['precip_ratio'][:,:,:])
+    bias_ratio[:,:,:] = args['precip_ratio'][:,:,:]
+    log10_bias_ratio[:,:,:] = np.log10(args['precip_ratio'][:,:,:])
 
     rootgrp.close()
 

From f0bfb64a4c7882fe18b3414602761952c309dbb7 Mon Sep 17 00:00:00 2001
From: Eric Kemp <eric.kemp@nasa.gov>
Date: Fri, 6 Oct 2023 12:15:57 -0400
Subject: [PATCH 12/16] Refactored to appease pylint.

---
 .../calc_gfsgalwem_biasratios_multiyear.py    | 185 +++++++++++-------
 1 file changed, 114 insertions(+), 71 deletions(-)

diff --git a/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py b/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
index d30c92080..4d2e12b45 100755
--- a/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
+++ b/lvt/utils/usaf/nafpa/calc_gfsgalwem_biasratios_multiyear.py
@@ -11,6 +11,7 @@
 REVISION HISTORY:
 04 Oct 2023:  Eric Kemp:  Initial specification.
 06 Oct 2023:  Eric Kemp:  Added log10 transform to output for plotting.
+                          Refactored to please pylint.
 """
 
 import configparser
@@ -90,88 +91,127 @@ def _create_arrays():
     precip_ratio = np.zeros([nmon,nlat,nlon])
     return sum_blended, sum_back, precip_ratio
 
-def _calc_biasratios(imergdir, backindir, startdate, enddate):
+def _read_imerg_file(imergdir, nextdate):
+    """Read the interpolated monthly IMERG file."""
+    filename = f"{imergdir}/SUM_TS."
+    filename += f"{nextdate.year:04d}{nextdate.month:02d}"
+    filename += f"{nextdate.day:02d}0000.d01.nc"
+    ncid_imerg = nc4_dataset(filename, mode='r', \
+                             format="NETCDF4_CLASSIC")
+    precip_imerg = ncid_imerg.variables['TotalPrecip'][:,:]
+    del ncid_imerg
+    return precip_imerg
+
+def _read_imerg_file_latlondims(imergdir, startdate):
+    """Read the interpolated monthly IMERG file."""
+    nextdate = _get_nextdate(startdate)
+    filename = f"{imergdir}/SUM_TS."
+    filename += f"{nextdate.year:04d}{nextdate.month:02d}"
+    filename += f"{nextdate.day:02d}0000.d01.nc"
+    ncid_imerg = nc4_dataset(filename, mode='r', \
+                             format="NETCDF4_CLASSIC")
+    lats_imerg = ncid_imerg.variables["latitude"][:,:]
+    lons_imerg = ncid_imerg.variables["longitude"][:,:]
+    north_south = ncid_imerg.dimensions['north_south'].size
+    east_west = ncid_imerg.dimensions['east_west'].size
+    del ncid_imerg
+    return lats_imerg, lons_imerg, north_south, east_west
+
+def _read_background_file(backindir, nextdate):
+    """Read interpolated monthly background file."""
+    filename = f"{backindir}/SUM_TS."
+    filename += f"{nextdate.year:04d}{nextdate.month:02d}"
+    filename += f"{nextdate.day:02d}0000.d01.nc"
+    ncid_back = nc4_dataset(filename, mode='r', \
+                            format="NETCDF4_CLASSIC")
+    precip_back = ncid_back.variables['TotalPrecip'][:,:]
+    del ncid_back
+    return precip_back
+
+def _calc_precip_blended(lats_imerg, precip_imerg, precip_back):
+    """Calculate weighted blend of IMERG and background precip"""
+
+    # Use IMERG from 40S to 71N; use linear tapers from 40S to 60S,
+    # and 51N to 71N; and don't use IMERG south of 60S and north of
+    # 71N.  Rationale:  No IR and little gauge data is available
+    # south of 60S, so we don't expect IMERG to be useful here for
+    # bias correction.  In the northern hemisphere, there is no IR data
+    # north of 60N, but there is good GPCC gauge density in Scandinavia
+    # up to 71N.  Finally, the 20 degree latitude linear tapers mimicks
+    # MERRA-2 usage of CPCU gauge analyses.
+    precip_imerg_weights = np.ones(np.shape(lats_imerg))
+    precip_imerg_weights = np.where(lats_imerg < -60,
+                                    0,
+                                    precip_imerg_weights)
+    precip_imerg_weights = np.where( (lats_imerg >= -60) &
+                                     (lats_imerg <  -40),
+                                     ( (lats_imerg + 60.) / 20.),
+                                     precip_imerg_weights)
+    precip_imerg_weights = np.where( (lats_imerg >  51) &
+                                     (lats_imerg <= 71),
+                                     ( (71. - lats_imerg) / 20.),
+                                     precip_imerg_weights)
+    precip_imerg_weights = np.where(lats_imerg > 71,
+                                    0,
+                                    precip_imerg_weights)
+
+    # Conservative alternative:  Linear taper from 40N to 60N, and
+    # screen out everything north of 60N.  We'll do this as a backup
+    # if we find unphysical IMERG patterns north of 60N
+    #precip_imerg_weights = np.where( (lats_imerg >  40) &
+    #                                (lats_imerg <= 60),
+    #                                ( (60. - lats_imerg) / 20.),
+    #                                precip_imerg_weights)
+    #precip_imerg_weights = np.where(lats_imerg > 60,
+    #                              0,
+    #                              precip_imerg_weights)
+
+    precip_blended = precip_imerg_weights[:,:]*precip_imerg[:,:] + \
+        (1. - precip_imerg_weights[:,:])*precip_back[:,:]
+
+    return precip_blended
+
+def _get_nextdate(curdate):
+    """Determine next date (first of next month)"""
+    # LVT output files have data from the prior month, so we must
+    # advance one month to find the appropriate file, e.g., data
+    # for May 2020 will be in SUM_TS.202006010000.d01.nc.  Python's
+    # timedelta object isn't smart enough to jump a whole month, so
+    # we loop through each day instead.
+    timedelta = datetime.timedelta(days=1)
+    nextdate = curdate + timedelta
+    while nextdate.day != 1:
+        nextdate += timedelta
+    return nextdate
+
+def _calc_biasratios(imergdir, backindir, startdate, enddate, lats_imerg):
     """Calculate the bias ratios from the input data files."""
 
     sum_blended, sum_back, precip_ratio = _create_arrays()
 
-    timedelta = datetime.timedelta(days=1)
-
     # Loop through each month
     curdate = startdate
     while curdate <= enddate:
 
         # LVT output files have data from the prior month, so we must
         # advance one month to find the appropriate file, e.g., data
-        # for May 2020 will be in SUM_TS.202006010000.d01.nc.  Python's
-        # timedelta object isn't smart enough to jump a whole month, so
-        # we loop through each day instead.
-        nextdate = curdate + timedelta
-        while nextdate.day != 1:
-            nextdate += timedelta
+        # for May 2020 will be in SUM_TS.202006010000.d01.nc.
+        nextdate = _get_nextdate(curdate)
 
         # First, the IMERG file
-        filename = f"{imergdir}/SUM_TS."
-        filename += f"{nextdate.year:04d}{nextdate.month:02d}"
-        filename += f"{nextdate.day:02d}0000.d01.nc"
-        ncid_imerg = nc4_dataset(filename, mode='r', \
-                                 format="NETCDF4_CLASSIC")
-        precip_imerg = ncid_imerg.variables['TotalPrecip'][:,:]
-        lats_imerg = ncid_imerg.variables["latitude"][:,:]
-        lons_imerg = ncid_imerg.variables["longitude"][:,:]
-        north_south = ncid_imerg.dimensions['north_south'].size
-        east_west = ncid_imerg.dimensions['east_west'].size
-        del ncid_imerg
+        precip_imerg = \
+            _read_imerg_file(imergdir, nextdate)
 
         # Next, the background data (GFS or GALWEM)
-        filename = f"{backindir}/SUM_TS."
-        filename += f"{nextdate.year:04d}{nextdate.month:02d}"
-        filename += f"{nextdate.day:02d}0000.d01.nc"
-        ncid_back = nc4_dataset(filename, mode='r', \
-                               format="NETCDF4_CLASSIC")
-        precip_back = ncid_back.variables['TotalPrecip'][:,:]
-        del ncid_back
-
-        # Use IMERG from 40S to 71N; use linear tapers from 40S to 60S,
-        # and 51N to 71N; and don't use IMERG south of 60S and north of
-        # 71N.  Rationale:  No IR and little gauge data is available
-        # south of 60S, so we don't expect IMERG to be useful here.
-        # In the northern hemisphere, there is no IR data north of 60N,
-        # but there is good GPCC gauge density in Scandinavia up to 71N.
-        # Finally, the 20 degree latitude linear taper mimicks MERRA-2
-        # usage of CPCU gauge analyses.
-        precip_imerg_weights = np.ones(np.shape(lats_imerg))
-        precip_imerg_weights = np.where(lats_imerg < -60,
-                                        0,
-                                        precip_imerg_weights)
-        precip_imerg_weights = np.where( (lats_imerg >= -60) &
-                                         (lats_imerg <  -40),
-                                         ( (lats_imerg + 60.) / 20.),
-                                         precip_imerg_weights)
-        precip_imerg_weights = np.where( (lats_imerg >  51) &
-                                          (lats_imerg <= 71),
-                                         ( (71. - lats_imerg) / 20.),
-                                          precip_imerg_weights)
-        precip_imerg_weights = np.where(lats_imerg > 71,
-                                        0,
-                                        precip_imerg_weights)
-
-        # Conservative alternative:  Linear taper from 40N to 60N, and
-        # screen out everything north of 60N.  We'll do this as a backup
-        # if we find unphysical IMERG patterns north of 60N
-        #precip_imerg_weights = np.where( (lats_imerg >  40) &
-        #                                (lats_imerg <= 60),
-        #                                ( (60. - lats_imerg) / 20.),
-        #                                precip_imerg_weights)
-        #precip_imerg_weights = np.where(lats_imerg > 60,
-        #                              0,
-        #                              precip_imerg_weights)
-
-        precip_blended = precip_imerg_weights[:,:]*precip_imerg[:,:] + \
-            (1. - precip_imerg_weights[:,:])*precip_back[:,:]
-
-        # Add trace precipitation every month to prevent undefined
-        # ratios in deserts.
+        precip_back = _read_background_file(backindir, nextdate)
+
+        # Calculate blended precipitation (weighted average of IMERG
+        # and background, varying by latitude).
+        precip_blended = \
+            _calc_precip_blended(lats_imerg, precip_imerg, precip_back)
+
+        # Updated precip sums.  Add trace precipitation every month to
+        # prevent undefined ratios in deserts.
         sum_blended[(curdate.month-1),:,:] += precip_blended[:,:] + 0.05
         sum_back[(curdate.month-1),:,:] += precip_back[:,:] + 0.05
 
@@ -182,7 +222,7 @@ def _calc_biasratios(imergdir, backindir, startdate, enddate):
     precip_ratio[:,:,:] = sum_blended[:,:,:] / sum_back[:,:,:]
     precip_ratio[:,:,:] = np.where(precip_ratio == 0, 1, precip_ratio)
 
-    return lats_imerg, lons_imerg, precip_ratio, east_west, north_south
+    return precip_ratio
 
 def _create_output_filename(outdir, backsource, startdate, enddate):
     """Create output netCDF file name"""
@@ -291,8 +331,11 @@ def _main():
     cfgfile, backsource, startdate, enddate = _process_cmd_line()
     backindir, imergdir, outdir = \
         _process_cfg_file(cfgfile, backsource)
-    lats_imerg, lons_imerg, precip_ratio, east_west, north_south = \
-        _calc_biasratios(imergdir, backindir, startdate, enddate)
+    lats_imerg, lons_imerg, north_south, east_west = \
+        _read_imerg_file_latlondims(imergdir, startdate)
+    precip_ratio = \
+        _calc_biasratios(imergdir, backindir, startdate, enddate, \
+                         lats_imerg)
     outfile = _create_output_filename(outdir, backsource,
                                       startdate, enddate)
     # To satisfy pylint....

From a18cca2df55f6f2cb2785bb34c24c804dd445654 Mon Sep 17 00:00:00 2001
From: Eric Kemp <eric.kemp@nasa.gov>
Date: Fri, 6 Oct 2023 16:27:30 -0400
Subject: [PATCH 13/16] Added print statements for when bias ratio data are
 read in.

Also, minor code clean up.
---
 lis/metforcing/usaf/USAF_bratsethMod.F90 | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/lis/metforcing/usaf/USAF_bratsethMod.F90 b/lis/metforcing/usaf/USAF_bratsethMod.F90
index 3995698ec..cc4754641 100644
--- a/lis/metforcing/usaf/USAF_bratsethMod.F90
+++ b/lis/metforcing/usaf/USAF_bratsethMod.F90
@@ -6464,9 +6464,7 @@ subroutine USAF_pcpBackBiasRatio_s2s(n, yyyymmddhh)
 
      ! Imports
      use AGRMET_forcingMod, only: agrmet_struc
-     use LIS_coreMod, only: LIS_rc!, LIS_localPet, &
-          !LIS_ews_halo_ind, LIS_ewe_halo_ind, &
-          !LIS_nss_halo_ind, LIS_nse_halo_ind
+     use LIS_coreMod, only: LIS_rc
      use LIS_logMod, only: LIS_logunit, LIS_endrun, LIS_verify
 #if ( defined USE_NETCDF3 || defined USE_NETCDF4 )
      use netcdf
@@ -6598,6 +6596,9 @@ subroutine USAF_pcpBackBiasRatio_nrt(n, yyyymmddhh)
           start=(/1,1,imonth/), &
           count=(/LIS_rc%gnc(n),LIS_rc%gnr(n),1/)), &
           'n90_get_var failed for biasRatio in LIS param file')
+     write(LIS_logunit,*) &
+          '[INFO] Read GFS bias ratios for month ', imonth, &
+          ' from ', trim(agrmet_struc(n)%gfs_nrt_bias_ratio_file)
      call LIS_verify(nf90_close(ncid),&
           'nf90_close failed in USAF_pcpBackBiasRatio_nrt')
 
@@ -6620,6 +6621,9 @@ subroutine USAF_pcpBackBiasRatio_nrt(n, yyyymmddhh)
           start=(/1,1,imonth/), &
           count=(/LIS_rc%gnc(n),LIS_rc%gnr(n),1/)),&
           'n90_get_var failed for biasRatio in LIS param file')
+     write(LIS_logunit,*) &
+          '[INFO] Read GALWEM bias ratios for month ', imonth, &
+          ' from ', trim(agrmet_struc(n)%galwem_nrt_bias_ratio_file)
      call LIS_verify(nf90_close(ncid),&
           'nf90_close failed in USAF_pcpBackBiasRatio_nrt')
 

From b6cd67f4d00649b80af64a38f68423ed66dba139 Mon Sep 17 00:00:00 2001
From: Eric Kemp <eric.kemp@nasa.gov>
Date: Tue, 10 Oct 2023 14:04:48 -0400
Subject: [PATCH 14/16] Updated LVT documentation to add "GPM IMERG Monthly"
 datastream.

---
 lvt/configs/lvt.config.adoc | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/lvt/configs/lvt.config.adoc b/lvt/configs/lvt.config.adoc
index 9ba94fea5..0c3f15a93 100644
--- a/lvt/configs/lvt.config.adoc
+++ b/lvt/configs/lvt.config.adoc
@@ -118,6 +118,7 @@ LVT output methodology:                 "2d gridspace"
 | "`ANSA SWE`"                      | ANSA SWE retrievals
 | "`CPC precipitation`"             | CPC unified precipitation product
 | "`GPM IMERG`"                     | IMERG precipitation data product
+| "`GPM IMERG Monthly`"             | IMERG-Monthly precipitation product
 | "`USGS streamflow`"               | USGS streamflow observations
 | "`Naturalized streamflow`"        | Naturalized streamflow estimates
 | "`FLUXNET MTE`"                   | Gridded FLUXNET MTE data from MPI
@@ -1276,6 +1277,21 @@ IMERG version:                   V06B
 IMERG product:                   final
 ....
 
+==== GPM IMERG Monthly precipitation data
+
+`IMERG monthly data directory:` specifies the location of the GPM IMERG Monthly precipitation product
+
+`IMERG monthly version:` specifies the version of the GPM IMERG Monthly precipitation product. Most current version is V07A.
+
+`IMERG monthly product:` specifies the GPM IMERG Monthly precipitation product. Current option is: final.  (early and late will be added later)
+
+.Example _lvt.config_ entry
+....
+IMERG monthly data directory:            input/IMERGF_V07A_Monthly
+IMERG monthly version:                   V07A
+IMERG monthly product:                   final
+....
+
 ==== ESA CCI soil moisture data
 
 `ESA CCI soil moisture data directory:` specifies the location of the ESA CCI soil moisture data

From 94aac20570b6ea83ef0fbfd0e8b201c2e65489e3 Mon Sep 17 00:00:00 2001
From: Eric Kemp <eric.kemp@nasa.gov>
Date: Mon, 11 Dec 2023 15:56:28 -0500
Subject: [PATCH 15/16] Tweaks to declaring SLEEP and MPI calls, to pacify GNU
 compiler.

---
 lis/metforcing/usaf/USAF_bratsethMod.F90 | 22 ----------------------
 lis/metforcing/usaf/readcrd_agrmet.F90   |  5 +----
 2 files changed, 1 insertion(+), 26 deletions(-)

diff --git a/lis/metforcing/usaf/USAF_bratsethMod.F90 b/lis/metforcing/usaf/USAF_bratsethMod.F90
index cc4754641..883e79acf 100644
--- a/lis/metforcing/usaf/USAF_bratsethMod.F90
+++ b/lis/metforcing/usaf/USAF_bratsethMod.F90
@@ -777,8 +777,6 @@ subroutine USAF_getBackNWP(nest,back4,pcp_src, use_twelve, j6hr, findex)
       integer :: c, r
 
       external :: AGRMET_julhr_date10
-      external :: MPI_Barrier
-      external :: sleep
 
       TRACE_ENTER("bratseth_getBackNWP")
       rc = 0
@@ -1872,8 +1870,6 @@ subroutine calc_invDataDensities(this,sigmaBSqr,nest,max_dist, &
       integer :: imax,jmax
       logical :: verbose
 
-      external :: MPI_Barrier, MPI_ALLREDUCE
-
       verbose = .true.
       if (present(silent)) then
          if (silent) verbose = .false.
@@ -2098,8 +2094,6 @@ subroutine calc_obsAnalysis(this,sigmaBSqr,nobs,invDataDensities,nest,&
       real :: O, A
       integer :: good_obs
 
-      external :: MPI_Barrier, MPI_ALLREDUCE
-
       verbose = .true.
       if (present(silent)) then
          if (silent) verbose = .false.
@@ -2551,8 +2545,6 @@ subroutine calc_gridAnalysis(this,nest,sigmaBSqr,nobs,invDataDensities,&
       integer :: r_local, c_local
       integer :: gindex
 
-      external :: MPI_Barrier, MPI_ALLREDUCE
-
       ! Sanity checks
       if (nobs .eq. 0) return
 
@@ -3254,8 +3246,6 @@ subroutine USAF_superstatQC(this,nest,new_name,network,silent_rejects)
       double precision :: t1, t2
       logical :: silent_rejects_local
 
-      external :: MPI_Barrier, MPI_Allreduce
-
       ! Sanity check
       nobs = this%nobs
       if (nobs .eq. 0) then
@@ -3670,8 +3660,6 @@ subroutine USAF_dupQC(this)
       integer :: ierr
       logical :: location_issue
 
-      external :: MPI_Barrier
-
       nobs = this%nobs
       if (nobs .eq. 0) then
          write(LIS_logunit,*)&
@@ -4061,8 +4049,6 @@ subroutine USAF_backQC(this,sigmaBSqr,silent_rejects)
       double precision :: t1, t2
       logical :: silent_rejects_local
 
-      external :: MPI_Barrier
-      
       nobs = this%nobs
       if (nobs .eq. 0) then
          write(LIS_logunit,*)&
@@ -4839,8 +4825,6 @@ subroutine find_LIS_cols_rows(this,nest,cols,rows)
       integer :: j,r,c,gindex
       integer :: ierr
 
-      external :: MPI_Barrier, MPI_ALLREDUCE
-
       nobs = this%nobs
       if (nobs .eq. 0) return
 
@@ -4987,8 +4971,6 @@ subroutine USAF_waterQC(this,nest,silent_rejects)
       integer :: reject_count
       logical :: silent_rejects_local
 
-      external :: MPI_Barrier
-
       ! Sanity check
       nobs = this%nobs
       if (nobs .eq. 0) then
@@ -5091,8 +5073,6 @@ subroutine USAF_snowQC(this,nest,hourindex,threshold,silent_rejects)
       real :: threshold_local
       logical :: silent_rejects_local
 
-      external :: MPI_Barrier, MPI_ALLREDUCE
-
       ! Sanity check
       nobs = this%nobs
       if (nobs .eq. 0) then
@@ -5276,8 +5256,6 @@ subroutine USAF_snowDepthQC(this,nest,silent_rejects)
       integer :: rglb,cglb
       integer :: gid
 
-      external :: MPI_Barrier, MPI_ALLREDUCE
-
       ! Sanity check
       nobs = this%nobs
       if (nobs .eq. 0) then
diff --git a/lis/metforcing/usaf/readcrd_agrmet.F90 b/lis/metforcing/usaf/readcrd_agrmet.F90
index 819905484..47635b67f 100644
--- a/lis/metforcing/usaf/readcrd_agrmet.F90
+++ b/lis/metforcing/usaf/readcrd_agrmet.F90
@@ -40,7 +40,7 @@ subroutine readcrd_agrmet()
   use LIS_logMod,     only : LIS_logunit, LIS_verify, LIS_abort, &
        LIS_endrun
 #if (defined SPMD)
-  use LIS_mpiMod, only: LIS_MPI_COMM
+  use LIS_mpiMod
 #endif
   use LIS_pluginIndices, only : LIS_agrmetrunId
   use AGRMET_forcingMod, only : agrmet_struc
@@ -64,9 +64,6 @@ subroutine readcrd_agrmet()
   integer :: ierr
   logical :: use_nrt_bias_files ! EMK
 
-  external :: MPI_Barrier
-  external :: sleep
-
   call ESMF_ConfigFindLabel(LIS_config,"AGRMET forcing directory:",rc=rc)
   do n=1,LIS_rc%nnest
      call ESMF_ConfigGetAttribute(LIS_config,agrmet_struc(n)%agrmetdir,rc=rc)

From 0b2b0727218f37a61398f2a9592e8f4441d540c1 Mon Sep 17 00:00:00 2001
From: Eric Kemp <eric.kemp@nasa.gov>
Date: Mon, 11 Dec 2023 16:45:53 -0500
Subject: [PATCH 16/16] Updated character string lengths to use
 LIS_CONST_PATH_LEN.

---
 lis/metforcing/usaf/AGRMET_forcingMod.F90 | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/lis/metforcing/usaf/AGRMET_forcingMod.F90 b/lis/metforcing/usaf/AGRMET_forcingMod.F90
index ad1f1f809..220ca1563 100644
--- a/lis/metforcing/usaf/AGRMET_forcingMod.F90
+++ b/lis/metforcing/usaf/AGRMET_forcingMod.F90
@@ -302,7 +302,9 @@ module AGRMET_forcingMod
 !   contains the previous and next forcing values, respectively
 !  \end{description}
 !
-! !USES: 
+! !USES:
+  use LIS_constantsMod, only : LIS_CONST_PATH_LEN
+
   implicit none
 
   PRIVATE
@@ -721,8 +723,8 @@ module AGRMET_forcingMod
      real, allocatable :: pcp_back_bias_ratio(:,:)
      ! EMK Add NRT bias correction toward IMERG-Final Run
      ! (back_bias_corr == 2)
-     character(255) :: gfs_nrt_bias_ratio_file
-     character(255) :: galwem_nrt_bias_ratio_file
+     character(LIS_CONST_PATH_LEN) :: gfs_nrt_bias_ratio_file
+     character(LIS_CONST_PATH_LEN) :: galwem_nrt_bias_ratio_file
      real, allocatable :: gfs_nrt_bias_ratio(:,:)
      real, allocatable :: galwem_nrt_bias_ratio(:,:)
      integer :: pcp_back_bias_ratio_month