vertransform.f90 occasionally fails to initialize

[adingwell]

When running with smaller domains verttransform might fail to initialize ixm and jym.

At initialization the following block is run in verttransform.f90:

  ! Search for a point with high surface pressure (i.e. not above significant topography)
  ! Then, use this point to construct a reference z profile, to be used at all times
  !**************************************************************************************

    do jy=0,nymin1
      do ix=0,nxmin1
        if (ps(ix,jy,1,n).gt.100000.) then
          ixm=ix
          jym=jy
          goto 3
        endif
      end do
    end do
3   continue

In order to work, there has to be a grid point within the domain with a surface pressure above 1000 hPa. Such a point might not be present in smaller domains, where there is only a small number of grid points over sea surface.

It would be good if this block at least produces an error if no reference point is found.

[pesei]

I suggest that we take the lowest grid point found if there is none >1000 hPa.
I'll take the ticket for v9.2.

[adingwell]

I just noticed that this bug found its way into 10.2.

I've added the following lines between the loop and the continue statement in my local version:

ij_min  = MINLOC(ps(:,:,1,n))  ! ij-index of minimum pressure
ixm = ij_min(1)
jym = ij_min(2)

but I haven't tested it for a problematic domain yet (as my problems this time were related to #197)

[pesei]

Replying to adingwell:

I just noticed that this bug found its way into 10.2.

I've added the following lines between the loop and the continue statement in my local version:

ij_min  = MINLOC(ps(:,:,1,n))  ! ij-index of minimum pressure

don't use min it is the opposite of what we do in the code! We have to use maxloc to find the point the ca lowest elevation. Also, we have to be aware that this will cause results to differ from precvous versions, so maybe either offer option for old one or do this only if the old condition fails.

[pesei]

This is now fixed in dev, see
verttransform_ecmwf: changeset:1a8fbee08b4dd8d4aab8d4f78b6bd0bb37e21f13/flexpart.git
verttransform_gfs: changeset:f251e576e36228966cea4e1a90c93e625f5cd2d9/flexpart.git

Currently, there is a parameter in subroutine header which can be set to 0 for the old selection of the reference grid cell, or 1 for the new method. This might be moved to par_mod.f90 later. The method selects the grid point whose sfc pressure is closest to the mean over the domain (skipping regions with pressure lower than first-guess mean plus 1 sigma).

For a global GFS data set, this is a significant difference. The old method takes a very cold Antarctic grid cell and has about 10% lower heights.

[pesei]

The fix was implemented in branch univie and it was not taken into a release. Therefore I reopen the ticket.