fix for 4diau with iau_filter_increments=T (#167)

* fix for 4diau with iau_filter_increments=T

* fix time interval for 3DIAU

* fix typo in comment

* fix bug found in review by @MingjingTong-NOAA

* change tnext to integer variable itnext

Co-authored-by: jswhit2 <[email protected]>

tools/fv_iau_mod.F90

@@ -253,6 +253,7 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm)
  allocate (iau_state%inc1%tracer_inc(is:ie, js:je, km,ntracers))
  iau_state%hr1=IPD_Control%iaufhrs(1)
  iau_state%wt = 1.0 ! IAU increment filter weights (default 1.0)
+ iau_state%wt_normfact = 1.0
  if (IPD_Control%iau_filter_increments) then
  ! compute increment filter weights, sum to obtain normalization factor
  dtp=IPD_control%dtp
@@ -298,26 +299,31 @@ subroutine getiauforcing(IPD_Control,IAU_Data)
  type (IPD_control_type), intent(in) :: IPD_Control
  type(IAU_external_data_type), intent(inout) :: IAU_Data
  real(kind=kind_phys) t1,t2,sx,wx,wt,dtp
- integer n,i,j,k,sphum,kstep,nstep
+ integer n,i,j,k,sphum,kstep,nstep,itnext
 
  IAU_Data%in_interval=.false.
  if (nfiles.LE.0) then
  return
  endif
 
- t1=iau_state%hr1 - IPD_Control%iau_delthrs*0.5
- t2=iau_state%hr1 + IPD_Control%iau_delthrs*0.5
+ if (nfiles .eq. 1) then 
+ t1 = IPD_Control%iaufhrs(1)-0.5*IPD_Control%iau_delthrs
+ t2 = IPD_Control%iaufhrs(1)+0.5*IPD_Control%iau_delthrs
+ else
+ t1 = IPD_Control%iaufhrs(1)
+ t2 = IPD_Control%iaufhrs(nfiles)
+ endif
  if (IPD_Control%iau_filter_increments) then
  ! compute increment filter weight
- ! t1 beginning of window, t2 end of window
+ ! t1 is beginning of window, t2 end of window
  ! IPD_Control%fhour current time
  ! in window kstep=-nstep,nstep (2*nstep+1 total)
  ! time step IPD_control%dtp
  dtp=IPD_control%dtp
  nstep = 0.5*IPD_Control%iau_delthrs*3600/dtp
  ! compute normalized filter weight
- kstep = (IPD_Control%fhour-(t1+IPD_Control%iau_delthrs*0.5))*3600./dtp
- if (kstep .ge. -nstep .and. kstep .le. nstep) then
+ kstep = ((IPD_Control%fhour-t1) - 0.5*IPD_Control%iau_delthrs)*3600./dtp
+ if (IPD_Control%fhour >= t1 .and. IPD_Control%fhour < t2) then
  sx = acos(-1.)*kstep/nstep
  wx = acos(-1.)*kstep/(nstep+1)
  if (kstep .ne. 0) then
@@ -326,7 +332,7 @@ subroutine getiauforcing(IPD_Control,IAU_Data)
  wt = 1.
  endif
  iau_state%wt = iau_state%wt_normfact*wt
- if (is_master()) print *,'filter wt',kstep,IPD_Control%fhour,iau_state%wt
+ !if (is_master()) print *,'kstep,t1,t,t2,filter wt=',kstep,t1,IPD_Control%fhour,t2,iau_state%wt/iau_state%wt_normfact
  else
  iau_state%wt = 0.
  endif
@@ -340,31 +346,32 @@ subroutine getiauforcing(IPD_Control,IAU_Data)
  IAU_Data%in_interval=.false.
  else
  if (IPD_Control%iau_filter_increments) call setiauforcing(IPD_Control,IAU_Data,iau_state%wt)
- if (is_master()) print *,'apply iau forcing',t1,IPD_Control%fhour,t2
+ if (is_master()) print *,'apply iau forcing t1,t,t2,filter wt=',t1,IPD_Control%fhour,t2,iau_state%wt/iau_state%wt_normfact
  IAU_Data%in_interval=.true.
  endif
  return
  endif
 
  if (nfiles > 1) then
- t2=2
- if (IPD_Control%fhour < IPD_Control%iaufhrs(1) .or. IPD_Control%fhour >= IPD_Control%iaufhrs(nfiles)) then
+ itnext=2
+ if (IPD_Control%fhour < t1 .or. IPD_Control%fhour >= t2) then
 ! if (is_master()) print *,'no iau forcing',IPD_Control%iaufhrs(1),IPD_Control%fhour,IPD_Control%iaufhrs(nfiles)
  IAU_Data%in_interval=.false.
  else
+ if (is_master()) print *,'apply iau forcing t1,t,t2,filter wt=',t1,IPD_Control%fhour,t2,iau_state%wt/iau_state%wt_normfact
  IAU_Data%in_interval=.true.
  do k=nfiles,1,-1
  if (IPD_Control%iaufhrs(k) > IPD_Control%fhour) then
- t2=k
+ itnext=k
  endif
  enddo
-! if (is_master()) print *,'t2=',t2
+! if (is_master()) print *,'itnext=',itnext
  if (IPD_Control%fhour >= iau_state%hr2) then ! need to read in next increment file
  iau_state%hr1=iau_state%hr2
- iau_state%hr2=IPD_Control%iaufhrs(int(t2))
+ iau_state%hr2=IPD_Control%iaufhrs(itnext)
  iau_state%inc1=iau_state%inc2
- if (is_master()) print *,'reading next increment file',trim(IPD_Control%iau_inc_files(int(t2)))
- call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(int(t2))))
+ if (is_master()) print *,'reading next increment file',trim(IPD_Control%iau_inc_files(itnext))
+ call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(itnext)))
  endif
  call updateiauforcing(IPD_Control,IAU_Data,iau_state%wt)
  endif