New routines to initialize the Obs1dRefrac and Obs1dBangle structures [ ?]

[Ian Culverwell]

https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_fm/build/Makefile.am https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_fm/common/ropp_fm.f90 https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_fm/refrac_1d/Makefile.am https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_fm/refrac_1d/ropp_fm_tdry.f90 https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_fm/refrac_1d/ropp_fm_tdry_1d.f90 https://trac.romsaf.org/ropp/changeset/5427/ropp_src/branches/dev/Share/dmi_trunk_9.0/ropp_fm/tools/ropp_fm_bg2ro_1d.f90

Original changesets:

New routines to initialize the Obs1dRefrac and Obs1dBangle structures: https://trac.romsaf.org/ropp/changeset/4777/ropp_src/branches/dev/Share/dmi_trunk_8.1/ropp_fm/tools/ropp_fm_bg2ro_1d.f90

[Ian Culverwell]

Made a preliminary stab at this in local branch. (Difficult because of the mismatch between the ROPP9.1 code and the ROPP8.1 DMI source.) It involves the introduction of new subroutines set_obs_levels_refrac_file and set_obs_levels_bangle_file in ropp_fm_bg2ro_1d.f90, which copy data from levels files to the output. These are basically the same as the routines in ropp_fm_roprof2obs.f90, except that refrac and bangle are nullified, and the covariances aren't copied.

Also rename set_obs_levels_refrac to set_obs_levels_refrac_uni, and set_obs_levels_bangle2 to set_obs_levels_bangle_uni. These define the uniform levels sets (and 247L for bangle).

Run through the attached test.sh program, which includes a -l test. Results are ~OK, apart from RoC and refrac_sigma, which are different.

Input levels file: WOP_Study_Cases_prof3_cntl1.nc

ncks -Q -H -ddim_lev2a,0,4 -vr_coc,refrac_sigma  WOP_Study_Cases_prof3_cntl1.nc
r_coc[0]=10.0354 
r_coc[1]=64.6754 
r_coc[2]=41989.8 

refrac_sigma[0]=-99999000 
refrac_sigma[1]=-99999000 
refrac_sigma[2]=-99999000 
refrac_sigma[3]=-99999000 
refrac_sigma[4]=-99999000 

Output file: WOP_Study_Cases_prof3_test6.nc

ncks -Q -H -ddim_lev2a,0,4 -vr_coc,refrac_sigma WOP_Study_Cases_prof3_test6.nc 
r_coc[0]=10035.4 
r_coc[1]=-9.9999e+07 
r_coc[2]=-9.9999e+07 

refrac_sigma[0]=0.0173205080757 
refrac_sigma[1]=0.0173205080757 
refrac_sigma[2]=0.0173205080757 
refrac_sigma[3]=0.0173205080757 
refrac_sigma[4]=0.0173205080757 

The fix was supposed to fix the 0.01732 problem, which it obviously hasn't. It has also apparently generated a problem with r_coc (multiplied by 1000, apparently).

To be continued.

[Ian Culverwell]

test.sh

[Ian Culverwell]

(In the input file we have

ncks -Q -H -ddim_lev2a,0,4 -vr_coc WOP_Study_Cases_prof3.nc            
r_coc[0]=10.0354 
r_coc[1]=64.6754 
r_coc[2]=41989.8

and no refrac_sigmas.)

[Ian Culverwell]

(And it has a unit of m:

	float r_coc(dim_unlim, xyz) ;
		r_coc:long_name = "Centre of curvature for the reference coordinate" ;
		r_coc:units = "metres" ;
		r_coc:valid_range = -50000., 50000. ;
		r_coc:reference_frame = "ECF" ;

)

[Ian Culverwell]

The reason is that there is a CALL ropp_fm_set_units(ro_data) statement in ropp_fm_bg2ro_1d.f90, which defines the units of r_coc to be km:

  CALL ropp_fm_set_units_range(rodata%GEOref%units%r_coc, "kilometres",  &
                               rodata%GEOref%range%r_coc)

(Because km is considered an appropriate unit for r_coc within the FM, presumably.)

But the current code doesn't do the same to the l_data structure read from the file, which means that we read in r_cocs of ~40 000 m. That would be fine, except that ropp_io_write also contains a call to ropp_fm_set_units(ro_data), which means that these r_cocs are interpreted as ~ 40 000 km. They therefore fail the implicit range checking on output, which is why we end up with

ncks -H -Q -vr_coc WOP_Study_Cases_prof3_test6.nc               
r_coc[0]=10035.4 
r_coc[1]=-9.9999e+07 
r_coc[2]=-9.9999e+07

The solution is to add a line

          CALL ropp_fm_set_units(l_data) ! Must use the same units as ro_data

in the appropriate part of the code (sec 9.2). This delivers

ncks -H -Q -vr_coc WOP_Study_Cases_prof3_test6.nc
r_coc[0]=10.0354 
r_coc[1]=64.6754 
r_coc[2]=41989.8 

as required.

[Ian Culverwell]

