Polar wind problem for ERA5 (CDS)

[douglowe]

Hi all,

I'm looking at moving from ERA-Interim to ERA5 input data for WRF. I've downloaded the lat-lon gridded data from the CDS service, and generated successfully(?) met_em* files from this data. However my domain encompasses the north pole, and there is an issue with the ERA5 lat/lon gridded data at the poles - U & V winds are all 0 at 90 N/S: https://confluence.ecmwf.int/pages/viewpage.action?pageId=129134800. The final met_em* files produced from these data are mostly okay, however there is a distinct drop in the U & V wind speeds at the pole (not zero, but a definite reduction).

The suggestion from ECMWF is to replace the U & V wind values at 90 North with those at 89.75 North. However I've never edited grib (1 or 2) files before, and will need to do this for 365 files (at least). Has anyone else done this, or can anyone recommend a good toolset for scripting this process?

Thanks,
Doug

 

[kwerner]

Hi Doug,
Have you ever edited NetCDF files? It may be easier to wait until after you run the metgrid program to edit those files (which are then in NetCDF format). There are several online tools or methods for making modifications to NetCDF files.

 

[douglowe]

Hi kwerner,

I have edited NetCDF files - coding-wise definitely I would find doing this easier than editing GRIB files. Numerically, it would be easier to edit the GRIB inputs before they are run through metgrid, as it seems metgrid.exe (very sensibly) is smoothing out the input data, and so the depression in U & V wind speeds covers several grid cells, not just the polar grid cell. Determining the best way to iron out that depression in the met_em files would be more work than determining how to edit GRIB files I fear.

 

[douglowe]

