Dynamic vs static error estimation in statistical optimization

[Stig Syndergaard]

Investigations of new WO (GS2 test data) and GO (NRT data) processsings, revealed several cases where the dynamic error estimation puts the transition between observation and background either too high or too low.

When it is too high, even moderate ionospheric noise can result in negative optimized bending angle near the top of profiles, which could produce negative refractivity and very noisy dry temperature at high altitudes, if it doens't totally fail to process to dry temperature.

When it is too low, it can result in overall biases and removal of stratospheric features, such as gravity waves.

Part of the problem might be the logic of the determination of the interval above 50 km where the observation noise is estimated.

As it is:

     ALLOCATE(m_diff(nh))
     m_diff = ABS( ba_diff(2,:) - ba_diff(1,:) - ba_diff0) > 6.0_wp*ba_diffS
     IF (ANY(m_diff(iil:nh) .AND.   &
            impact_LH(iil:nh)-config%r_curve > 75000.0_wp)) THEN
       iiu = SUM(MaxLoc(impact_LH(:), impact_LH(:)-config%r_curve < 70000.0_wp))
     ELSE IF (ANY(m_diff(iil:nh) .AND. &
            impact_LH(iil:nh)-config%r_curve > 70000.0_wp)) THEN
     ELSE IF (ANY(m_diff(iil:nh) .AND. &
            impact_LH(iil:nh)-config%r_curve > 65000.0_wp)) THEN
       iiu = SUM(MaxLoc(impact_LH(:), impact_LH(:)-config%r_curve < 60000.0_wp))
     ELSE                                                             
       iiu = SUM(MaxLoc(impact_LH(:), impact_LH(:)-config%r_curve < 80000.0_wp))                                                                                                                               
     ENDIF

This seems odd, because if there are, e.g., scintillations above 75 km, the cut is made at 70 km, but there is then no check on if there are also scintillations below 70 km.

A better logic could be:

     ALLOCATE(m_diff(nh))
     m_diff = ABS( ba_diff(2,:) - ba_diff(1,:) - ba_diff0) > 6.0_wp*ba_diffS
     IF (ANY(m_diff(iil:nh) .AND.   &
            impact_LH(iil:nh)-config%r_curve < 66000.0_wp)) THEN
       iiu = SUM(MaxLoc(impact_LH(:), impact_LH(:)-config%r_curve < 60000.0_wp))
     ELSE IF (ANY(m_diff(iil:nh) .AND. &
            impact_LH(iil:nh)-config%r_curve < 71000.0_wp)) THEN
       iiu = SUM(MaxLoc(impact_LH(:), impact_LH(:)-config%r_curve < 65000.0_wp))
     ELSE IF (ANY(m_diff(iil:nh) .AND. & 
            impact_LH(iil:nh)-config%r_curve < 76000.0_wp)) THEN
       iiu = SUM(MaxLoc(impact_LH(:), impact_LH(:)-config%r_curve < 70000.0_wp))  
     ELSE                                                                                                                                                                                                      
       iiu = SUM(MaxLoc(impact_LH(:), impact_LH(:)-config%r_curve < 80000.0_wp))                                                                                                                               
     ENDIF

This seems to give better results, but more investigation is needed.

Alternatively, one could totally skip the dynamic error estimation such that the transition is about the same for all occultations. This would mean that noisy profiles are noisy, also in refractivity and dry temperature, but at least that would reflect the data, and not as now, where sometimes noisy profiles are smoothed too much and less noisy ones end up with excessive noise.

The attached files show some results:

XXXX: processed with GPAC 0.4.1

ROPP: processed with ROPP

PFIX: processed with ROPP after fixing a bug at top sample

HACK: setting a priori ionospheric bending very high (do not smooth L1 - L2 bending)

SAME: setting fixed error estimates on both observation and background giving about the same transition height for all profiles (relative error in background means slightly different absolute errors because all profiles are a bit different, resulting in slightly different transition zones).

More investigation is needed.