[Wrf-users] Bug in wrfpostproc 2.2.1 NMM centre lon/lat code

[Jaakko Hyvätti]

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
-------------- next part --------------
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= '