Possible bug in obs_write.F90

[hjsong]

Hi, I would like to report the possible bug in computing the misfit in obs_write.F90, which is

Code: Select all

!-----------------------------------------------------------------------
!  Compute and write initial and final model-observation misfit
!  (innovation) vector for output purposes only. Write also initial
!  nonlinear model at observation locations.
!-----------------------------------------------------------------------
!
        IF (wrtMisfit(ng)) THEN
          DO iobs=Mstr,Mend
#  if defined IS4DVAR
            misfit(iobs)=ObsScale(iobs)*SQRT(ObsErr(iobs))*             &
     &                   (NLmodVal(iobs)+TLmodVal(iobs)-ObsVal(iobs))

The computation for the misfit looks right for the inner loops.
But I think the 'TLmodVal' term should be taken out for the outer loops, especially for the very last nonlinear run.

I did a very simple test, which is just a one day assimilation using a single mooring data.
In order to fit the data perfectly, I inflated the model standard devation on purpose.
So the ROMS should yield very small misfit.
When I looked at the variable named 'misfit_final', I was very surprised and thought that something is very wrong.
But it turned out the way how ROMS compute the misfit.
After taking out the 'TLmodVal' from the misfit, everything looks perfect.

Please see this figure.



(Blue, red and black lines represent the initial misfit, final misfit from the calculation above (ROMS) and corrected misfit in which the 'TLmodVal' is removed.)

[arango]

Indeed, this is a bug. Good catch I corrected the problem and made few updates. Please update and check the following ticket for more details.

Thank you for reporting this bug and shown it so convincingly with the plots.

[hjsong]

Hi Hernan,

Thank you for the update.
I updated the code and test it, and I think I found another bug regarding this.

The misfit is not computed as "NLmodVal(iobs)-ObsVal(iobs)".

According to the 'is4dvar_ocean.h', misfit is saved when

Code: Select all

(line 548 - 555)
!
!  Activate switch to write out initial and final misfit between
!  model and observations.
!
            wrtMisfit(ng)=.FALSE.
            IF (((outer.eq.1).and.(inner.eq.0)).or.                     &
     &          ((outer.eq.Nouter).and.(inner.eq.Ninner))) THEN
              wrtMisfit(ng)=.TRUE.
            END IF

This will save the misfit twice - at the very first TL model run, and the very last TL model run.

After all inner loops, the flags are reset before the last NL model run.

Code: Select all

(line 1089 - 1098)
!
!  Initialize nonlinear model with estimated initial conditions.
!
        wrtNLmod(ng)=.TRUE.
        wrtTLmod(ng)=.FALSE.
        wrtMisfit(ng)=.FALSE.
        tRSTindx(ng)=0
        NrecRST(ng)=0
        CALL initial (ng)
        IF (exit_flag.ne.NoError) RETURN

Thus, the misfit is not save at the very last NL model run, and the 'misfit_final' in MOD file is not the final misfit.

Especially when Nouter=1, the "misfit_initial" and the "misfit_final" are the same when they are "NLmodVal(iobs)-ObsVal(iobs)" because "NLmodVal(iobs)" is the very first NL model value.

So I think the line 1094

Code: Select all

        wrtMisfit(ng)=.FALSE.

should be commented out unless this causes other problems.

Thank you.

Best,
Hajoon