g2spec.F90

Go to the documentation of this file.

 
 
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
 
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