fix for zero total ozone pressure in EnKF (issue #625) (#626)

The fix involves setting the pressure to 0.001Pa for ozone obs that have
zero pressure (to avoid Inf when log(p) calculated), and turning off
ob-space vertical localization. Has no effect on current operational
setup which uses model-space vertical localization (modelspace_vloc=T).

---------

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

src/enkf/enkf_obsmod.f90

@@ -262,7 +262,6 @@ subroutine readobs()
 allocate(corrlengthsq(nobstot),lnsigl(nobstot),obtimel(nobstot))
 lnsigl=1.e10
 do nob=1,nobstot
-   oblnp(nob) = -log(obpress(nob)) ! distance measured in log(p) units
    if (obloclon(nob) < zero) obloclon(nob) = obloclon(nob) + 360._r_single
    radlon=deg2rad*obloclon(nob)
    radlat=deg2rad*obloclat(nob)
@@ -283,6 +282,13 @@ subroutine readobs()
       lnsigl(nob)=latval(deglat,lnsigcutoffnh,lnsigcutofftr,lnsigcutoffsh)
    end if
    endif
+   ! total column ozone has pressure set to zero, set to 0.001Pa
+   ! and turn vertical localization off (no effect if modelspace_vloc=T)
+   if (obpress(nob) < 0.001 .and. obtype(nob)(1:3) .eq. ' oz') then
+      lnsigl(nob) = 1.e30    ! turn ob-space vert localization off
+      obpress(nob) = 0.001   ! set to a non-zero value
+   endif
+   oblnp(nob) = -log(obpress(nob)) ! distance measured in log(p) units
    corrlengthsq(nob)=latval(deglat,corrlengthnh,corrlengthtr,corrlengthsh)**2
    if ( (obtype(nob)(1:3) == 'dbz' .or. obtype(nob)(1:3) == ' rw') .and. l_use_enkf_directZDA ) then
        corrlengthsq(nob)=latval(deglat,corrlengthrdrnh,corrlengthrdrtr,corrlengthrdrsh)**2