[2759] | 1 | subroutine specunpack(cpack,len,idrstmpl,ndpts,JJ,KK,MM,fld) |
---|
| 2 | !$$$ SUBPROGRAM DOCUMENTATION BLOCK |
---|
| 3 | ! . . . . |
---|
| 4 | ! SUBPROGRAM: specunpack |
---|
| 5 | ! PRGMMR: Gilbert ORG: W/NP11 DATE: 2002-12-19 |
---|
| 6 | ! |
---|
| 7 | ! ABSTRACT: This subroutine unpacks a spectral data field that was packed |
---|
| 8 | ! using the complex packing algorithm for spherical harmonic data as |
---|
| 9 | ! defined in the GRIB2 documention, |
---|
| 10 | ! using info from the GRIB2 Data Representation Template 5.51. |
---|
| 11 | ! |
---|
| 12 | ! PROGRAM HISTORY LOG: |
---|
| 13 | ! 2002-12-19 Gilbert |
---|
| 14 | ! |
---|
| 15 | ! USAGE: CALL specunpack(cpack,len,idrstmpl,ndpts,JJ,KK,MM,fld) |
---|
| 16 | ! INPUT ARGUMENT LIST: |
---|
| 17 | ! cpack - The packed data field (character*1 array) |
---|
| 18 | ! len - length of packed field cpack(). |
---|
| 19 | ! idrstmpl - Contains the array of values for Data Representation |
---|
| 20 | ! Template 5.51 |
---|
| 21 | ! ndpts - The number of data values to unpack |
---|
| 22 | ! JJ - J - pentagonal resolution parameter |
---|
| 23 | ! KK - K - pentagonal resolution parameter |
---|
| 24 | ! MM - M - pentagonal resolution parameter |
---|
| 25 | ! |
---|
| 26 | ! OUTPUT ARGUMENT LIST: |
---|
| 27 | ! fld() - Contains the unpacked data values |
---|
| 28 | ! |
---|
| 29 | ! REMARKS: None |
---|
| 30 | ! |
---|
| 31 | ! ATTRIBUTES: |
---|
| 32 | ! LANGUAGE: XL Fortran 90 |
---|
| 33 | ! MACHINE: IBM SP |
---|
| 34 | ! |
---|
| 35 | !$$$ |
---|
| 36 | |
---|
| 37 | character(len=1),intent(in) :: cpack(len) |
---|
| 38 | integer,intent(in) :: ndpts,len,JJ,KK,MM |
---|
| 39 | integer,intent(in) :: idrstmpl(*) |
---|
| 40 | real,intent(out) :: fld(ndpts) |
---|
| 41 | |
---|
| 42 | integer :: ifld(ndpts),Ts |
---|
| 43 | integer(4) :: ieee |
---|
| 44 | real :: ref,bscale,dscale,unpk(ndpts) |
---|
| 45 | real,allocatable :: pscale(:) |
---|
| 46 | |
---|
| 47 | ieee = idrstmpl(1) |
---|
| 48 | call rdieee(ieee,ref,1) |
---|
| 49 | bscale = 2.0**real(idrstmpl(2)) |
---|
| 50 | dscale = 10.0**real(-idrstmpl(3)) |
---|
| 51 | nbits = idrstmpl(4) |
---|
| 52 | Js=idrstmpl(6) |
---|
| 53 | Ks=idrstmpl(7) |
---|
| 54 | Ms=idrstmpl(8) |
---|
| 55 | Ts=idrstmpl(9) |
---|
| 56 | |
---|
| 57 | if (idrstmpl(10).eq.1) then ! unpacked floats are 32-bit IEEE |
---|
| 58 | !call g2lib_gbytes(cpack,ifld,0,32,0,Ts) |
---|
| 59 | call rdieee(cpack,unpk,Ts) ! read IEEE unpacked floats |
---|
| 60 | iofst=32*Ts |
---|
| 61 | call g2lib_gbytes(cpack,ifld,iofst,nbits,0,ndpts-Ts) ! unpack scaled data |
---|
| 62 | ! |
---|
| 63 | ! Calculate Laplacian scaling factors for each possible wave number. |
---|
| 64 | ! |
---|
| 65 | allocate(pscale(JJ+MM)) |
---|
| 66 | tscale=real(idrstmpl(5))*1E-6 |
---|
| 67 | do n=Js,JJ+MM |
---|
| 68 | pscale(n)=real(n*(n+1))**(-tscale) |
---|
| 69 | enddo |
---|
| 70 | ! |
---|
| 71 | ! Assemble spectral coeffs back to original order. |
---|
| 72 | ! |
---|
| 73 | inc=1 |
---|
| 74 | incu=1 |
---|
| 75 | incp=1 |
---|
| 76 | do m=0,MM |
---|
| 77 | Nm=JJ ! triangular or trapezoidal |
---|
| 78 | if ( KK .eq. JJ+MM ) Nm=JJ+m ! rhombodial |
---|
| 79 | Ns=Js ! triangular or trapezoidal |
---|
| 80 | if ( Ks .eq. Js+Ms ) Ns=Js+m ! rhombodial |
---|
| 81 | do n=m,Nm |
---|
| 82 | if (n.le.Ns .AND. m.le.Ms) then ! grab unpacked value |
---|
| 83 | fld(inc)=unpk(incu) ! real part |
---|
| 84 | fld(inc+1)=unpk(incu+1) ! imaginary part |
---|
| 85 | inc=inc+2 |
---|
| 86 | incu=incu+2 |
---|
| 87 | else ! Calc coeff from packed value |
---|
| 88 | fld(inc)=((real(ifld(incp))*bscale)+ref)* |
---|
| 89 | & dscale*pscale(n) ! real part |
---|
| 90 | fld(inc+1)=((real(ifld(incp+1))*bscale)+ref)* |
---|
| 91 | & dscale*pscale(n) ! imaginary part |
---|
| 92 | inc=inc+2 |
---|
| 93 | incp=incp+2 |
---|
| 94 | endif |
---|
| 95 | enddo |
---|
| 96 | enddo |
---|
| 97 | |
---|
| 98 | deallocate(pscale) |
---|
| 99 | |
---|
| 100 | else |
---|
| 101 | print *,'specunpack: Cannot handle 64 or 128-bit floats.' |
---|
| 102 | fld=0.0 |
---|
| 103 | return |
---|
| 104 | endif |
---|
| 105 | |
---|
| 106 | return |
---|
| 107 | end |
---|