1 : #include <stdio.h>
2 : #include <stdlib.h>
3 : #include <math.h>
4 : #include "grib2.h"
5 :
6 :
7 0 : void specpack(g2float *fld,g2int ndpts,g2int JJ,g2int KK,g2int MM,
8 : g2int *idrstmpl,unsigned char *cpack,g2int *lcpack)
9 : //$$$ SUBPROGRAM DOCUMENTATION BLOCK
10 : // . . . .
11 : // SUBPROGRAM: specpack
12 : // PRGMMR: Gilbert ORG: W/NP11 DATE: 2002-12-19
13 : //
14 : // ABSTRACT: This subroutine packs a spectral data field using the complex
15 : // packing algorithm for spherical harmonic data as
16 : // defined in the GRIB2 Data Representation Template 5.51.
17 : //
18 : // PROGRAM HISTORY LOG:
19 : // 2002-12-19 Gilbert
20 : //
21 : // USAGE: void specpack(g2float *fld,g2int ndpts,g2int JJ,g2int KK,g2int MM,
22 : // g2int *idrstmpl,insigned char *cpack,g2int *lcpack)
23 : // INPUT ARGUMENT LIST:
24 : // fld[] - Contains the packed data values
25 : // ndpts - The number of data values to pack
26 : // JJ - J - pentagonal resolution parameter
27 : // KK - K - pentagonal resolution parameter
28 : // MM - M - pentagonal resolution parameter
29 : // idrstmpl - Contains the array of values for Data Representation
30 : // Template 5.51
31 : //
32 : // OUTPUT ARGUMENT LIST:
33 : // cpack - The packed data field (character*1 array)
34 : // lcpack - length of packed field cpack().
35 : //
36 : // REMARKS: None
37 : //
38 : // ATTRIBUTES:
39 : // LANGUAGE: C
40 : // MACHINE: IBM SP
41 : //
42 : //$$$
43 : {
44 :
45 : g2int *ifld,tmplsim[5];
46 : g2float bscale,dscale,*unpk,*tfld;
47 : g2float *pscale,tscale;
48 : g2int Js,Ks,Ms,Ts,Ns,inc,incu,incp,n,Nm,m,ipos;
49 :
50 0 : bscale = int_power(2.0,-idrstmpl[1]);
51 0 : dscale = int_power(10.0,idrstmpl[2]);
52 0 : Js=idrstmpl[5];
53 0 : Ks=idrstmpl[6];
54 0 : Ms=idrstmpl[7];
55 0 : Ts=idrstmpl[8];
56 :
57 : //
58 : // Calculate Laplacian scaling factors for each possible wave number.
59 : //
60 0 : pscale=(g2float *)malloc((JJ+MM)*sizeof(g2float));
61 0 : tscale=(g2float)idrstmpl[4]*1E-6;
62 0 : for (n=Js;n<=JJ+MM;n++)
63 0 : pscale[n]=pow((g2float)(n*(n+1)),tscale);
64 : //
65 : // Separate spectral coeffs into two lists; one to contain unpacked
66 : // values within the sub-spectrum Js, Ks, Ms, and the other with values
67 : // outside of the sub-spectrum to be packed.
68 : //
69 0 : tfld=(g2float *)malloc(ndpts*sizeof(g2float));
70 0 : unpk=(g2float *)malloc(ndpts*sizeof(g2float));
71 0 : ifld=(g2int *)malloc(ndpts*sizeof(g2int));
72 0 : inc=0;
73 0 : incu=0;
74 0 : incp=0;
75 0 : for (m=0;m<=MM;m++) {
76 0 : Nm=JJ; // triangular or trapezoidal
77 0 : if ( KK == JJ+MM ) Nm=JJ+m; // rhombodial
78 0 : Ns=Js; // triangular or trapezoidal
79 0 : if ( Ks == Js+Ms ) Ns=Js+m; // rhombodial
80 0 : for (n=m;n<=Nm;n++) {
81 0 : if (n<=Ns && m<=Ms) { // save unpacked value
82 0 : unpk[incu++]=fld[inc++]; // real part
83 0 : unpk[incu++]=fld[inc++]; // imaginary part
84 : }
85 : else { // Save value to be packed and scale
86 : // Laplacian scale factor
87 0 : tfld[incp++]=fld[inc++]*pscale[n]; // real part
88 0 : tfld[incp++]=fld[inc++]*pscale[n]; // imaginary part
89 : }
90 : }
91 : }
92 :
93 0 : free(pscale);
94 :
95 0 : if (incu != Ts) {
96 0 : printf("specpack: Incorrect number of unpacked values %d given:\n",(int)Ts);
97 0 : printf("specpack: Resetting idrstmpl[8] to %d\n",(int)incu);
98 0 : Ts=incu;
99 : }
100 : //
101 : // Add unpacked values to the packed data array in 32-bit IEEE format
102 : //
103 0 : mkieee(unpk,(g2int *)cpack,Ts);
104 0 : ipos=4*Ts;
105 : //
106 : // Scale and pack the rest of the coefficients
107 : //
108 0 : tmplsim[1]=idrstmpl[1];
109 0 : tmplsim[2]=idrstmpl[2];
110 0 : tmplsim[3]=idrstmpl[3];
111 0 : simpack(tfld,ndpts-Ts,tmplsim,cpack+ipos,lcpack);
112 0 : *lcpack=(*lcpack)+ipos;
113 : //
114 : // Fill in Template 5.51
115 : //
116 0 : idrstmpl[0]=tmplsim[0];
117 0 : idrstmpl[1]=tmplsim[1];
118 0 : idrstmpl[2]=tmplsim[2];
119 0 : idrstmpl[3]=tmplsim[3];
120 0 : idrstmpl[8]=Ts;
121 0 : idrstmpl[9]=1; // Unpacked spectral data is 32-bit IEEE
122 :
123 0 : free(tfld);
124 0 : free(unpk);
125 0 : free(ifld);
126 :
127 : return;
128 : }
|
