-
Notifications
You must be signed in to change notification settings - Fork 17
/
Copy pathg2spec.F90
185 lines (166 loc) · 6.22 KB
/
g2spec.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
!> @file
!> @brief Pack/unpack a spectral data field using the complex
!> packing algorithm for spherical harmonic data as defined in
!> [Data Representation Template
!> 5.51](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-51.shtml).
!> @author Stephen Gilbert @date 2002-12-19
!> Pack a spectral data field using the complex
!> packing algorithm for spherical harmonic data as defined in
!> [Data Representation Template
!> 5.51](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-51.shtml).
!>
!> @param[in] fld Contains the data values to pack.
!> @param[in] ndpts The number of data values in array fld.
!> @param[in] JJ J pentagonal resolution parameter.
!> @param[in] KK K pentagonal resolution parameter.
!> @param[in] MM M pentagonal resolution parameter.
!> @param[in] idrstmpl Contains the array of values for Data
!> Representation Template 5.51.
!> @param[out] cpack The packed data field (character*1 array).
!> @param[out] lcpack length of packed field cpack.
!>
!> @author Stephen Gilbert @date 2002-12-19
subroutine specpack(fld,ndpts,JJ,KK,MM,idrstmpl,cpack,lcpack)
real,intent(in) :: fld(ndpts)
integer,intent(in) :: ndpts,JJ,KK,MM
integer,intent(inout) :: idrstmpl(*)
character(len=1),intent(out) :: cpack(*)
integer,intent(out) :: lcpack
integer :: Ts,tmplsim(5)
real :: bscale,dscale,unpk(ndpts),tfld(ndpts)
real,allocatable :: pscale(:)
bscale = 2.0**real(-idrstmpl(2))
dscale = 10.0**real(idrstmpl(3))
nbits = idrstmpl(4)
Js=idrstmpl(6)
Ks=idrstmpl(7)
Ms=idrstmpl(8)
Ts=idrstmpl(9)
! Calculate Laplacian scaling factors for each possible wave number.
allocate(pscale(JJ+MM))
tscale=real(idrstmpl(5))*1E-6
do n=Js,JJ+MM
pscale(n)=real(n*(n+1))**(tscale)
enddo
! 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.
inc=1
incu=1
incp=1
do m=0,MM
Nm=JJ ! triangular or trapezoidal
if (KK .eq. JJ+MM) Nm=JJ+m ! rhombodial
Ns=Js ! triangular or trapezoidal
if (Ks .eq. Js+Ms) Ns=Js+m ! rhombodial
do n=m,Nm
if (n.le.Ns .AND. m.le.Ms) then ! save unpacked value
unpk(incu)=fld(inc) ! real part
unpk(incu+1)=fld(inc+1) ! imaginary part
inc=inc+2
incu=incu+2
else
! Save value to be packed and scale Laplacian scale factor.
tfld(incp)=fld(inc)*pscale(n) ! real part
tfld(incp+1)=fld(inc+1)*pscale(n) ! imaginary part
inc=inc+2
incp=incp+2
endif
enddo
enddo
deallocate(pscale)
incu=incu-1
if (incu .ne. Ts) then
print *,'specpack: Incorrect number of unpacked values ', 'given:',Ts
print *,'specpack: Resetting idrstmpl(9) to ',incu
Ts=incu
endif
! Add unpacked values to the packed data array in 32-bit IEEE format.
call mkieee(unpk,cpack,Ts)
ipos=4*Ts
! Scale and pack the rest of the coefficients.
tmplsim(2)=idrstmpl(2)
tmplsim(3)=idrstmpl(3)
tmplsim(4)=idrstmpl(4)
call simpack(tfld,ndpts-Ts,tmplsim,cpack(ipos+1),lcpack)
lcpack=lcpack+ipos
! Fill in Template 5.51.
idrstmpl(1)=tmplsim(1)
idrstmpl(2)=tmplsim(2)
idrstmpl(3)=tmplsim(3)
idrstmpl(4)=tmplsim(4)
idrstmpl(9)=Ts
idrstmpl(10)=1 ! Unpacked spectral data is 32-bit IEEE
end subroutine specpack
!> Unpack a spectral data field using the complex packing algorithm
!> for spherical harmonic data, [Data Representation Template
!> 5.51](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-51.shtml).
!>
!> @param[in] cpack The packed data field (character*1 array).
!> @param[in] len length of packed field cpack.
!> @param[in] idrstmpl Contains the array of values for Data
!> Representation Template 5.51.
!> @param[in] ndpts The number of data values in array fld.
!> @param[in] JJ J pentagonal resolution parameter.
!> @param[in] KK K pentagonal resolution parameter.
!> @param[in] MM M pentagonal resolution parameter.
!> @param[out] fld Contains the unpacked data values.
!>
!> @author Stephen Gilbert @date 2002-12-19
subroutine specunpack(cpack, len, idrstmpl, ndpts, JJ, KK, MM, fld)
character(len = 1), intent(in) :: cpack(len)
integer, intent(in) :: ndpts, len, JJ, KK, MM
integer, intent(in) :: idrstmpl(*)
real, intent(out) :: fld(ndpts)
integer :: ifld(ndpts), Ts
integer(4) :: ieee
real :: ref, bscale, dscale, unpk(ndpts)
real, allocatable :: pscale(:)
ieee = idrstmpl(1)
call rdieee(ieee, ref, 1)
bscale = 2.0**real(idrstmpl(2))
dscale = 10.0**real(-idrstmpl(3))
nbits = idrstmpl(4)
Js = idrstmpl(6)
Ks = idrstmpl(7)
Ms = idrstmpl(8)
Ts = idrstmpl(9)
if (idrstmpl(10) .eq. 1) then ! unpacked floats are 32-bit IEEE
call rdieeec(cpack, unpk, Ts) ! read IEEE unpacked floats
iofst = 32 * Ts
call g2_gbytesc(cpack, ifld, iofst, nbits, 0, ndpts - Ts) ! unpack scaled data
! Calculate Laplacian scaling factors for each possible wave number.
allocate(pscale(JJ + MM))
tscale = real(idrstmpl(5)) * 1E-6
do n = Js, JJ + MM
pscale(n) = real(n * (n + 1))**(-tscale)
enddo
! Assemble spectral coeffs back to original order.
inc = 1
incu = 1
incp = 1
do m = 0, MM
Nm = JJ ! triangular or trapezoidal
if (KK .eq. JJ + MM) Nm = JJ + m ! rhombodial
Ns = Js ! triangular or trapezoidal
if (Ks .eq. Js + Ms) Ns = Js + m ! rhombodial
do n = m, Nm
if (n .le. Ns .AND. m .le. Ms) then ! grab unpacked value
fld(inc) = unpk(incu) ! real part
fld(inc + 1) = unpk(incu + 1) ! imaginary part
inc = inc + 2
incu = incu + 2
else ! Calc coeff from packed value
fld(inc) = ((real(ifld(incp))*bscale) + ref)*dscale*pscale(n) ! real part
fld(inc + 1) = ((real(ifld(incp + 1)) * bscale) + ref) * dscale * pscale(n) ! imaginary part
inc = inc + 2
incp = incp + 2
endif
enddo
enddo
deallocate(pscale)
else
print *, 'specunpack: Cannot handle 64 or 128-bit floats.'
fld = 0.0
endif
end subroutine specunpack