Re: Fwd: wrf_map() doesn't put terrain in right place

I was searching through some old WRF e-mails and came across this one
from last year. This error seems to correspond exactly with the
troubles we are seeing, a westward shift in the data fields and terrain
by about 1 dx. This problem seems to be in the raw WRF output which
then gets into grib output as well.

++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Hi all,

I just stumbled on a bug that with certain grid sizes results in
incorrect coordinates to be stored in output grib gds. The centre
longitude is shifted one grid point to east or west, resulting in all
data being shifted in all model grib output. Or more specifically, the
resulting values are contradictory with the grid corner coordinates,
results may vary depending on the reading program.

This occurs with current code if (JM+1)/2 = e_sn/2 is even.

Attached is a patch to correct this problem, against
wrfpostproc_v2.2.1.tar.gz files.

Another bug in the code is, that if the centre point of the domain is
close to the poles, and the centre point is not a mass point, the
interpolated latitude coordinate of the centre point is incorrect.

The best solution would be to carry the centre point from global
attributes in the netcdf file, CEN_LAT and CEN_LON. I have not yet
checked how to do that or if it will work. Maybe later. And I do not
know how this relates to nested domains.

Regards,
Jaakko Hyvätti

------------------------------------------------------------------------

diff -ru old/wrfpostprocV2/sorc/wrfpost/INITPOST_NMM.f wrfpostprocV2/sorc/wrfpost/INITPOST_NMM.f
--- old/wrfpostprocV2/sorc/wrfpost/INITPOST_NMM.f 2007-02-21 05:07:49.000000000 +0200
+++ wrfpostprocV2/sorc/wrfpost/INITPOST_NMM.f 2007-06-27 18:26:15.000000000 +0300
@@ -1740,17 +1740,34 @@
 ! pos east
        call collect_loc(gdlat,dummy)
        if(me.eq.0)then
         latstart=nint(dummy(1,1)*1000.)
         latlast=nint(dummy(im,jm)*1000.)
 ! temporary patch for nmm wrf for moving nest. gopal's doing
+ jcen = (jm+1)/2
         if(mod(im,2).ne.0) then
- icen=(im+1)/2
- jcen=(jm+1)/2
- cenlat=nint(dummy(icen,jcen)*1000.)
+ if (mod((jm+1)/2,2).ne.0) then
+ icen=(im+1)/2
+ cenlat=nint(dummy(icen,jcen)*1000.)
+ print*,'case 1',im,jm,icen,jcen,cenlat
+ else
+ icen=(im-1)/2
+ cenlat=nint(0.5*(dummy(icen,jcen)+dummy(icen+1,jcen))*1000.)
+ print*,'WARNING! GRIB CENTRE LATITUDE ONLY APPROXIMATE'
+ print*,'WARNING! ERROR INCREASES NEAR POLES'
+ print*,'case 2',im,jm,icen,jcen,cenlat
+ end if
         else
- icen=im/2
- jcen=(jm+1)/2
- cenlat=nint(0.5*(dummy(icen,jcen)+dummy(icen+1,jcen))*1000.)
+ if (mod((jm+1)/2,2).ne.0) then
+ icen=im/2
+ cenlat=nint(0.5*(dummy(icen,jcen)+dummy(icen+1,jcen))*1000.)
+ print*,'WARNING! GRIB CENTRE LATITUDE ONLY APPROXIMATE'
+ print*,'WARNING! ERROR INCREASES NEAR POLES'
+ print*,'case 3',im,jm,icen,jcen,cenlat
+ else
+ icen=im/2
+ cenlat=nint(dummy(icen,jcen)*1000.)
+ print*,'case 4',im,jm,icen,jcen,cenlat
+ end if
         end if
        end if
        write(6,*) 'laststart,latlast B calling bcast= '
@@ -1761,17 +1778,30 @@
      +, latstart,latlast
        call collect_loc(gdlon,dummy)
        if(me.eq.0)then
         lonstart=nint(dummy(1,1)*1000.)
         lonlast=nint(dummy(im,jm)*1000.)
 ! temporary patch for nmm wrf for moving nest. gopal's doing
+ jcen = (jm+1)/2
         if(mod(im,2).ne.0) then
