Hi there,
I ran into some problem when trying to estimate the posterior errors in data assimilation WC13/PSAS case. I'm using the latest version of ROMS while adding two cpp options "POSTERIOR_EOFS" and "POSTERIOR_ERROR_I" to estimate posterior error. An error occurred after the DA run when trying to compute the Posterior Error Covariance Matrix as follow:
"
<<<< Full Posterior Error Covariance Matrix >>>>
Found Error: ** Line: 849 Source: ROMS/Utility/nf_fread2d.F
Found Error: ** Line: 2167 Source: ROMS/Utility/posterior.F
READ_STATE - error while reading variable: zeta at time record = 25
in NetCDF file: ocean_hss.nc
Found Error: 03 Line: 2155 Source: ROMS/Utility/posterior.F
"
I checked the useast_hss.nc file which contains 24 time records instead of 25 (Nouter=1; Ninner=25). But the Posterior Error analysis code is trying to read the 25th time record. After checking the model output log file, I found the model only started to write Hessian matrix values from the 2nd iteration showing below:
"
TL_WRT_INI - wrote inner-loop initial fields (Iter=0002, Index=1,1, Rec=0001)
WRT_HESSIAN - wrote Hessian fields (Index=1,1) in record = 0000001
"
Therefore, the useast_hss.nc file only has 24 records instead of 25. While trying to fix this problem, I found some suspicious code in ROMS/Utility/convolve.F as follow (also attached the whole file):
"
430 ! Write out tangent linear model initial conditions and tangent
431 ! linear surface forcing adjustments for next inner loop into
432 ! ITL(ng)%name (record Rec1). The tangent model initial
433 ! conditions are set to the convolved adjoint solution.
434 !
435 IF (model.eq.iTLM) THEN
436 DO ng=1,Ngrids
437 CALL tl_wrt_ini (ng, Rold(ng), Rec1)
438 IF (FoundError(exit_flag, NoError, __LINE__, &
439 & __FILE__)) RETURN
440 END DO
441 # if defined ARRAY_MODES || defined W4DVAR || \
442 defined W4DVAR_SENSITIVITY
443 ELSE IF (model.eq.iRPM) THEN
444 DO ng=1,Ngrids
445 CALL rp_wrt_ini (ng, Rold(ng), Rec2)
446 IF (FoundError(exit_flag, NoError, __LINE__, &
447 & __FILE__)) RETURN
448 END DO
449 # elif (defined W4DPSAS || defined W4DPSAS_SENSITIVITY) && \
450 !defined RPCG
451 ELSE IF (model.eq.iNLM) THEN
452 DO ng=1,Ngrids
453 CALL wrt_ini (ng, Rnew(ng))
454 IF (FoundError(exit_flag, NoError, __LINE__, &
455 & __FILE__)) RETURN
456 # if defined ADJUST_STFLUX || defined ADJUST_WSTRESS || \
457 defined ADJUST_BOUNDARY
458 CALL wrt_frc_AD (ng, Rold(ng), INI(ng)%Rindex)
459 IF (FoundError(exit_flag, NoError, __LINE__, &
460 & __FILE__)) RETURN
461 # endif
462 END DO
463 # endif
464 END IF
465 # ifdef RPCG
466 !
467 ! Write out tangent linear model initial conditions and tangent
468 ! linear surface forcing and obc adjustments for next outer
469 ! loop into ITLname (record Rec2). The tangent model initial
470 ! conditions are set to the convolved adjoint solution.
471 !
472 IF (model.eq.iNLM) THEN
473 DO ng=1,Ngrids
474 CALL tl_wrt_ini (ng, Rold(ng), Rec2)
475 IF (FoundError(exit_flag, NoError, __LINE__, &
476 & __FILE__)) RETURN
477 END DO
478 END IF
479 # endif
480
481 # ifdef POSTERIOR_ERROR_I
482 !
483 ! If weak constraint, convolve records 2:Nrec in ADM(ng)%name and
484 ! Write convolved adjoint solution into Hessian NetCDF file for use
485 ! later.
486 !
487 IF (Lposterior.and.(inner.ne.0)) THEN
488 DO ng=1,Ngrids
489 CALL wrt_hessian (ng, Rold(ng), Rold(ng))
490 IF (FoundError(exit_flag, NoError, __LINE__, &
491 & __FILE__)) RETURN
492 END DO
493 END IF
494 # endif
"
In line 472, shouldn't it be "IF (model.eq.iTLM)" while using tl_wrt_ini in line 474? Or any specific reason for the current setting? I am not 100% sure about whether this causes the time record mismatch in the Hessian file. But the time record mismatch in the Hessian file does occur.