I cannot reproduce the earlier problems with the odd refrac_sigmas: I now get

ncks -H -Q -ddim_lev2a,0,2 -vrefrac_sigma WOP_Study_Cases_prof3_test6.nc
refrac_sigma[0]=-99999000 
refrac_sigma[1]=-99999000 
refrac_sigma[2]=-99999000 

[Ian Culverwell]

I do recover the odd values of I revert to calling ropp_fm_roprof2obs rather than set_obs_levels_refrac_file. The reason for this is that when the obs1drefac structure is being populated in rropp_fm_roprof2obs.f90, it passes through the following code:

  IF (ASSOCIATED(y%cov%d)) DEALLOCATE(y%cov%d)
  CALL callocate(y%cov%d, (n*(n+1)/2))

  DO i = 1, n
    IF (ro_data%Lev2a%refrac(i) > 0.0_wp) THEN
      IF (ro_data%Lev2a%refrac_sigma(i) > 0.0_wp) THEN
       ! matrix_pp type, uplo = 'U'
        y%cov%d(i + i*(i-1)/2) = ro_data%Lev2a%refrac_sigma(i)**2
      ELSE
        y%cov_ok = .FALSE.
      END IF
    ELSE 
      y%cov%d(i + i*(i-1)/2) = 0.0003_wp
    ENDIF
  END DO

The default refractivity covariance appears to be a straight copy of the default bending angle covariance of 3e-4 rad²:

    DO i = 1, n
      IF (ro_data%Lev1b%bangle(i) > 0.0_wp) THEN
        IF (ro_data%Lev1b%bangle_sigma(i) > 0.0_wp) THEN
         ! matrix_pp type, uplo = 'U'
          y%cov%d(i + i*(i-1)/2) = ro_data%Lev1b%bangle_sigma(i)**2
        ELSE
          y%cov_ok = .FALSE.
        END IF
      ELSE
        y%cov%d(i + i*(i-1)/2) = 0.0003_wp  ! = 1 deg**2
      ENDIF
    END DO

  ELSE

So it's probably a mistake. I don't know what a good default refrac might be - ask DMI for ideas. But anyway, when the refractivity covariance of 3e-4 N-units² is square-rooted in ropp_fm_obs2roprof.f90, we get a refrac_sigma of 1.732e-2 N-units.

[Ian Culverwell]

But why isn't it an issue for bangle? Because its covs aren't associated, for some reason:

 12.1a ASSOCIATED(obs_refrac%cov%d) =  T
 13.1a ASSOCIATED(obs_bangle%cov%d) =  F
 in ropp_fm_obs1drefrac2roprof: ASSOCIATED(y%cov%d) =  T
 14a ASSOCIATED(obs_refrac%cov%d) =  T
 in ropp_fm_obs1dbangle2roprof: ASSOCIATED(y%cov%d) =  F
 14b ASSOCIATED(obs_bangle%cov%d) =  F

Why?

[Ian Culverwell]

Because I was using the new code for bangle, so forget that last bit. But in fact the problem remains when I use ropp_fm_roprof2obs for both refrac and bangle. This causes bangle_sigma to be set to 1.732e-2, as for refrac_sigma, but this is greater than the max allowed value of 1e-2, so it gets set to ropp_MDFV by the automatic range checking of the output. (Max refrac_sigma = 50, so it's not an issue for that.)

[Ian Culverwell]

While doing this, it became clear that the refrac and bangle portions of roprof2obs (attached) were slightly inconsistent. In particular, 'bangle' returns if the number of lev1b levels is zero, while 'refrac' carries on if n_lev2a = 0. The latter seems a bit weird, so update ropp_fm_roprof2obs.f90 to make refrac do the same as bangle. This, and some cosmetic changes to this routine and ropp_fm_obs2roprof.f90 have been made at r5927.

[Ian Culverwell]

ropp_fm_roprof2obs.f90_refrac

[Ian Culverwell]

ropp_fm_roprof2obs.f90_bangle

[Ian Culverwell]

The required changes to ropp_fm_bg2ro_1d.90 have been committed at r5928. These comprise:

• a line CALL ropp_fm_set_units(l_data), to fix the r_coc problem above;
• renaming set_obs_levels_bangle/refrac to set_obs_levels_bangle/refrac_uni
• the introduction of new routines set_obs_levels_bangle/refrac_file to take only what's needed from the 'levels' file.

Checks out OK in test.sh and standard 'make tests'

Leaving ticket open, as there may be something to discuss with DMI.

[Ian Culverwell]

(Note that, for some reason, even if ro_data%lev2a%npoints = 0, then you can still print*,'9.1 ro_data%lev2a%refrac_sigma = ', ro_data%lev2a%refrac_sigma etc on the first profile (you just get 9.1 ro_data%lev2a%refrac_sigma =), but that this print* statement gives a seg fault on later profiles (i.e. after call ropp_io_free(ro_data) has been called on the first profile). Something to bear in mind and live with, rather than try to fix, I think.)

[Ian Culverwell]

That's more than enough on this. Closing ticket.