From 5193c6b60c02c744b1ffe27078eccbeed2a22ad8 Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Fri, 4 Feb 2022 08:56:14 -0700 Subject: [PATCH] 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 --- tools/fv_iau_mod.F90 | 37 ++++++++++++++++++++++--------------- 1 file changed, 22 insertions(+), 15 deletions(-) diff --git a/tools/fv_iau_mod.F90 b/tools/fv_iau_mod.F90 index 12b87e70f..703640a4f 100644 --- a/tools/fv_iau_mod.F90 +++ b/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