Alternatively - the CDS system does enable the downloading of the ERA5 data as netCDF, rather than GRIB. Is it possible to use netCDF files as inputs for WPS? If so (and if it's not too much work to setup) then I could do that instead?

 

[kwerner]

Hi,
Yes, you can use netcdf files as input. In this case, you would need to convert the files to intermediate format and then would skip the ungrib step of WPS. There is information found in Chapter 3 of our Users' Guide that can help with converting the files to intermediate format:
https://www2.mmm.ucar.edu/wrf/users/docs/user_guide_v4/v4.1/users_guide_chap3.html#_Writing_Meteorological_Data
But we also have a couple of different scripts available that were given to us by outside groups, specifically for converting netcdf files to intermediate files. You can find those scripts here:
https://www2.mmm.ucar.edu/wrf/users/special_code.html
Just keep in mind that these were created for the use of their particular projects, so some modifications will be necessary, and because they are passed down from groups outside of the WRF Support group, we aren't able to support the scripts.

 

[bartbrashers]

Douglowe,

It might be even easier to write a Fortran program to fix the data in the WPS intermediate files produced by UNGRIB, and run the fixer before running METGRID. The file format of WPS Intermediate Files is well-known, and there's a code chunk in the docs you can use. See https://www2.mmm.ucar.edu/wrf/users...guide_chap3.html#_Writing_Meteorological_Data.

So your workflow would be:
ungrib.exe # creates many FILE:YYYY-MM-DD_HH files
fix_era5_pole_hole.exe FILE*
metgrid.exe

I'm about to use ERA5 for polar regions too, and could benefit from your fix - please post it here, if you choose this method.

 

[douglowe]

hi kwerner & bertbrashers,

Thanks for your advice & links to documentation. Maybe editing the WPS intermediate files will be the easiest solution, if I can pull out information on the location of that data from the intermediate file. I'll have a look at the metgrid source files too - perhaps I can repurpose some of the subroutines from those for this task.

I'll let you know if I get anywhere with this & will post links to any tools that I might manage to make for this.

cheers,
Doug

 

[bartbrashers]

Some code snippets to contribute, from my program intermediate-interp.f90, which interpolates in time between two WPS intermediate files. Useful for data in netCDF format that is not available every interval_seconds. UNGRIB would interpolate in time, but my netCDF-to-intermediate programs don't.

  real, allocatable,dimension(:,:) :: input1         ! input data
  real, allocatable,dimension(:,:) :: output         ! output data

  real :: bad,bad1                    ! missing/bad value flags
  real :: earth_radius                ! radius of the earth in the files
  real :: nlats                       ! Number of latitudes north of equator
  real :: xfcst                       ! Forecast hour of data
  real :: xlvl                        ! Vertical level of data in 2-d array
  real :: startlat, startlon          ! Lat/lon of point in array indicated by startloc
  real :: deltalat, deltalon          ! Grid spacing, degrees
  real :: dx, dy                      ! Grid spacing, km
  real :: xlonc                       ! Standard longitude of projection
  real :: tlat1, tlat2                ! True latitudes of projection

  real :: startlat2, startlon2, deltalat2, deltalon2
  real :: dx2, dy2, xlonc2, tlat12, tlat22, nlats2

  logical :: ok
  logical :: debug = .false.          ! set to .true. to debug
  logical :: is_wind_grid_rel         ! Flag indicating whether winds are
                                      !   relative to source grid (TRUE) or
                                      !   relative to earth (FALSE)
  character (len=8)   :: startloc     ! Which point in array is given by
                                      !   startlat/startlon; set either
                                      !   to 'SWCORNER' or 'CENTER  '
  character (len=8)   :: startloc2
  character (len=9)   :: field        ! Name of the field
  character (len=24)  :: hdate        ! Valid date for data YYYY:MM:DD_HH:00:00
  character (len=25)  :: units        ! Units of data
  character (len=32)  :: map_source   ! Source model / originating center
  character (len=46)  :: desc         ! Short description of data

! Make sure the file exists and open it

  inquire(file=infile1, exist=ok)
  if (.not. ok) then
     write(*,*) "Error: file not found: ",trim(infile1)
     stop
  endif
  open(10,file=infile1,form='unformatted',convert='big_endian', &
       status='old',err=99)

! Open the output file

  open(12,file=outfile,form='unformatted',convert='big_endian',err=99)


!*********************************************************************
! LOOP OVER FIELDS FOUND IN INPUT FILES
!*********************************************************************

  ios = 0 ! a leap of faith

  do while (ios == 0)

! Read the 1st input file

     read(10, iostat=ios) version

     if (ios /= 0) exit ! should be triggered by end-of-file

     if (version /= 5) then
        write(*,*) "This code only works for Intermediate files version 5 (WPS)."
        write(*,*) trim(infile1)," is version ",version
        stop
     endif
     
     read(10) hdate,xfcst,map_source,field,units,desc,xlvl,xdim,ydim,iproj

     if (iproj .eq. 0) then        ! Cylindrical equidistant
        read(10) startloc,startlat,startlon,deltalat,deltalon,earth_radius
     else if (iproj .eq. 1) then   ! Mercator
        read(10) startloc,startlat,startlon,dx,dy,tlat1,earth_radius
     else if (iproj .eq. 3) then   ! Lambert conformal
        read(10) startloc,startlat,startlon,dx,dy,xlonc,tlat1,tlat2,earth_radius
     else if (iproj .eq. 4) then   ! Gaussian
        read(10) startloc,startlat,startlon,nlats,deltalon,earth_radius
     else if (iproj .eq. 5) then   ! Polar stereographic
        read(10) startloc,startlat,startlon,dx,dy,xlonc,tlat1,earth_radius
     end if

     read(10) is_wind_grid_rel     ! rotation flag

     if (debug) then
        write(*,*) "Read infile:"
        write(*,*) "  hdate = ",hdate
        write(*,*) "  field = ",field
        write(*,*) "  units = ",units
        write(*,*) "  desc  = ",desc
        write(*,*) "  dims  = ",xdim,ydim
        write(*,*) "  iproj = ",iproj
     endif

     if (.not. allocated(input1)) allocate( input1(xdim,ydim) )
     read(10) input1               ! array of data

### INSERT FIXER CODE HERE ###

! Write the output file

     write(12) version
     write(12) hdate,xfcst,map_source,field,units,desc,xlvl,xdim,ydim,iproj

     if (iproj .eq. 0) then ! Cylindrical equidistant
        write(12) startloc,startlat,startlon,deltalat,deltalon,earth_radius
     else if (iproj .eq. 1) then ! Mercator
        write(12) startloc,startlat,startlon,dx,dy,tlat1,earth_radius
     else if (iproj .eq. 3) then ! Lambert conformal
        write(12) startloc,startlat,startlon,dx,dy,xlonc,tlat1,tlat2,earth_radius
     else if (iproj .eq. 4) then ! Gaussian
        write(12) startloc,startlat,startlon,nlats,deltalon,earth_radius
     else if (iproj .eq. 5) then ! Polar stereographic
        write(12) startloc,startlat,startlon,dx,dy,xlonc,tlat1,earth_radius
     end if

     write(12) is_wind_grid_rel   !  3) WRITE WIND ROTATION FLAG

     if (debug) then
        write(*,*) "Writing to outfile:"
        write(*,*) "  hdate = ",hdate
        write(*,*) "  field = ",field
        write(*,*) "  units = ",units
        write(*,*) "  desc  = ",desc
        write(*,*) "  dims  = ",xdim,ydim
        write(*,*) "  iproj = ",iproj
     endif

     write(12) output             !  4) WRITE 2-D ARRAY OF DATA

  end do ! main loop over blocks in the input files

  close(10) ! infile1
  close(12) ! outfile

  deallocate( input1, input2, output )