Create GSI covariance model in SABER for regional JEDI analysis

[Masanori-NOAA]

Although there is GSI covariance model in SABER for global analysis, there is no one for regional analysis. Therefore, we want to introduce GSI covariance model in SABER for regional analysis.

[Masanori-NOAA]

It is planned to create the regional GSI covariance model by compiling the full GSI code and link it to SABER interface. This is because the codes of GSI B model for global analysis (GSIbec) deleted the parts related to regional analysis and it is too complicated to rewrite it for use in regional case.

It has been confirmed that full GSI code can be linked to SABER interface on Orion after recompiling some NCEPLIBS (sigio, sfcio, nemsio, ncio, ncdiag, wrf_io) with “-fPIC” option by myself because of the following error.

“relocation R_X86_64_32 against `__STRLITPACK_14.0.2' can not be used when making a shared object; recompile with –fPIC”

　

[Masanori-NOAA]

It ran to the end but the minimization process failed. The first calculated value of the residual norm is “NaN” as follows. This occurs when including subroutines to convert variables from physical space to control space (control2state_ad, control2state).

 0: Initial RHS squared norm = 0.00242799
 0: normr0 0.00242799
 0:   Residual norm ( 0) = -nan
 0:   Norm reduction ( 0) = 1
 0:
 0:   Quadratic cost function: J   ( 0) = 0.09641414870677294
 0:   Quadratic cost function: Jb  ( 0) = 0
 0:   Quadratic cost function: JoJc( 0) = 0.09641414870677294
 0:
 0:
 0:  DRPCG Starting Iteration 1
 0: rdots nan      iteration    0
 0: DRPCG end of iteration 1
 0:   Residual norm ( 1) = nan
 0:   Norm reduction ( 1) = nan
 0:
 0:   Quadratic cost function: J   ( 1) = -nan
 0:   Quadratic cost function: Jb  ( 1) = -nan
 0:   Quadratic cost function: JoJc( 1) = -nan

[Masanori-NOAA]

When commented out the variable “pres” from the state_vector in the namelist (anavinfo), “NaN” no longer appeared (because there is no “pres” in the atlas_field). I calculated the iterations five times, including the multiplication of the B matrix, but the cost dropped too rapidly compared to the case using BUMP. The interpolation process from the JEDI model grid to the GSI analysis grid may not be working well now.

Cost (case using GSI B)

| | J | JoJc | | - | - | - | | 0 | 0.1035626576144349 | 0.1035626576144349 | | 1 | 0.0001853076338450693 | 3.315762645833954e-07 | | 2 | 0.0001853076338450693 | 3.315762645833683e-07 | | 3 | 0.0001853076338450554 | 3.315762645694363e-07 | | 4 | 0.0001853076338450416 | 3.315762645555856e-07 | | 5 | 0.0001853076338450416 | 3.315762645555856e-07 |

Cost (case using BUMP)

| | J | JoJc | | - | - | - | | 0 | 0.09641414870677294 | 0.09641414870677294 | | 1 | 0.09527945473022176 | 0.09415811491836208 | | 2 | 0.09527945473022176 | 0.09415811491836208 | | 3 | 0.09527945473022176 | 0.09415811491836208 | | 4 | 0.09527945473022176 | 0.09415811491836208 | | 5 | 0.09527945473022176 | 0.09415811491836208 |

[Masanori-NOAA]

There was a mistake in calculating the units of lat/lon of the gsi analysis grid used for interpolation calculations in radians instead of degrees. This part was corrected and the analysis grid was made to cover the JEDI’s model grid. With this setting, a single-obs test of surface temperature was conducted but still see that the analysis result got too close to the observation value.

Cost

| | J | JoJc | | - | - | - | | 0 | 0.1035626576144349 | 0.1035626576144349 | | 1 | 0.0003379927161903523 | 1.103091392495859e-06 | | 2 | 0.0003379927161902829 | 1.10309139242609e-06 |

[Masanori-NOAA]

Decided to first check whether it works properly with the dirac test of air temperature, based on the advice of Shun and Ting. After modifying interface part related to retrieving “Tv” from “T” of the atlas filed, the dirac test of fv3jedi case result became as follows.

When bypassing the central block, the peak value was 3.8 ( This should be less than 1). When bypassing both the central and outer (interpolation) block, the peak value became 1. This means that I still need to check whether the interpolation part is working properly.

[Masanori-NOAA]

After the modification in the previous comment, the problem in which the analysis result got too close to the observation value in the single-obs test was solved.

Cost

| | J | JoJc | | - | - | - | | 0 | 0.1035626576144349 | 0.1035626576144349 | | 1 | 0.1004734898498509 | 0.09747646878850486 | | 2 | 0.1004734898498508 | 0.09747646878850484 |

[Masanori-NOAA]

Thanks to a comment from Ting, it turned out that when running a Dirac test, the cause of the incorrect interpolation value using saber::gsi::interpolater in the dirac test was because of a bug in the handling of atlas' halo exchange (JCSDA-internal/saber#767). The fv3-jedi interface has already been updated to solve this (JCSDA-internal/fv3-jedi/pull/1170) but fv3jedi’s version I used in the RDASApp (https://github.com/NOAA-EMC/RDASApp/tree/f2548724bd5f38d84eabd085b372946c1f1f3c89) was older than this revision. When I use the old version of fv3jedi, haloExchange created additional dirac impulses in the halo when the location of the impulse was over halo region(x=198,y=116,lay=50), resulting in unnaturally large interpolated value. Yesterday, the versions of the submodules of RDASApp was updated (#121) and the correction of fv3jedi's handling of halo exchange was incorporated. When I use this RDASApp’s version as a base, the problem with interpolation calculation was no longer observed.

When bypassing the multiplication of the B matrix (namely only doing ajoint and forward interpolation), the previous maximum value was 3.8, but in the re-experiment, the value is 0.9.

[Masanori-NOAA]

Compared the single-obs test results of GSI and the recursive filter linked to saber/jedi in case of fv3jedi’s surface airTemp assimilation of the RDASApp. As for GSI’s obs data, created prepbufr based on the wiki page (https://github.com/NOAA-EMC/RDASApp/wiki/Running-FV3%E2%80%90JEDI-Hofx-Validation). There is small difference in the maximum value of the increment.

[Masanori-NOAA]

Conducted the experiment again after matching the hofx values by feeding GSI-geovals into vertInterp in RDASApp, thanks to comments from Donnie. In addtion, changed the error inflation factor so that the final oberr values become same.