id	summary	reporter	owner	description	type	status	priority	milestone	component	version	resolution	keywords	cc
461	ropp_fm_bg2ro_1d -new_op can generate NaNs above model top	Ian Culverwell	Ian Culverwell cburrows sti	"Hans Gleisner (DMI) reports that the new interpolation scheme can generate !NaNs above the model top:
{{{
          i         geop(i)                  refrac(i)
         773   77200.0000000000       4.595168840727702E-003
         774   77300.0000000000       4.513892752756471E-003
         775   77400.0000000000       4.433896771776462E-003
         776   77500.0000000000      NaN                    
         777   77600.0000000000      NaN                    
         778   77700.0000000000      NaN                    

}}}

This is not ideal.

I have traced this to the section of '''ropp_fm/refrac_1d/ropp_fm_refrac_1d_new.f90''' which says
{{{
! 3.4 Find model layer (j to j+1) in which ob lies
! ------------------------------------------------
    DO j=1, x%n_lev-1
      IF (y%geop(i) < x%geop(1)) THEN
        below_lowest_lev=.TRUE.
        EXIT
      END IF
      IF (y%geop(i) > x%geop(j) .AND. y%geop(i) < x%geop(j+1)) THEN
        EXIT
      END IF
    END DO
}}}

The problem is that if y%geop(i) > x%geop(x%n_lev), ie, the ob is above the model top, then the last IF ... EXIT is never executed and j is updated to the 'next' index in the loop, ie j becomes x%n_lev.  All later references to variables at index j+1 are therefore outside the array bounds.  This causes, in Hans' case, alpha (=-dlog q/dz) to be minus infinity, and generates all the !NaNs.

I see two possible solutions.

(1) Extrapolate the top layer variables, between j=n_lev-1 and n_lev
above the model top.  In other words, set j=x%n_lev-1 if the ob is above the model top.  This is fine unless one is extrapolating wildly, as Hans was (up to 750 km in his case).  You can then get negative temperature t_ob etc.  This causes !NaNs for a different reason.

(2) Safer would be to introduce a new logical {{{above_highest_lev}}}, analagous to {{{below_lowest_lev}}}, which would set refracs above the model top to ropp_MDFV.  This fixes Hans's problem, but is it what we want to do?

(Attached is Hans' test file, which he ran with
ropp_fm_bg2ro_1d  singlefile.bg.nc  -f  -comp  -new_op  --no-ranchk  -d -o singlefile.bg_out.nc.)

"	defect	closed	normal	9.0	ROPP (all)	8.0	fixed