128 lines
4 KiB
C
128 lines
4 KiB
C
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <math.h>
|
|
#include "grib2.h"
|
|
|
|
|
|
void specpack(g2float *fld,g2int ndpts,g2int JJ,g2int KK,g2int MM,
|
|
g2int *idrstmpl,unsigned char *cpack,g2int *lcpack)
|
|
//$$$ SUBPROGRAM DOCUMENTATION BLOCK
|
|
// . . . .
|
|
// SUBPROGRAM: specpack
|
|
// PRGMMR: Gilbert ORG: W/NP11 DATE: 2002-12-19
|
|
//
|
|
// ABSTRACT: This subroutine packs a spectral data field using the complex
|
|
// packing algorithm for spherical harmonic data as
|
|
// defined in the GRIB2 Data Representation Template 5.51.
|
|
//
|
|
// PROGRAM HISTORY LOG:
|
|
// 2002-12-19 Gilbert
|
|
//
|
|
// USAGE: void specpack(g2float *fld,g2int ndpts,g2int JJ,g2int KK,g2int MM,
|
|
// g2int *idrstmpl,insigned char *cpack,g2int *lcpack)
|
|
// INPUT ARGUMENT LIST:
|
|
// fld[] - Contains the packed data values
|
|
// ndpts - The number of data values to pack
|
|
// JJ - J - pentagonal resolution parameter
|
|
// KK - K - pentagonal resolution parameter
|
|
// MM - M - pentagonal resolution parameter
|
|
// idrstmpl - Contains the array of values for Data Representation
|
|
// Template 5.51
|
|
//
|
|
// OUTPUT ARGUMENT LIST:
|
|
// cpack - The packed data field (character*1 array)
|
|
// lcpack - length of packed field cpack().
|
|
//
|
|
// REMARKS: None
|
|
//
|
|
// ATTRIBUTES:
|
|
// LANGUAGE: C
|
|
// MACHINE: IBM SP
|
|
//
|
|
//$$$
|
|
{
|
|
|
|
g2int *ifld,tmplsim[5];
|
|
g2float bscale,dscale,*unpk,*tfld;
|
|
g2float *pscale,tscale;
|
|
g2int Js,Ks,Ms,Ts,Ns,inc,incu,incp,n,Nm,m,ipos;
|
|
|
|
bscale = int_power(2.0,-idrstmpl[1]);
|
|
dscale = int_power(10.0,idrstmpl[2]);
|
|
Js=idrstmpl[5];
|
|
Ks=idrstmpl[6];
|
|
Ms=idrstmpl[7];
|
|
Ts=idrstmpl[8];
|
|
|
|
//
|
|
// Calculate Laplacian scaling factors for each possible wave number.
|
|
//
|
|
pscale=(g2float *)malloc((JJ+MM)*sizeof(g2float));
|
|
tscale=(g2float)idrstmpl[4]*1E-6;
|
|
for (n=Js;n<=JJ+MM;n++)
|
|
pscale[n]=pow((g2float)(n*(n+1)),tscale);
|
|
//
|
|
// Separate spectral coeffs into two lists; one to contain unpacked
|
|
// values within the sub-spectrum Js, Ks, Ms, and the other with values
|
|
// outside of the sub-spectrum to be packed.
|
|
//
|
|
tfld=(g2float *)malloc(ndpts*sizeof(g2float));
|
|
unpk=(g2float *)malloc(ndpts*sizeof(g2float));
|
|
ifld=(g2int *)malloc(ndpts*sizeof(g2int));
|
|
inc=0;
|
|
incu=0;
|
|
incp=0;
|
|
for (m=0;m<=MM;m++) {
|
|
Nm=JJ; // triangular or trapezoidal
|
|
if ( KK == JJ+MM ) Nm=JJ+m; // rhombodial
|
|
Ns=Js; // triangular or trapezoidal
|
|
if ( Ks == Js+Ms ) Ns=Js+m; // rhombodial
|
|
for (n=m;n<=Nm;n++) {
|
|
if (n<=Ns && m<=Ms) { // save unpacked value
|
|
unpk[incu++]=fld[inc++]; // real part
|
|
unpk[incu++]=fld[inc++]; // imaginary part
|
|
}
|
|
else { // Save value to be packed and scale
|
|
// Laplacian scale factor
|
|
tfld[incp++]=fld[inc++]*pscale[n]; // real part
|
|
tfld[incp++]=fld[inc++]*pscale[n]; // imaginary part
|
|
}
|
|
}
|
|
}
|
|
|
|
free(pscale);
|
|
|
|
if (incu != Ts) {
|
|
printf("specpack: Incorrect number of unpacked values %d given:\n",(int)Ts);
|
|
printf("specpack: Resetting idrstmpl[8] to %d\n",(int)incu);
|
|
Ts=incu;
|
|
}
|
|
//
|
|
// Add unpacked values to the packed data array in 32-bit IEEE format
|
|
//
|
|
mkieee(unpk,(g2int *)cpack,Ts);
|
|
ipos=4*Ts;
|
|
//
|
|
// Scale and pack the rest of the coefficients
|
|
//
|
|
tmplsim[1]=idrstmpl[1];
|
|
tmplsim[2]=idrstmpl[2];
|
|
tmplsim[3]=idrstmpl[3];
|
|
simpack(tfld,ndpts-Ts,tmplsim,cpack+ipos,lcpack);
|
|
*lcpack=(*lcpack)+ipos;
|
|
//
|
|
// Fill in Template 5.51
|
|
//
|
|
idrstmpl[0]=tmplsim[0];
|
|
idrstmpl[1]=tmplsim[1];
|
|
idrstmpl[2]=tmplsim[2];
|
|
idrstmpl[3]=tmplsim[3];
|
|
idrstmpl[8]=Ts;
|
|
idrstmpl[9]=1; // Unpacked spectral data is 32-bit IEEE
|
|
|
|
free(tfld);
|
|
free(unpk);
|
|
free(ifld);
|
|
|
|
return;
|
|
}
|