- icen=(im+1)/2
- jcen=(jm+1)/2
- cenlon=nint(dummy(icen,jcen)*1000.)
+ if (mod((jm+1)/2,2).ne.0) then
+ icen=(im+1)/2
+ cenlon=nint(dummy(icen,jcen)*1000.)
+ print*,'case 1',im,jm,icen,jcen,cenlon
+ else
+ icen=(im-1)/2
+ cenlon=nint(0.5*(dummy(icen,jcen)+dummy(icen+1,jcen))*1000.)
+ print*,'case 2',im,jm,icen,jcen,cenlon
+ end if
         else
- icen=im/2
- jcen=(jm+1)/2
- cenlon=nint(0.5*(dummy(icen,jcen)+dummy(icen+1,jcen))*1000.)
+ if (mod((jm+1)/2,2).ne.0) then
+ icen=im/2
+ cenlon=nint(0.5*(dummy(icen,jcen)+dummy(icen+1,jcen))*1000.)
+ print*,'case 3',im,jm,icen,jcen,cenlon
+ else
+ icen=im/2
+ cenlon=nint(dummy(icen,jcen)*1000.)
+ print*,'case 4',im,jm,icen,jcen,cenlon
+ end if
         end if
        end if
        write(6,*)'lonstart,lonlast B calling bcast= '

------------------------------------------------------------------------

Mary Haley wrote:
>
> Hi Don,
>
> I remember this thread, but I'm afraid I don't remember where we left
> it. (Or did we completely drop the ball on it?!)
>
> Anyway, one thing I have seen occasionally is that the WRF data file
> has the wrong center longitude value, and hence the map projection
> and the data don't line up quite right. I'm not sure NCL can do anything
> to detect this.
>
> Can you provide a sample dataset and script for us to look at?
>
> --Mary
>
>
> On Thu, 10 Jul 2008, Don Morton wrote:
>
>> Hi all,
>>
>> I decided to politely re-surface the following thread :) It seems like
>> colleagues are running into this issue more and more, and I'm afraid
>> I don't
>> know enough to say much more than that for some reason the WRF domains
>> aren't plotted correctly by a number of software packages, including
>> the NCL
>> WRF functions (and, we've learned recently, GEMPAK). Is there a way
>> that we
>> out here in the trenches might be wanting to define our original
>> domains to
>> avoid these problems?
>>
>> Thanks,
>>
>> Don
>>
>> ---------- Forwarded message ----------
>> From: Don Morton <donaldjmorton_at_gmail.com>
>> Date: Thu, May 1, 2008 at 10:16 PM
>> Subject: wrf_map() doesn't put terrain in right place
>> To: ncl forum <ncl-talk_at_ucar.edu>, wrfhelp <wrfhelp_at_ucar.edu>
>>
>>
>> Hi, here's another problem I've been running into lately -
>> it's very similar to one I had last summer, that you all helped
>> me find a fix for.
>>
>> I'm trying to use the newer wrf NCL functions, as I really like
>> their potential ease of use. However, in trying to do a simple
>> terrain plot on a WRF nest, the terrain ends up in the wrong
>> place. The script is available at
>>
>> <http://weather.arsc.edu/Miscellaneous/NCL/Maps/wrfHGT-1.ncl>
>>
>> and the input data, "wrfoutd02.nc" is in the same directory.
>> The resulting image is at
>>
>> <http://weather.arsc.edu/Miscellaneous/NCL/Maps/Image1.jpg>
>>
>> You can see that the Alaska Range has moved to the Gulf of Alaska.
>>
>> This is presumably because the WRF output file represents an inner
>> nest of a larger domain, and in this nest the Standard Lon and the
>> Center Lon are different. This "seems" like a common occurrence to me,
>> unless you place the nest smack in the middle of the parent?
>>
>> Last summer, you all helped me fix this issue, and such a
>> script exists at
>>
>> <http://weather.arsc.edu/Miscellaneous/NCL/Maps/wrfHGT-2.ncl>
>>
>> with resulting image at
>>
>> <http://weather.arsc.edu/Miscellaneous/NCL/Maps/Image2.jpg>
>>
>> I tried playing around with some resource settings for the wrf_map()
>> functions, but couldn't get anything to work.
>>
>> Thanks,
>>
>> Don
>> --
>> University of Montana
>> http://www.cs.umt.edu/~morton/ <http://www.cs.umt.edu/%7Emorton/>
>> Arctic Region Supercomputing Center
>> http://weather.arsc.edu/
>>
>>
>>
>> --
>> Arctic Region Supercomputing Center
>> http://weather.arsc.edu/
>> University of Montana
>> http://www.cs.umt.edu/~morton/
>>
Received on Thu Jul 10 2008 - 16:40:42 MDT