GRIB-format

På SMHI (och andra meteorologiska institut) sparas modelldata i något som kallas GRIB-format och är definierat av WMO. Formatet har nyligen reviderats. Den nya versionen kallas GRIB2 och den gamla GRIB1. Historiska analysdata som SMHI levererar är i GRIB1-format.

GRIB är ett binärt dataformat för data på ett beräkningsbart gridnät. Filen innehåller mestadels kod som kombineras med externa metadata för att rätt kunna läsa ut och tolka datasetet i GRIB-filen. Här finner du en översiktlig beskrivning av GRIB-formatet.

Nedladdningsbara meteorologiska modeller har ett roterat latitud/longitud gridnät medan oceanografiska modeller har ett reguljärt latitud/longitud gridnät.

För att kunna läsa data på GRIB-format behöver du:

  • Program för att öppna och avkoda GRIB-fil
  • Översätta möjliga roterade koordinater
  • Rotera möjliga roterade vindar
  • Kodtabell för att tolka parameternummer

Program för läsande av GRIB-fil

För att kunna läsa GRIB-packad data behövs program för att packa upp GRIB-filen och tolka data. Vid SMHI används mestadels program framtagna av ECMWF (http://www.ecmwf.int).

Dokumentation av den senaste varianten av programpaket (GRIB-API) samt möjlighet att ladda ned källkoden finns på ECMWF:s hemsida för GRIB-API,
https://software.ecmwf.int/wiki/display/GRIB/Home.

Andra program för läsning, skrivning och visualisering av data på GRIB-format finns att tillgå.

Roterade koordinater

De meteorologiska modeller som kan hämtas som ”Öppna data” har data lagrat på ett roterat gridnät. Roterade koordinater innebär att man flyttar sydpolen från lat/lon -90/0. Skälet till rotationen är att man genom att flytta ekvatorn till Sveriges latitud får mindre skillnad i avståndet mellan gridpunkter i norra o södra delen av modellområdet än vid oroterat gridnät.

I de fall då gridnätet är roterat kan koordinaterna i gridnätet behöva beräknas om till reguljär latitud och longitud. Detta kan göras med hjälp av Fortran-rutinen REGROT, se nedan, eller de ekvationer den innehåller. All indata som behövs för REGROT finns som metadata i GRIB-filen.

Rotera vinden

Vinden beskrivs i de nedladdningsbara GRIB-filerna som U/V (ost/nord) relaterat till det underliggande gridnätet. Är gridnätet roterat är även vinden roterad och för att bli korrekt i ett reguljärt gridnät behöver en omräkning ske.

Denna omräkning kan ske med hjälp av Fortran-rutinen TURNWI, se nedan, eller de ekvationer den innehåller.

Parameter-koder

Metadata för att beskriva vad en parameterkod betyder finns beskrivet i en extern fil. Olika tabeller översätter olika modellers parametrar. Vilken tabellversion som används i en fil anges av metadata-nyckeln table2Version i sektion 1 i GRIB-filen.

Kodtabeller för de modeller som finns nedladdningsbara kan laddas ner som en textfil; Kodlista Arome (3 kB, txt), Kodlista Arome (3 kB, txt), Kodlista Mesan (1 kB, txt), Kodlista PMP (2 kB, txt), Kodlista Hiromb (1 kB, txt).

Gridnät för vindberäkningar

Det finns en skillnad i hur beräkningar av vinden sker på olika nivåer. U och V beräknas på olika typer av ARAKAWA-grid beroende på typ av nivå.

Vinden på 10-meters nivå beräknas enligt ARAKAWA C medan vinden på modellnivåer (indicatorOfTypeOfLevel=109) beräknas enligt ARAKAWA A.

Exempel på rotation och område

En modells hörnkoordinater (beräknat i modellens projektion) och rotation finns beskrivet i form av metadata i GRIB-meddelandet. I GRIB edition 1 anges detta (värden tagna från PMP) bl.a. med metadata ur sektion 2:

latitudeOfFirstGridPoint = -6749 [milligrader]

longitudeOfFirstGridPoint = -7775 [milligrader]

latitudeOfSouthernPole = -30000 [milligrader]

longitudeOfSouthernPole = 15000 [milligrader]

De två första värdena anger i PMP:s fall sydvästra hörnet i modellens projektion (roterad latitud/longitud i PMP). De två sista talar om till vilka koordinater sydpolen har flyttats till i samband med rotationen.

Omräkning mellan roterat och reguljär projektion är komplicerat. Som en förenklad sett att förklara områdets storlek kan man tänka sig att skärningen mellan nollmeridianen och ekvatorn genom rotation flyttats till 15° ost och 60° nord. Utifrån denna punkt kan man grovt ange sydvästra hörnet av området till 6,749° syd och 7,775° väst.

Nordöstra hörnet kan sedan grovt beräknat med hjälp av avståndet mellan gridpunkterna och antalet punkter i x och y-led.

Kod för regrot

      CALL REGROT(ORTLON,ORTLAT,
      +     ALON_R,ALAT_R,
      +      1,1,1,1,ALONSP,ALATSP,1)
c ALON_R och ALAT_R är positionen i det roterade gridet

      YY = 1. + (ALAT_R - ALATS)/DTHETA
      XX = 1. + (ALON_R - ALONW)/DLAMDA

      CALL INTER(NLON,NLAT,XX,YY,GRDVAL,PROG)

      BIAS     = GRDVAL-OBSVAL
      WRITE(*,'(2f7.2,3f10.3)') ORTLAT,ORTLON,BIAS,
 

Fortrandkod för regrot

subroutine regrot(pxreg,pyreg,pxrot,pyrot,kxdim,kydim,kx,ky,
      +           pxcen,pycen,kcall)
c
       implicit none
c
c-------------------------------------------------------------------------------
c
c*    conversion between regular and rotated spherical coordinates.
c*
C*    pxreg       longitudes of the regular coordinates
c*     pyreg       latitudes of the regular coordinates
c*     pxrot        longitudes of the rotated coordinates
c*     pyrot        latitudes of the rotated coordinates
c*                     all coordinates given in degrees N (negative values for S)
c*                     and degrees E (negative values for W)
c*     kxdim       dimension of the gridpoint fields in the x-direction
c*     kydim       dimension of the gridpoint fieLds in the y-direction
c*     kx            number of gridpoint in the x-direction
c*     ky            number of gridpoints in the y-direction
c*     pxcen       regular longitude of the south pole of the rotated grid
c*     pycen        regular latitude of the south pole of the rotated grid
c*
C*    kcall=-1: find regular as functions of rotated coordinates.
c*     kcall= 1: find rotated as functions of regular coordinates.
c*
c*     j.e. haugen hirlam june -92
c
c-------------------------------------------------------------------------------
c
        integer kxdim,kydim,kx,ky,kcall
        real pxreg(kxdim,kydim),pyreg(kxdim,kydim),

        +               pxrot(kxdim,kydim),pyrot(kxdim,kydim),
        +               pxcen,pycen
c
c-------------------------------------------------------------------------------
c
        real pi,zrad,zsycen,zcycen,zxmxc,zsxmxc,zcxmxc,zsyreg,zcyreg,
        x               zsyrot,zcyrot,zcxrot,zsxrot,zradi
        integer jy,jx
c
c-------------------------------------------------------------------------------
c
        pi = 4.*atan(1.)
        zrad = pi/180.
        zradi = 1./zrad
        zsycen = sin(zrad*(pycen+90.))
        zcycen = cos(zrad*(pycen+90.))
c
        if (kcall.eq.1) then
c
        do jy = 1,ky
        do jx = 1,kx
c
        zxmxc = zrad*(pxreg(jx,jy) - pxcen)
        zsxmxc = sin(zxmxc)
        zcxmxc = cos(zxmxc)
        zsyreg = sin(zrad*pyreg(jx,jy))
        zcyreg = cos(zrad*pyreg(jx,jy))
        zsyrot = zcycen*zsyreg - zsycen*zcyreg*zcxmxc
        zsyrot = amax1(zsyrot,-1.0)
        zsyrot = amin1(zsyrot,+1.0)
c
        pyrot(jx,jy) = asin(zsyrot)*zradi
c
        zcyrot = cos(pyrot(jx,jy)*zrad)
        zcxrot = (zcycen*zcyreg*zcxmxc +
     +               zsycen*zsyreg)/zcyrot
        zcxrot = amax1(zcxrot,-1.0)
        zcxrot = amin1(zcxrot,+1.0)
        zsxrot = zcyreg*zsxmxc/zcyrot
c
        pxrot(jx,jy) = acos(zcxrot)*zradi
c
        if (zsxrot.lt.0.0) pxrot(jx,jy) = -pxrot(jx,jy)
c
        enddo
        enddo
c
        elseif (kcall.eq.-1) then
c
        do jy = 1,ky
        do jx = 1,kx
c
        zsxrot = sin(zrad*pxrot(jx,jy))
        zcxrot = cos(zrad*pxrot(jx,jy))
        zsyrot = sin(zrad*pyrot(jx,jy))
        zcyrot = cos(zrad*pyrot(jx,jy))
        zsyreg = zcycen*zsyrot + zsycen*zcyrot*zcxrot
        zsyreg = amax1(zsyreg,-1.0)
        zsyreg = amin1(zsyreg,+1.0)
c
        pyreg(jx,jy) = asin(zsyreg)*zradi
c
        zcyreg = cos(pyreg(jx,jy)*zrad)
        zcxmxc = (zcycen*zcyrot*zcxrot -
      +                zsycen*zsyrot)/zcyreg
        zcxmxc = amax1(zcxmxc,-1.0)
        zcxmxc = amin1(zcxmxc,+1.0)
        zsxmxc = zcyrot*zsxrot/zcyreg
        zxmxc = acos(zcxmxc)*zradi
        if (zsxmxc.lt.0.0) zxmxc = -zxmxc
c
        pxreg(jx,jy) = zxmxc + pxcen
c
        enddo
        enddo
c
        else
        write(6,'(1x,''invalid kcall in regrot'')')
        stop
        endif
c
        return
        end

Fortrandkod för TURNWI

               subroutine turnwi(puarg,pvarg,pures,pvres,pa,pb,pc,pd,
          +           pxreg,pyreg,pxrot,pyrot,
          +           pxcen,pycen,kcall)
c
c*           turn horizontal velocity components between regular and
c*           rotated spherical coordinates.
c*           kcall=-2 or 2: preparation of coefficients.
c*           kcall=-1 or 1: multiplication after preparations.
c*
c*           j.e. haugen hirlam june -92
c
              real puarg,pvarg,
         +            pures,pvres,
         +            pa,pb,
         +            pc,pd,
         +            pxreg,pyreg,
         +            pxrot,pyrot
c
              if (kcall.eq.1 .or. kcall.eq.-1) then
c
c*          multiplication between regular and rotated spherical grid
c
             pures = pa*puarg + pb*pvarg
             pvres = pc*puarg + pd*pvarg
c
             elseif (kcall.eq.2) then
c
c*          precalculate matrixes from regular to rotated grid
c
             zsyc = sin(pycen)
             zcyc = cos(pycen)
c
             zsxreg = sin(pxreg)
             zcxreg = cos(pxreg)
             zsyreg = sin(pyreg)
             zcyreg = cos(pyreg)
c
             zxmxc = pxreg - pxcen
             zsxmxc = sin(zxmxc)
             zcxmxc = cos(zxmxc)
c
             zsxrot = sin(pxrot)
             zcxrot = cos(pxrot)
             zsyrot = sin(pyrot)
             zcyrot = cos(pyrot)
c
             pa = zcyc*zsxmxc*zsxrot + zcxmxc*zcxrot
             pb = zcyc*zcxmxc*zsyreg*zsxrot - zsyc*zcyreg*zsxrot -
        +                   zsxmxc*zsyreg*zcxrot
             pc = zsyc*zsxmxc/zcyrot
             pd = (zsyc*zcxmxc*zsyreg + zcyc*zcyreg)/zcyrot
c
             elseif (kcall.eq.-2) then
c
c*          precalculate matrixes from rotated to regular grid
c
             zsyc = sin(pycen)
             zcyc = cos(pycen)
c
             zsxreg = sin(pxreg)
             zcxreg = cos(pxreg)
             zsyreg = sin(pyreg)
             zcyreg = cos(pyreg)
c
             zxmxc = pxreg - pxcen
             zsxmxc = sin(zxmxc)
             zcxmxc = cos(zxmxc)
c
             zsxrot = sin(pxrot)
             zcxrot = cos(pxrot)
             zsyrot = sin(pyrot)
             zcyrot = cos(pyrot)
c
             pa = zcxmxc*zcxrot + zcyc*zsxmxc*zsxrot
             pb = zcyc*zsxmxc*zcxrot*zsyrot + zsyc*zsxmxc*zcyrot -
         +                  zcxmxc*zsxrot*zsyrot
             pc =-zsyc*zsxrot/zcyreg
             pd = (zcyc*zcyrot - zsyc*zcxrot*zsyrot)/zcyreg
c
             else
             write(6,'(1x,''invalid kcall in turnwi'')')
             stop
             endif
c
            return
            end