Indices in ropp_fm_hybrid_phi_flv.f90

[frcm]

In the routine for calculating geopotential height on hybrid vertical levels (ropp_fm_hybrid_phi_flv.f90), the following construction

! Calculate the integral (eq. 2.21)
! ---------------------------------

  do lvl = 1, size(p_hlv)
     phi_hlv(lvl) = phi_sfc - sum(del_phi(lvl:size(Tv_flv)))
  enddo

exists. however, the half levels (currently) have 61 elements, but there are only 60 full levels. Thus, for the calculation of phi_hlv(lvl), the result is

phi_hlv(61) = phi_sfc - sum(del_phi(61:60))

where del_phi(61) isn't defined. Is that a bug in the implementation, or does it use a particular feature of Fortran 90?

Chris.

[Huw Lewis]

To be fixed in new calculation of geopotential height ropp_fm_ecmwf_geop, which assumes that array indices increase away from the surface (i.e. index 1 towards surface). In this case, the code reads

! 7.5 Calculate geopotential height integral
 
  do lvl = 1, n_hlv
      geop_hlv(lvl) = x%geop_sfc + sum(del_geop(1:max(lvl-1,1)))
  enddo

At the lowest level, lvl=1, the sum is evaluated as

   geop_hlv(1) = x%geop_sfc + sum(del_geop(1:1))

source:ropp_src/branches/dev/Share/r1182_frhl/ropp_fm/model_ecmwf/ropp_fm_ecmwf_geop.f90@1221

[Huw Lewis]

Solution in ropp_fm_ecmwf_geop.f90 is not appropriate, leading to numerical errors relative to solutions obtained using v1.0 software. Require geop_hlv(1) = x%geop_sfc (with sum(X(1:0))=0). Possible to use ECMWF implementation (see commented section in code) to avoid this ambiguity.

[Huw Lewis]

ROPP v1.2 code (model_ecmwf/ropp_fm_state2state_ecmwf) follows same algorithm as v1.0 (accounting for updated array ordering convention - increasing away from surface).

This reads:

do lvl = 1, n_hlv
      geop_hlv(lvl) = x%geop_sfc + sum(del_geop(1:lvl-1))
  enddo

Note, this involves the sum(del_geop(1:0)) for lvl=0. Testing shows this sum is equal to zero (i.e. geop_hlv(1) = x%geop_sfc), rather than being a bug. Ticket closed.