subroutine general_sptez_s(sp,wave,grid,idir)
!$$$ subprogram documentation block
! . . . .
! subprogram: sptez_s perform a simple scalar spherical transform
! prgmmr: iredell org: np23 date: 1996-02-29
!
! absract: this subprogram performs a spherical transform
! between spectral coefficients of a scalar quantity
! and a field on a global cylindrical grid.
! the wave-space can be either triangular or rhomboidal.
! the grid-space can be either an equally-spaced grid
! (with or without pole points) or a gaussian grid.
! the wave field is in sequential 'ibm order'.
! the grid field is indexed east to west, then north to south.
! for more flexibility and efficiency, call sptran.
! subprogram can be called from a multiprocessing environment.
!
! This routine differs from splib routine sptez in that
! 1) the calling list only contains the in/out arrays and
! flag for the direction in which to transform
! 2) it calls a version of sptranf that does not invoke
! initialization routines on each entry
! 3) some generality built into the splib version is
! removed in the code below
!
! program history log:
! 1996-02-29 iredell
! 2004-08-23 treadon - adapt splib routine sptez for gsi use
! 2007-04-25 errico - replace use of duplicate arguments in sptranf_s
! 2010-02-18 parrish - copy to general_sptez_s, and pass specmod vars through
! input variable sp of type(spec_vars)
! input arguments:
! wave - real (2*mx) wave field if idir>0
! where mx=(jcap+1)*((iromb+1)*jcap+2)/2
! grid - real (imax,jmax) grid field (e->w,n->s) if idir<0
! idir - integer transform flag
! (idir>0 for wave to grid, idir<0 for grid to wave)
!
! output arguments:
! wave - real (2*mx) wave field if idir<0
! where mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
! grid - real (imax,jmax) grid field (e->w,n->s) if idir>0
!
! subprograms called:
! sptranf_s - perform a scalar spherical transform
!
! remarks: minimum grid dimensions for unaliased transforms to spectral:
! dimension linear quadratic
! ----------------------- --------- -------------
! imax 2*maxwv+2 3*maxwv/2*2+2
! jmax (idrt=4,iromb=0) 1*maxwv+1 3*maxwv/2+1
! jmax (idrt=4,iromb=1) 2*maxwv+1 5*maxwv/2+1
! jmax (idrt=0,iromb=0) 2*maxwv+3 3*maxwv/2*2+3
! jmax (idrt=0,iromb=1) 4*maxwv+3 5*maxwv/2*2+3
! jmax (idrt=256,iromb=0) 2*maxwv+1 3*maxwv/2*2+1
! jmax (idrt=256,iromb=1) 4*maxwv+1 5*maxwv/2*2+1
! ----------------------- --------- -------------
!
! attributes:
! language: fortran 77
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: zero
use general_specmod, only: spec_vars
implicit none
! Declare passed variables
type(spec_vars) ,intent(in ) :: sp
integer(i_kind) ,intent(in ) :: idir
real(r_kind),dimension(sp%nc) ,intent(inout) :: wave
real(r_kind),dimension(sp%ijmax),intent(inout) :: grid
! Declare local variables
integer(i_kind) i
! Zero appropriate output array based on direction of transform
if (idir<0) then
do i=1,sp%nc
wave(i)=zero
end do
elseif (idir>0) then
do i=1,sp%ijmax
grid(i)=zero
end do
endif
! Call spectral <--> grid transform
call general_sptranf_s(sp,wave,grid,idir)
return
end subroutine general_sptez_s
subroutine general_sptranf_s(sp_a,wave,grid,idir)
!$$$ subprogram documentation block
! . . . .
! subprogram: sptranf_s perform a scalar spherical transform
! prgmmr: iredell org: np23 date: 1996-02-29
!
! abstract: this subprogram performs a spherical transform
! between spectral coefficients of scalar quantities
! and fields on a global cylindrical grid.
! the wave-space can be either triangular or rhomboidal.
! the grid-space can be either an equally-spaced grid
! (with or without pole points) or a gaussian grid.
! the wave and grid fields may have general indexing,
! but each wave field is in sequential 'ibm order',
! i.e. with zonal wavenumber as the slower index.
! transforms are done in latitude pairs for efficiency;
! thus grid arrays for each hemisphere must be passed.
! if so requested, just a subset of the latitude pairs
! may be transformed in each invocation of the subprogram.
! the transforms are all multiprocessed over latitude except
! the transform from fourier to spectral is multiprocessed
! over zonal wavenumber to ensure reproducibility.
! transform several fields at a time to improve vectorization.
! subprogram can be called from a multiprocessing environment.
!
! This routine differs from splib routine sptranf in that
! it does not call sptranf0 (an initialization routine).
!
! program history log:
! 1996-02-29 iredell
! 1998-12-15 iredell generic fft used
! 2004-08-23 treadon - adapt splib routine sptranf for gsi use
! 2006-05-03 treadon - remove jc from specmod list since not used
! 2007-04-25 errico - replace use of duplicate arguments
! 2008-04-03 safford - rm unused vars and uses
! 2010-02-18 parrish - copy to general_sptranf_s, and pass specmod vars through
! input variable sp of type(spec_vars)
! 2010-03-31 treadon - add double FFT capabilty (from Joe Sela, Mark Iredell)
!
! input arguments:
! wave - real (*) wave fields if idir>0
! grid - real (*) grid fields (starting at jb) if idir<0
! idir - integer transform flag
! (idir>0 for wave to grid, idir<0 for grid to wave)
!
! output arguments:
! wave - real (*) wave fields if idir<0
! grid - real (*) grid fields (starting at jb) if idir>0
!
! subprograms called:
! sptranf1 sptranf spectral transform
!
! remarks:
! This routine assumes that splib routine sptranf0 has been
! previously called. sptranf0 initializes arrays needed in
! the transforms.
!
! minimum grid dimensions for unaliased transforms to spectral:
! dimension linear quadratic
! ----------------------- --------- -------------
! imax 2*maxwv+2 3*maxwv/2*2+2
! jmax (idrt=4,iromb=0) 1*maxwv+1 3*maxwv/2+1
! jmax (idrt=4,iromb=1) 2*maxwv+1 5*maxwv/2+1
! jmax (idrt=0,iromb=0) 2*maxwv+3 3*maxwv/2*2+3
! jmax (idrt=0,iromb=1) 4*maxwv+3 5*maxwv/2*2+3
! jmax (idrt=256,iromb=0) 2*maxwv+1 3*maxwv/2*2+1
! jmax (idrt=256,iromb=1) 4*maxwv+1 5*maxwv/2*2+1
! ----------------------- --------- -------------
!
! attributes:
! language: fortran 77
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: zero
use general_specmod, only: spec_vars
implicit none
! Declare passed variables
type(spec_vars) ,intent(in ) :: sp_a
integer(i_kind) ,intent(in ) :: idir
real(r_kind),dimension(sp_a%nc) ,intent(inout) :: wave
real(r_kind),dimension(sp_a%ijmax),intent(inout) :: grid
! Declare local variables
integer(i_kind) i,j,jj,ijn,ijs,mp,kw,kwtop,imaxp2
real(r_kind),dimension(2*(sp_a%jcap+1)):: wtop
real(r_kind),dimension(sp_a%imax,2):: g
real(r_kind),dimension(sp_a%imax+2,2):: f
real(r_kind),dimension(50000+4*sp_a%imax):: tmpafft
! Initialize local variables
mp=0
do i=1,2*(sp_a%jcap+1)
wtop(i)=zero
end do
! Transform wave to grid
! ***NOTE***
! The FFT used in the transform below has been generalized to
! allow for projection of spectral coefficients onto double
! the desired number of longitudinal grid points. This
! approach is needed when transforming high wavenumber spectral
! coefficients to a coarser resoultion grid. For example, using
! splib to directly transform T878 spectral coefficients to an
! 1152 x 576 grid does not use Fourier modes above wavenumber 576.
! Joe Sela insightfully suggested doubling the number of points
! in the FFT and using every other point in the output grid.
! Mark Iredell coded up Joe's idea below.
tmpafft(:)=sp_a%afft(:)
kw=(sp_a%jcap+1)*((sp_a%iromb+1)*sp_a%jcap+2)
kwtop=2*(sp_a%jcap+1)
imaxp2=sp_a%imax+2
if(idir>0) then
do j=sp_a%jb,sp_a%je
call spsynth(sp_a%iromb,sp_a%jcap,sp_a%imax,imaxp2,kw,kwtop,1, &
sp_a%clat(j),sp_a%pln(1,j),sp_a%plntop(1,j),0,wave,wtop,f)
call spffte(sp_a%imax,imaxp2/2,sp_a%imax,2,f,g,1,tmpafft)
! call sptranf1(sp_a%iromb,sp_a%jcap,sp_a%idrt,sp_a%imax,sp_a%jmax,j,j, &
! sp_a%eps,sp_a%epstop,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! sp_a%pln(1,j),sp_a%plntop(1,j),mp, &
! wave,wtop,g,idir)
jj = j-sp_a%jb
ijn = jj*sp_a%jn
ijs = jj*sp_a%js + sp_a%ioffset
do i=1,sp_a%imax
grid(ijn+i)=g(i,1)
grid(ijs+i)=g(i,2)
enddo
enddo
! Transform grid to wave
! ***WARNING***
! The code above has been generalized to handle the transform of
! high spectral representation fields to coarse physical space
! grids. The code below should not be used to transform coarse
! resolution grids to high spectral representation. Since this
! functionality is not yet needed in the GSI, the prudent action
! to take here is to print an ERROR message and terminate program
! execution if such a transform is requested.
else
do j=sp_a%jb,sp_a%je
if(sp_a%wlat(j)>zero) then
jj = j-sp_a%jb
ijn = jj*sp_a%jn
ijs = jj*sp_a%js + sp_a%ioffset
do i=1,sp_a%imax
g(i,1)=grid(ijn+i)
g(i,2)=grid(ijs+i)
enddo
call spffte(sp_a%imax,imaxp2/2,sp_a%imax,2,f,g,-1,tmpafft)
call spanaly(sp_a%iromb,sp_a%jcap,sp_a%imax,imaxp2,kw,kwtop,1, &
sp_a%wlat(j),sp_a%clat(j),sp_a%pln(1,j),sp_a%plntop(1,j),0,f,wave,wtop)
! call sptranf1(sp_a%iromb,sp_a%jcap,sp_a%idrt,sp_a%imax,sp_a%jmax,j,j, &
! sp_a%eps,sp_a%epstop,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! sp_a%pln(1,j),sp_a%plntop(1,j),mp, &
! wave,wtop,g,idir)
endif
enddo
endif
end subroutine general_sptranf_s
subroutine general_sptranf_s_b(sp_a,sp_b,wave,grid,idir)
!$$$ subprogram documentation block
! . . . .
! subprogram: sptranf_s perform a scalar spherical transform
! prgmmr: iredell org: np23 date: 1996-02-29
!
! abstract: this subprogram performs a spherical transform
! between spectral coefficients of scalar quantities
! and fields on a global cylindrical grid.
! the wave-space can be either triangular or rhomboidal.
! the grid-space can be either an equally-spaced grid
! (with or without pole points) or a gaussian grid.
! the wave and grid fields may have general indexing,
! but each wave field is in sequential 'ibm order',
! i.e. with zonal wavenumber as the slower index.
! transforms are done in latitude pairs for efficiency;
! thus grid arrays for each hemisphere must be passed.
! if so requested, just a subset of the latitude pairs
! may be transformed in each invocation of the subprogram.
! the transforms are all multiprocessed over latitude except
! the transform from fourier to spectral is multiprocessed
! over zonal wavenumber to ensure reproducibility.
! transform several fields at a time to improve vectorization.
! subprogram can be called from a multiprocessing environment.
!
! This routine differs from splib routine sptranf in that
! it does not call sptranf0 (an initialization routine).
!
! program history log:
! 1996-02-29 iredell
! 1998-12-15 iredell generic fft used
! 2004-08-23 treadon - adapt splib routine sptranf for gsi use
! 2006-05-03 treadon - remove jc from specmod list since not used
! 2007-04-25 errico - replace use of duplicate arguments
! 2008-04-03 safford - rm unused vars and uses
! 2010-02-18 parrish - copy to general_sptranf_s, and pass specmod vars through
! input variable sp of type(spec_vars)
! 2010-03-31 treadon - add double FFT capabilty (from Joe Sela, Mark Iredell)
!
! input arguments:
! wave - real (*) wave fields if idir>0
! grid - real (*) grid fields (starting at jb) if idir<0
! idir - integer transform flag
! (idir>0 for wave to grid, idir<0 for grid to wave)
!
! output arguments:
! wave - real (*) wave fields if idir<0
! grid - real (*) grid fields (starting at jb) if idir>0
!
! subprograms called:
! sptranf1 sptranf spectral transform
!
! remarks:
! This routine assumes that splib routine sptranf0 has been
! previously called. sptranf0 initializes arrays needed in
! the transforms.
!
! minimum grid dimensions for unaliased transforms to spectral:
! dimension linear quadratic
! ----------------------- --------- -------------
! imax 2*maxwv+2 3*maxwv/2*2+2
! jmax (idrt=4,iromb=0) 1*maxwv+1 3*maxwv/2+1
! jmax (idrt=4,iromb=1) 2*maxwv+1 5*maxwv/2+1
! jmax (idrt=0,iromb=0) 2*maxwv+3 3*maxwv/2*2+3
! jmax (idrt=0,iromb=1) 4*maxwv+3 5*maxwv/2*2+3
! jmax (idrt=256,iromb=0) 2*maxwv+1 3*maxwv/2*2+1
! jmax (idrt=256,iromb=1) 4*maxwv+1 5*maxwv/2*2+1
! ----------------------- --------- -------------
!
! attributes:
! language: fortran 77
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: zero
use general_specmod, only: spec_vars
implicit none
! Declare passed variables
type(spec_vars) ,intent(in ) :: sp_a
type(spec_vars) ,intent(in ) :: sp_b
integer(i_kind) ,intent(in ) :: idir
real(r_kind),dimension(sp_b%nc) ,intent(inout) :: wave
real(r_kind),dimension(sp_a%ijmax),intent(inout) :: grid
! Declare local variables
integer(i_kind) i,j,ii,jj,ijn,ijs,mp,ifact,kw,kwtop,imaxp2
real(r_kind),dimension(2*(sp_b%jcap+1)):: wtop
real(r_kind),dimension(sp_b%imax,2):: g
real(r_kind),dimension(sp_b%imax+2,2):: f
real(r_kind),dimension(50000+4*sp_b%imax):: tmpafft
real(r_kind),dimension(sp_b%ncd2)::tmppln
real(r_kind),dimension(sp_b%jcap+1)::tmpplntop
! Initialize local variables
mp=0
ifact = sp_b%imax/sp_a%imax
do i=1,2*(sp_b%jcap+1)
wtop(i)=zero
end do
kw=(sp_b%jcap+1)*((sp_b%iromb+1)*sp_b%jcap+2)
kwtop=2*(sp_b%jcap+1)
imaxp2=sp_b%imax+2
! Transform wave to grid
! ***NOTE***
! The FFT used in the transform below has been generalized to
! allow for projection of spectral coefficients onto double
! the desired number of longitudinal grid points. This
! approach is needed when transforming high wavenumber spectral
! coefficients to a coarser resoultion grid. For example, using
! splib to directly transform T878 spectral coefficients to an
! 1152 x 576 grid does not use Fourier modes above wavenumber 576.
! Joe Sela insightfully suggested doubling the number of points
! in the FFT and using every other point in the output grid.
! Mark Iredell coded up Joe's idea below.
if(idir>0) then
if(sp_b%precalc_pln)then
!$omp parallel do schedule(dynamic,1) private(j,tmpafft,f,g,i,ii,jj,ijn,ijs)
do j=sp_a%jb,sp_a%je
tmpafft(:)=sp_b%afft(:)
call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, &
sp_a%clat(j),sp_b%pln(1,j),sp_b%plntop(1,j),0,wave,wtop,f)
call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft)
! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, &
! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! sp_b%pln(1,j),sp_b%plntop(1,j),mp, &
! wave,wtop,g,idir)
jj = j-sp_a%jb
ijn = jj*sp_a%jn
ijs = jj*sp_a%js + sp_a%ioffset
do i=1,sp_a%imax
ii = ifact*(i-1)+1
grid(ijn+i)=g(ii,1)
grid(ijs+i)=g(ii,2)
enddo
enddo
else
!$omp parallel do schedule(dynamic,1) private(j,tmpafft,g,f,i,ii,jj,ijn,ijs,tmppln,tmpplntop)
do j=sp_a%jb,sp_a%je
tmpafft(:)=sp_b%afft(:)
call splegend(sp_b%iromb,sp_b%jcap,sp_b%slat(j),sp_b%clat(j),sp_b%eps,&
sp_b%epstop,tmppln,tmpplntop)
call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, &
sp_a%clat(j),tmppln,tmpplntop,0,wave,wtop,f)
call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft)
! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, &
! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! tmppln,tmpplntop,mp, &
! wave,wtop,g,idir)
jj = j-sp_a%jb
ijn = jj*sp_a%jn
ijs = jj*sp_a%js + sp_a%ioffset
do i=1,sp_a%imax
ii = ifact*(i-1)+1
grid(ijn+i)=g(ii,1)
grid(ijs+i)=g(ii,2)
enddo
enddo
end if
! Transform grid to wave
! ***WARNING***
! The code above has been generalized to handle the transform of
! high spectral representation fields to coarse physical space
! grids. The code below should not be used to transform coarse
! resolution grids to high spectral representation. Since this
! functionality is not yet needed in the GSI, the prudent action
! to take here is to print an ERROR message and terminate program
! execution if such a transform is requested.
else
if (sp_a%imax /= sp_b%imax) then
write(6,*)'GENERAL_SPTRANF_S: ***ERROR*** grid --> spectral transform NOT SAFE'
call stop2(330)
else
tmpafft(:)=sp_b%afft(:)
if(sp_a%precalc_pln)then
do j=sp_a%jb,sp_a%je
if(sp_a%wlat(j)>zero) then
jj = j-sp_a%jb
ijn = jj*sp_a%jn
ijs = jj*sp_a%js + sp_a%ioffset
do i=1,sp_a%imax
g(i,1)=grid(ijn+i)
g(i,2)=grid(ijs+i)
enddo
call spffte(sp_a%imax,imaxp2/2,sp_a%imax,2,f,g,-1,tmpafft)
call spanaly(sp_a%iromb,sp_a%jcap,sp_a%imax,imaxp2,kw,kwtop,1, &
sp_a%wlat(j),sp_a%clat(j),sp_a%pln(1,j),sp_a%plntop(1,j),0,f,wave,wtop)
! call sptranf1(sp_a%iromb,sp_a%jcap,sp_a%idrt,sp_a%imax,sp_a%jmax,j,j, &
! sp_a%eps,sp_a%epstop,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! sp_a%pln(1,j),sp_a%plntop(1,j),mp, &
! wave,wtop,g,idir)
endif
enddo
else
do j=sp_a%jb,sp_a%je
if(sp_a%wlat(j)>zero) then
jj = j-sp_a%jb
ijn = jj*sp_a%jn
ijs = jj*sp_a%js + sp_a%ioffset
do i=1,sp_a%imax
g(i,1)=grid(ijn+i)
g(i,2)=grid(ijs+i)
enddo
call splegend(sp_a%iromb,sp_a%jcap,sp_a%slat(j),sp_a%clat(j),sp_a%eps,&
sp_a%epstop,tmppln,tmpplntop)
call spffte(sp_a%imax,imaxp2/2,sp_a%imax,2,f,g,-1,tmpafft)
call spanaly(sp_a%iromb,sp_a%jcap,sp_a%imax,imaxp2,kw,kwtop,1, &
sp_a%wlat(j),sp_a%clat(j),sp_a%pln(1,j),sp_a%plntop(1,j),0,f,wave,wtop)
! call sptranf1(sp_a%iromb,sp_a%jcap,sp_a%idrt,sp_a%imax,sp_a%jmax,j,j, &
! sp_a%eps,sp_a%epstop,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! tmppln,tmpplntop,mp, &
! wave,wtop,g,idir)
end if
enddo
endif
endif
endif
end subroutine general_sptranf_s_b
subroutine general_sptez_v(sp,waved,wavez,gridu,gridv,idir)
!$$$ subprogram documentation block
! . . . .
! subprogram: sptez_v perform a simple vector spherical transform
! prgmmr: iredell org: np23 date: 1996-02-29
!
! abstract: this subprogram performs a spherical transform
! between spectral coefficients of divergence and curl
! and a vector field on a global cylindrical grid.
! the wave-space can be either triangular or rhomboidal.
! the grid-space can be either an equally-spaced grid
! (with or without pole points) or a gaussian grid.
! the wave field is in sequential 'ibm order'.
! the grid fiels is indexed east to west, then north to south.
! for more flexibility and efficiency, call sptran.
! subprogram can be called from a multiprocessing environment.
!
! This routine differs from splib routine sptezv in that
! 1) the calling list only contains the in/out arrays and
! flag for the direction in which to transform
! 2) it calls a version of sptranfv that does not invoke
! initialization routines on each entry
! 3) some generality built into the splib version is
! removed in the code below
!
! program history log:
! 1996-02-29 iredell
! 2004-08-23 treadon - adapt splib routine sptezv for gsi use
! 2007-04-25 errico - replace use of duplicate arguments in sptranf_v
! 2008-04-03 safford - rm unused vars
! 2010-02-18 parrish - copy to general_sptez_v, and pass specmod vars through
! input variable sp of type(spec_vars)
!
! input arguments:
! waved - real (2*mx) wave divergence field if idir>0
! where mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
! wavez - real (2*mx) wave vorticity field if idir>0
! where mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
! gridu - real (imax,jmax) grid u-wind (e->w,n->s) if idir<0
! gridv - real (imax,jmax) grid v-wind (e->w,n->s) if idir<0
! idir - integer transform flag
! (idir>0 for wave to grid, idir<0 for grid to wave)
!
! output arguments:
! waved - real (2*mx) wave divergence field if idir<0
! where mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
! wavez - real (2*mx) wave vorticity field if idir>0
! where mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
! gridu - real (imax,jmax) grid u-wind (e->w,n->s) if idir>0
! gridv - real (imax,jmax) grid v-wind (e->w,n->s) if idir>0
!
! subprograms called:
! sptranf_v - perform a vector spherical transform
!
! remarks: minimum grid dimensions for unaliased transforms to spectral:
! dimension linear quadratic
! ----------------------- --------- -------------
! imax 2*maxwv+2 3*maxwv/2*2+2
! jmax (idrt=4,iromb=0) 1*maxwv+1 3*maxwv/2+1
! jmax (idrt=4,iromb=1) 2*maxwv+1 5*maxwv/2+1
! jmax (idrt=0,iromb=0) 2*maxwv+3 3*maxwv/2*2+3
! jmax (idrt=0,iromb=1) 4*maxwv+3 5*maxwv/2*2+3
! jmax (idrt=256,iromb=0) 2*maxwv+1 3*maxwv/2*2+1
! jmax (idrt=256,iromb=1) 4*maxwv+1 5*maxwv/2*2+1
! ----------------------- --------- -------------
!
! attributes:
! language: fortran 77
!
!$$$
use kinds, only: r_kind,i_kind
use general_specmod, only: spec_vars
use constants, only: zero
implicit none
! Declare passed variables
type(spec_vars) ,intent(in ) :: sp
integer(i_kind) ,intent(in ) :: idir
real(r_kind),dimension(sp%nc) ,intent(inout) :: waved,wavez
real(r_kind),dimension(sp%ijmax),intent(inout) :: gridu,gridv
! Declare local variables
integer(i_kind) i
! Zero appropriate output array based on direction of transform
if (idir<0) then
do i=1,sp%nc
waved(i)=zero
wavez(i)=zero
end do
elseif (idir>0) then
do i=1,sp%ijmax
gridu(i)=zero
gridv(i)=zero
end do
endif
! Call spectral <--> grid transform
call general_sptranf_v(sp,sp,waved,wavez,gridu,gridv,idir)
end subroutine general_sptez_v
subroutine general_sptranf_v(sp_a,sp_b,waved,wavez,gridu,gridv,idir)
!$$$ subprogram documentation block
! . . . .
! subprogram: sptranf_v perform a vecor spherical transform
! prgmmr: iredell org: np23 date: 1996-02-29
!
! abstract: this subprogram performs a spherical transform
! between spectral coefficients of divergences and curls
! and vector fields on a global cylindrical grid.
! the wave-space can be either triangular or rhomboidal.
! the grid-space can be either an equally-spaced grid
! (with or without pole points) or a gaussian grid.
! the wave and grid fields may have general indexing,
! but each wave field is in sequential 'ibm order',
! i.e. with zonal wavenumber as the slower index.
! transforms are done in latitude pairs for efficiency;
! thus grid arrays for each hemisphere must be passed.
! if so requested, just a subset of the latitude pairs
! may be transformed in each invocation of the subprogram.
! the transforms are all multiprocessed over latitude except
! the transform from fourier to spectral is multiprocessed
! over zonal wavenumber to ensure reproducibility.
! transform several fields at a time to improve vectorization.
! subprogram can be called from a multiprocessing environment.
!
! This routine differs from splib routine sptranfv in that
! it does not call sptranf0 (an initialization routine).
!
! program history log:
! 1996-02-29 iredell
! 1998-12-15 iredell generic fft used
! 2004-08-23 treadon - adapt splib routine sptranfv for gsi use
! 2006-05-03 treadon - remove jc from specmod list since not used
! 2006-07-07 kleist - correct bug in indexing of j=1,2*ncd2 loop
! 2007-04-25 errico - replace use of duplicate arguments
! 2010-02-18 parrish - copy to general_sptranf_v, and pass specmod vars through
! input variable sp of type(spec_vars)
! 2010-03-31 treadon - add double FFT capabilty (from Joe Sela, Mark Iredell)
!
! input arguments:
! waved - real (*) wave divergence fields if idir>0
! wavez - real (*) wave vorticity fields if idir>0
! gridu - real (*) grid u-winds (starting at jb) if idir<0
! gridv - real (*) grid v-winds (starting at jb) if idir<0
! idir - integer transform flag
! (idir>0 for wave to grid, idir<0 for grid to wave)
!
! output arguments:
! waved - real (*) wave divergence fields if idir<0
! [waved=(d(gridu)/dlam+d(clat*gridv)/dphi)/(clat*rerth)]
! wavez - real (*) wave vorticity fields if idir<0
! [wavez=(d(gridv)/dlam-d(clat*gridu)/dphi)/(clat*rerth)]
! gridu - real (*) grid u-winds (starting at jb) if idir>0
! gridv - real (*) grid v-winds (starting at jb) if idir>0
!
! subprograms called:
! sptranf1 sptranf spectral transform
! spdz2uv compute winds from divergence and vorticity
! spuv2dz compute divergence and vorticity from winds
!
! remarks:
! This routine assumes that splib routine sptranf0 has been
! previously called. sptranf0 initializes arrays needed in
! the transforms.
!
! minimum grid dimensions for unaliased transforms to spectral:
! dimension linear quadratic
! ----------------------- --------- -------------
! imax 2*maxwv+2 3*maxwv/2*2+2
! jmax (idrt=4,iromb=0) 1*maxwv+1 3*maxwv/2+1
! jmax (idrt=4,iromb=1) 2*maxwv+1 5*maxwv/2+1
! jmax (idrt=0,iromb=0) 2*maxwv+3 3*maxwv/2*2+3
! jmax (idrt=0,iromb=1) 4*maxwv+3 5*maxwv/2*2+3
! jmax (idrt=256,iromb=0) 2*maxwv+1 3*maxwv/2*2+1
! jmax (idrt=256,iromb=1) 4*maxwv+1 5*maxwv/2*2+1
! ----------------------- --------- -------------
!
! attributes:
! language: fortran 77
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: zero
use general_specmod, only: spec_vars
implicit none
! Declare passed variables
type(spec_vars) ,intent(in ) :: sp_a
type(spec_vars) ,intent(in ) :: sp_b
integer(i_kind) ,intent(in ) :: idir
real(r_kind),dimension(sp_b%nc) ,intent(inout) :: waved,wavez
real(r_kind),dimension(sp_a%ijmax),intent(inout) :: gridu,gridv
! Declare local variables
integer(i_kind) i,j,ii,jj,ijn,ijs,ifact,kw,kwtop,imaxp2
integer(i_kind),dimension(2):: mp
real(r_kind),dimension(sp_b%ncd2*2,2):: w
real(r_kind),dimension(2*(sp_b%jcap+1),2):: wtop
real(r_kind),dimension(sp_b%imax,2):: g
real(r_kind),dimension(sp_b%imax+2,2):: f
real(r_kind),dimension(sp_b%ncd2*2,2):: winc
real(r_kind),dimension(50000+4*sp_b%imax):: tmpafft
real(r_kind),dimension(sp_b%ncd2)::tmppln
real(r_kind),dimension(sp_b%jcap+1)::tmpplntop
! Set parameters
mp=1
ifact = sp_b%imax/sp_a%imax
kw=(sp_b%jcap+1)*((sp_b%iromb+1)*sp_b%jcap+2)
kwtop=2*(sp_b%jcap+1)
imaxp2=sp_b%imax+2
! Transform wave to grid
! ***NOTE***
! The FFT used in the transform below has been generalized to
! allow for projection of spectral coefficients onto double
! the desired number of longitudinal grid points. This
! approach is needed when transforming high wavenumber spectral
! coefficients to a coarser resoultion grid. For example, using
! splib to directly transform T878 spectral coefficients to an
! 1152 x 576 grid does not use Fourier modes above wavenumber 576.
! Joe Sela insightfully suggested doubling the number of points
! in the FFT and using every other point in the output grid.
! Mark Iredell coded up Joe's idea below.
if(idir>0) then
call spdz2uv(sp_b%iromb,sp_b%jcap,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
waved,wavez, &
w(1,1),w(1,2),wtop(1,1),wtop(1,2))
if(sp_b%precalc_pln)then
!$omp parallel do schedule(dynamic,1) private(j,tmpafft,f,g,i,ii,jj,ijn,ijs)
do j=sp_a%jb,sp_a%je
tmpafft(:)=sp_b%afft(:)
call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, &
sp_a%clat(j),sp_b%pln(1,j),sp_b%plntop(1,j),mp,w(1,1),wtop(1,1),f)
call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft)
! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, &
! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! sp_b%pln(1,j),sp_b%plntop(1,j),mp, &
! w(1,1),wtop(1,1),g,idir)
jj = j-sp_a%jb
ijn = jj*sp_a%jn
ijs = jj*sp_a%js + sp_a%ioffset
do i=1,sp_a%imax
ii = ifact*(i-1)+1
gridu(ijn+i)=g(ii,1)
gridu(ijs+i)=g(ii,2)
enddo
call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, &
sp_a%clat(j),sp_b%pln(1,j),sp_b%plntop(1,j),mp,w(1,2),wtop(1,2),f)
call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft)
! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, &
! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! sp_b%pln(1,j),sp_b%plntop(1,j),mp, &
! w(1,2),wtop(1,2),g,idir)
do i=1,sp_a%imax
ii = ifact*(i-1)+1
gridv(ijn+i)=g(ii,1)
gridv(ijs+i)=g(ii,2)
enddo
enddo
else
!$omp parallel do schedule(dynamic,1) private(j,tmpafft,f,g,i,ii,jj,ijn,ijs,tmppln,tmpplntop)
do j=sp_a%jb,sp_a%je
tmpafft(:)=sp_b%afft(:)
call splegend(sp_b%iromb,sp_b%jcap,sp_b%slat(j),sp_b%clat(j),sp_b%eps,&
sp_b%epstop,tmppln,tmpplntop)
call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, &
sp_a%clat(j),tmppln,tmpplntop,mp,w(1,1),wtop(1,1),f)
call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft)
! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, &
! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! tmppln,tmpplntop,mp, &
! w(1,1),wtop(1,1),g,idir)
jj = j-sp_a%jb
ijn = jj*sp_a%jn
ijs = jj*sp_a%js + sp_a%ioffset
do i=1,sp_a%imax
ii = ifact*(i-1)+1
gridu(ijn+i)=g(ii,1)
gridu(ijs+i)=g(ii,2)
enddo
call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, &
sp_a%clat(j),tmppln,tmpplntop,mp,w(1,2),wtop(1,2),f)
call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft)
! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, &
! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! sp_b%pln(1,1),sp_b%plntop(1,1),mp, &
! w(1,2),wtop(1,2),g,idir)
do i=1,sp_a%imax
ii = ifact*(i-1)+1
gridv(ijn+i)=g(ii,1)
gridv(ijs+i)=g(ii,2)
enddo
enddo
end if
! Transform grid to wave
! ***WARNING***
! The code above has been generalized to handle the transform of
! high spectral representation fields to coarse physical space
! grids. The code below should not be used to transform coarse
! resolution grids to high spectral representation. Since this
! functionality is not yet needed in the GSI, the prudent action
! to take here is to print an ERROR message and terminate program
! execution if such a transform is requested.
else
if (sp_a%imax /= sp_b%imax) then
write(6,*)'GENERAL_SPTRANF_V ***ERROR*** grid --> spectral transform NOT SAFE'
call stop2(330)
else
w=zero
wtop=zero
tmpafft(:)=sp_b%afft(:)
if(sp_a%precalc_pln)then
do j=sp_a%jb,sp_a%je
if(sp_a%wlat(j)>zero) then
jj = j-sp_a%jb
ijn = jj*sp_a%jn
ijs = jj*sp_a%js + sp_a%ioffset
do i=1,sp_a%imax
g(i,1)=gridu(ijn+i)/sp_a%clat(j)**2
g(i,2)=gridu(ijs+i)/sp_a%clat(j)**2
enddo
call spffte(sp_a%imax,imaxp2/2,sp_a%imax,2,f,g,-1,tmpafft)
call spanaly(sp_a%iromb,sp_a%jcap,sp_a%imax,imaxp2,kw,kwtop,1, &
sp_a%wlat(j),sp_a%clat(j),sp_a%pln(1,j),sp_a%plntop(1,j),mp,f, &
w(1,1),wtop(1,1))
! call sptranf1(sp_a%iromb,sp_a%jcap,sp_a%idrt,sp_a%imax,sp_a%jmax,j,j, &
! sp_a%eps,sp_a%epstop,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! sp_a%pln(1,j),sp_a%plntop(1,j),mp, &
! w(1,1),wtop(1,1),g,idir)
do i=1,sp_a%imax
g(i,1)=gridv(ijn+i)/sp_a%clat(j)**2
g(i,2)=gridv(ijs+i)/sp_a%clat(j)**2
enddo
call spffte(sp_a%imax,imaxp2/2,sp_a%imax,2,f,g,-1,tmpafft)
call spanaly(sp_a%iromb,sp_a%jcap,sp_a%imax,imaxp2,kw,kwtop,1, &
sp_a%wlat(j),sp_a%clat(j),sp_a%pln(1,j),sp_a%plntop(1,j),mp,f, &
w(1,2),wtop(1,2))
! call sptranf1(sp_a%iromb,sp_a%jcap,sp_a%idrt,sp_a%imax,sp_a%jmax,j,j, &
! sp_a%eps,sp_a%epstop,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! sp_a%pln(1,j),sp_a%plntop(1,j),mp, &
! w(1,2),wtop(1,2),g,idir)
endif
enddo
else
do j=sp_a%jb,sp_a%je
if(sp_a%wlat(j)>zero) then
call splegend(sp_a%iromb,sp_a%jcap,sp_a%slat(j),sp_a%clat(j),sp_a%eps,&
sp_a%epstop,sp_a%pln(1,1),sp_a%plntop(1,1))
jj = j-sp_a%jb
ijn = jj*sp_a%jn
ijs = jj*sp_a%js + sp_a%ioffset
do i=1,sp_a%imax
g(i,1)=gridu(ijn+i)/sp_a%clat(j)**2
g(i,2)=gridu(ijs+i)/sp_a%clat(j)**2
enddo
call spffte(sp_a%imax,imaxp2/2,sp_a%imax,2,f,g,-1,tmpafft)
call spanaly(sp_a%iromb,sp_a%jcap,sp_a%imax,imaxp2,kw,kwtop,1, &
sp_a%wlat(j),sp_a%clat(j),sp_a%pln(1,1),sp_a%plntop(1,1),mp,f, &
w(1,1),wtop(1,1))
! call sptranf1(sp_a%iromb,sp_a%jcap,sp_a%idrt,sp_a%imax,sp_a%jmax,j,j, &
! sp_a%eps,sp_a%epstop,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! sp_a%pln(1,1),sp_a%plntop(1,1),mp, &
! w(1,1),wtop(1,1),g,idir)
do i=1,sp_a%imax
g(i,1)=gridv(ijn+i)/sp_a%clat(j)**2
g(i,2)=gridv(ijs+i)/sp_a%clat(j)**2
enddo
call spffte(sp_a%imax,imaxp2/2,sp_a%imax,2,f,g,-1,tmpafft)
call spanaly(sp_a%iromb,sp_a%jcap,sp_a%imax,imaxp2,kw,kwtop,1, &
sp_a%wlat(j),sp_a%clat(j),sp_a%pln(1,1),sp_a%plntop(1,1),mp,f, &
w(1,2),wtop(1,2))
! call sptranf1(sp_a%iromb,sp_a%jcap,sp_a%idrt,sp_a%imax,sp_a%jmax,j,j, &
! sp_a%eps,sp_a%epstop,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! sp_a%pln(1,1),sp_a%plntop(1,1),mp, &
! w(1,2),wtop(1,2),g,idir)
end if
enddo
endif
call spuv2dz(sp_a%iromb,sp_a%jcap,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, &
w(1,1),w(1,2),wtop(1,1),wtop(1,2), &
winc(1,1),winc(1,2))
do j=1,2*sp_a%ncd2
waved(j)=waved(j)+winc(j,1)
wavez(j)=wavez(j)+winc(j,2)
end do
endif
endif
end subroutine general_sptranf_v
subroutine general_sptranf_v_u(sp_a,sp_b,waved,wavez,gridu,gridv)
!$$$ subprogram documentation block
! . . . .
! subprogram: sptranf_v perform a vecor spherical transform
! prgmmr: iredell org: np23 date: 1996-02-29
!
! abstract: this subprogram performs a spherical transform
! between spectral coefficients of divergences and curls
! and vector fields on a global cylindrical grid.
! the wave-space can be either triangular or rhomboidal.
! the grid-space can be either an equally-spaced grid
! (with or without pole points) or a gaussian grid.
! the wave and grid fields may have general indexing,
! but each wave field is in sequential 'ibm order',
! i.e. with zonal wavenumber as the slower index.
! transforms are done in latitude pairs for efficiency;
! thus grid arrays for each hemisphere must be passed.
! if so requested, just a subset of the latitude pairs
! may be transformed in each invocation of the subprogram.
! the transforms are all multiprocessed over latitude except
! the transform from fourier to spectral is multiprocessed
! over zonal wavenumber to ensure reproducibility.
! transform several fields at a time to improve vectorization.
! subprogram can be called from a multiprocessing environment.
!
! This routine differs from splib routine sptranfv in that
! it does not call sptranf0 (an initialization routine).
!
! program history log:
! 1996-02-29 iredell
! 1998-12-15 iredell generic fft used
! 2004-08-23 treadon - adapt splib routine sptranfv for gsi use
! 2006-05-03 treadon - remove jc from specmod list since not used
! 2006-07-07 kleist - correct bug in indexing of j=1,2*ncd2 loop
! 2007-04-25 errico - replace use of duplicate arguments
! 2010-02-18 parrish - copy to general_sptranf_v, and pass specmod vars through
! input variable sp of type(spec_vars)
! 2010-03-31 treadon - add double FFT capabilty (from Joe Sela, Mark Iredell)
!
! input arguments:
! waved - real (*) wave divergence fields if idir>0
! wavez - real (*) wave vorticity fields if idir>0
! idir - integer transform flag (assumed > 0)
! (idir>0 for wave to grid, idir<0 for grid to wave)
!
! output arguments:
! gridu - real (*) grid u-winds (starting at jb) if idir>0
! gridv - real (*) grid v-winds (starting at jb) if idir>0
!
! subprograms called:
! sptranf1 sptranf spectral transform
! spdz2uv compute winds from divergence and vorticity
! spuv2dz compute divergence and vorticity from winds
!
! remarks:
! This routine assumes that splib routine sptranf0 has been
! previously called. sptranf0 initializes arrays needed in
! the transforms.
!
! minimum grid dimensions for unaliased transforms to spectral:
! dimension linear quadratic
! ----------------------- --------- -------------
! imax 2*maxwv+2 3*maxwv/2*2+2
! jmax (idrt=4,iromb=0) 1*maxwv+1 3*maxwv/2+1
! jmax (idrt=4,iromb=1) 2*maxwv+1 5*maxwv/2+1
! jmax (idrt=0,iromb=0) 2*maxwv+3 3*maxwv/2*2+3
! jmax (idrt=0,iromb=1) 4*maxwv+3 5*maxwv/2*2+3
! jmax (idrt=256,iromb=0) 2*maxwv+1 3*maxwv/2*2+1
! jmax (idrt=256,iromb=1) 4*maxwv+1 5*maxwv/2*2+1
! ----------------------- --------- -------------
!
! attributes:
! language: fortran 77
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: zero
use general_specmod, only: spec_vars
implicit none
! Declare passed variables
type(spec_vars) ,intent(in ) :: sp_a
type(spec_vars) ,intent(in ) :: sp_b
real(r_kind),dimension(sp_b%nc) ,intent(inout) :: waved,wavez
real(r_kind),dimension(sp_a%ijmax),intent(inout) :: gridu,gridv
! Declare local variables
integer(i_kind) i,j,ii,jj,ijn,ijs,ifact,kw,kwtop,imaxp2
integer(i_kind),dimension(2):: mp
real(r_kind),dimension(sp_b%ncd2*2,2):: w
real(r_kind),dimension(2*(sp_b%jcap+1),2):: wtop
real(r_kind),dimension(sp_b%imax,2):: g
real(r_kind),dimension(sp_b%imax+2,2):: f
real(r_kind),dimension(50000+4*sp_b%imax):: tmpafft
real(r_kind),dimension(sp_b%ncd2)::tmppln
real(r_kind),dimension(sp_b%jcap+1)::tmpplntop
! Set parameters
mp=1
ifact = sp_b%imax/sp_a%imax
kw=(sp_b%jcap+1)*((sp_b%iromb+1)*sp_b%jcap+2)
kwtop=2*(sp_b%jcap+1)
imaxp2=sp_b%imax+2
! Transform wave to grid
! ***NOTE***
! The FFT used in the transform below has been generalized to
! allow for projection of spectral coefficients onto double
! the desired number of longitudinal grid points. This
! approach is needed when transforming high wavenumber spectral
! coefficients to a coarser resoultion grid. For example, using
! splib to directly transform T878 spectral coefficients to an
! 1152 x 576 grid does not use Fourier modes above wavenumber 576.
! Joe Sela insightfully suggested doubling the number of points
! in the FFT and using every other point in the output grid.
! Mark Iredell coded up Joe's idea below.
call spdz2uv(sp_b%iromb,sp_b%jcap,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
waved,wavez, &
w(1,1),w(1,2),wtop(1,1),wtop(1,2))
if(sp_b%precalc_pln)then
!$omp parallel do schedule(dynamic,1) private(j,tmpafft,f,g,i,ii,jj,ijn,ijs)
do j=sp_a%jb,sp_a%je
tmpafft(:)=sp_b%afft(:)
call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, &
sp_a%clat(j),sp_b%pln(1,j),sp_b%plntop(1,j),mp,w(1,1),wtop(1,1),f)
call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft)
! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, &
! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! sp_b%pln(1,j),sp_b%plntop(1,j),mp, &
! w(1,1),wtop(1,1),g,1)
jj = j-sp_a%jb
ijn = jj*sp_a%jn
ijs = jj*sp_a%js + sp_a%ioffset
do i=1,sp_a%imax
ii = ifact*(i-1)+1
gridu(ijn+i)=g(ii,1)
gridu(ijs+i)=g(ii,2)
enddo
if(j == sp_a%jb)then
call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, &
sp_a%clat(j),sp_b%pln(1,j),sp_b%plntop(1,j),mp,w(1,2),wtop(1,2),f)
call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft)
! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, &
! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! sp_b%pln(1,j),sp_b%plntop(1,j),mp, &
! w(1,2),wtop(1,2),g,1)
do i=1,sp_a%imax
ii = ifact*(i-1)+1
gridv(ijn+i)=g(ii,1)
gridv(ijs+i)=g(ii,2)
enddo
end if
enddo
else
!$omp parallel do schedule(dynamic,1) private(j,tmpafft,f,g,i,ii,jj,ijn,ijs,tmppln,tmpplntop)
do j=sp_a%jb,sp_a%je
call splegend(sp_b%iromb,sp_b%jcap,sp_b%slat(j),sp_b%clat(j),sp_b%eps,&
sp_b%epstop,tmppln,tmpplntop)
tmpafft(:)=sp_b%afft(:)
call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, &
sp_a%clat(j),tmppln,tmpplntop,mp,w(1,1),wtop(1,1),f)
call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft)
! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, &
! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! tmppln,tmpplntop,mp, &
! w(1,1),wtop(1,1),g,1)
jj = j-sp_a%jb
ijn = jj*sp_a%jn
ijs = jj*sp_a%js + sp_a%ioffset
do i=1,sp_a%imax
ii = ifact*(i-1)+1
gridu(ijn+i)=g(ii,1)
gridu(ijs+i)=g(ii,2)
enddo
if(j == sp_a%jb)then
call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, &
sp_a%clat(j),tmppln,tmpplntop,mp,w(1,2),wtop(1,2),f)
call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft)
! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, &
! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! tmppln,tmpplntop,mp, &
! w(1,2),wtop(1,2),g,1)
do i=1,sp_a%imax
ii = ifact*(i-1)+1
gridv(ijn+i)=g(ii,1)
gridv(ijs+i)=g(ii,2)
enddo
endif
enddo
end if
end subroutine general_sptranf_v_u
subroutine general_sptranf_v_v(sp_a,sp_b,waved,wavez,gridu,gridv)
!$$$ subprogram documentation block
! . . . .
! subprogram: sptranf_v perform a vecor spherical transform
! prgmmr: iredell org: np23 date: 1996-02-29
!
! abstract: this subprogram performs a spherical transform
! between spectral coefficients of divergences and curls
! and vector fields on a global cylindrical grid.
! the wave-space can be either triangular or rhomboidal.
! the grid-space can be either an equally-spaced grid
! (with or without pole points) or a gaussian grid.
! the wave and grid fields may have general indexing,
! but each wave field is in sequential 'ibm order',
! i.e. with zonal wavenumber as the slower index.
! transforms are done in latitude pairs for efficiency;
! thus grid arrays for each hemisphere must be passed.
! if so requested, just a subset of the latitude pairs
! may be transformed in each invocation of the subprogram.
! the transforms are all multiprocessed over latitude except
! the transform from fourier to spectral is multiprocessed
! over zonal wavenumber to ensure reproducibility.
! transform several fields at a time to improve vectorization.
! subprogram can be called from a multiprocessing environment.
!
! This routine differs from splib routine sptranfv in that
! it does not call sptranf0 (an initialization routine).
!
! program history log:
! 1996-02-29 iredell
! 1998-12-15 iredell generic fft used
! 2004-08-23 treadon - adapt splib routine sptranfv for gsi use
! 2006-05-03 treadon - remove jc from specmod list since not used
! 2006-07-07 kleist - correct bug in indexing of j=1,2*ncd2 loop
! 2007-04-25 errico - replace use of duplicate arguments
! 2010-02-18 parrish - copy to general_sptranf_v, and pass specmod vars through
! input variable sp of type(spec_vars)
! 2010-03-31 treadon - add double FFT capabilty (from Joe Sela, Mark Iredell)
!
! input arguments:
! waved - real (*) wave divergence fields if idir>0
! wavez - real (*) wave vorticity fields if idir>0
! idir - integer transform flag (assumed > 0)
! (idir>0 for wave to grid, idir<0 for grid to wave)
!
! output arguments:
! gridu - real (*) grid u-winds (starting at jb) if idir>0
! gridv - real (*) grid v-winds (starting at jb) if idir>0
!
! subprograms called:
! sptranf1 sptranf spectral transform
! spdz2uv compute winds from divergence and vorticity
! spuv2dz compute divergence and vorticity from winds
!
! remarks:
! This routine assumes that splib routine sptranf0 has been
! previously called. sptranf0 initializes arrays needed in
! the transforms.
!
! minimum grid dimensions for unaliased transforms to spectral:
! dimension linear quadratic
! ----------------------- --------- -------------
! imax 2*maxwv+2 3*maxwv/2*2+2
! jmax (idrt=4,iromb=0) 1*maxwv+1 3*maxwv/2+1
! jmax (idrt=4,iromb=1) 2*maxwv+1 5*maxwv/2+1
! jmax (idrt=0,iromb=0) 2*maxwv+3 3*maxwv/2*2+3
! jmax (idrt=0,iromb=1) 4*maxwv+3 5*maxwv/2*2+3
! jmax (idrt=256,iromb=0) 2*maxwv+1 3*maxwv/2*2+1
! jmax (idrt=256,iromb=1) 4*maxwv+1 5*maxwv/2*2+1
! ----------------------- --------- -------------
!
! attributes:
! language: fortran 77
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: zero
use general_specmod, only: spec_vars
implicit none
! Declare passed variables
type(spec_vars) ,intent(in ) :: sp_a
type(spec_vars) ,intent(in ) :: sp_b
real(r_kind),dimension(sp_b%nc) ,intent(inout) :: waved,wavez
real(r_kind),dimension(sp_a%ijmax),intent(inout) :: gridu,gridv
! Declare local variables
integer(i_kind) i,j,ii,jj,ijn,ijs,ifact,kw,kwtop,imaxp2
integer(i_kind),dimension(2):: mp
real(r_kind),dimension(sp_b%ncd2*2,2):: w
real(r_kind),dimension(2*(sp_b%jcap+1),2):: wtop
real(r_kind),dimension(sp_b%imax,2):: g
real(r_kind),dimension(sp_b%imax+2,2):: f
real(r_kind),dimension(50000+4*sp_b%imax):: tmpafft
real(r_kind),dimension(sp_b%ncd2)::tmppln
real(r_kind),dimension(sp_b%jcap+1)::tmpplntop
! Set parameters
mp=1
ifact = sp_b%imax/sp_a%imax
kw=(sp_b%jcap+1)*((sp_b%iromb+1)*sp_b%jcap+2)
kwtop=2*(sp_b%jcap+1)
imaxp2=sp_b%imax+2
! Transform wave to grid
! ***NOTE***
! The FFT used in the transform below has been generalized to
! allow for projection of spectral coefficients onto double
! the desired number of longitudinal grid points. This
! approach is needed when transforming high wavenumber spectral
! coefficients to a coarser resoultion grid. For example, using
! splib to directly transform T878 spectral coefficients to an
! 1152 x 576 grid does not use Fourier modes above wavenumber 576.
! Joe Sela insightfully suggested doubling the number of points
! in the FFT and using every other point in the output grid.
! Mark Iredell coded up Joe's idea below.
call spdz2uv(sp_b%iromb,sp_b%jcap,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
waved,wavez, &
w(1,1),w(1,2),wtop(1,1),wtop(1,2))
if(sp_b%precalc_pln)then
!$omp parallel do schedule(dynamic,1) private(j,tmpafft,f,g,i,ii,jj,ijn,ijs)
do j=sp_a%jb,sp_a%je
jj = j-sp_a%jb
ijn = jj*sp_a%jn
ijs = jj*sp_a%js + sp_a%ioffset
tmpafft(:)=sp_b%afft(:)
if(j == sp_a%jb)then
call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, &
sp_a%clat(j),sp_b%pln(1,j),sp_b%plntop(1,j),mp,w(1,1),wtop(1,1),f)
call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft)
! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, &
! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! sp_b%pln(1,j),sp_b%plntop(1,j),mp, &
! w(1,1),wtop(1,1),g,1)
do i=1,sp_a%imax
ii = ifact*(i-1)+1
gridu(ijn+i)=g(ii,1)
gridu(ijs+i)=g(ii,2)
enddo
end if
call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, &
sp_a%clat(j),sp_b%pln(1,j),sp_b%plntop(1,j),mp,w(1,2),wtop(1,2),f)
call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft)
! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, &
! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! sp_b%pln(1,j),sp_b%plntop(1,j),mp, &
! w(1,2),wtop(1,2),g,1)
do i=1,sp_a%imax
ii = ifact*(i-1)+1
gridv(ijn+i)=g(ii,1)
gridv(ijs+i)=g(ii,2)
enddo
enddo
else
!$omp parallel do schedule(dynamic,1) private(j,tmpafft,f,g,i,ii,jj,ijn,ijs,tmppln,tmpplntop)
do j=sp_a%jb,sp_a%je
call splegend(sp_b%iromb,sp_b%jcap,sp_b%slat(j),sp_b%clat(j),sp_b%eps,&
sp_b%epstop,tmppln,tmpplntop)
jj = j-sp_a%jb
ijn = jj*sp_a%jn
ijs = jj*sp_a%js + sp_a%ioffset
tmpafft(:)=sp_b%afft(:)
if(j == sp_a%jb)then
call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, &
sp_a%clat(j),tmppln,tmpplntop,mp,w(1,1),wtop(1,1),f)
call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft)
! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, &
! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! tmppln,tmpplntop,mp, &
! w(1,1),wtop(1,1),g,1)
do i=1,sp_a%imax
ii = ifact*(i-1)+1
gridu(ijn+i)=g(ii,1)
gridu(ijs+i)=g(ii,2)
enddo
end if
call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, &
sp_a%clat(j),tmppln,tmpplntop,mp,w(1,2),wtop(1,2),f)
call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft)
! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, &
! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, &
! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), &
! tmppln,tmpplntop,mp, &
! w(1,2),wtop(1,2),g,1)
do i=1,sp_a%imax
ii = ifact*(i-1)+1
gridv(ijn+i)=g(ii,1)
gridv(ijs+i)=g(ii,2)
enddo
enddo
end if
end subroutine general_sptranf_v_v
subroutine general_sptez_s_b(sp_a,sp_b,wave,grid,idir)
!$$$ subprogram documentation block
! . . . .
! subprogram: sptez_s perform a simple scalar spherical transform
! prgmmr: iredell org: np23 date: 1996-02-29
!
! absract: this subprogram performs a spherical transform
! between spectral coefficients of a scalar quantity
! and a field on a global cylindrical grid.
! the wave-space can be either triangular or rhomboidal.
! the grid-space can be either an equally-spaced grid
! (with or without pole points) or a gaussian grid.
! the wave field is in sequential 'ibm order'.
! the grid field is indexed east to west, then north to south.
! for more flexibility and efficiency, call sptran.
! subprogram can be called from a multiprocessing environment.
!
! This routine differs from splib routine sptez in that
! 1) the calling list only contains the in/out arrays and
! flag for the direction in which to transform
! 2) it calls a version of sptranf that does not invoke
! initialization routines on each entry
! 3) some generality built into the splib version is
! removed in the code below
!
! program history log:
! 1996-02-29 iredell
! 2004-08-23 treadon - adapt splib routine sptez for gsi use
! 2007-04-25 errico - replace use of duplicate arguments in sptranf_s
! 2010-02-18 parrish - copy to general_sptez_s, and pass specmod vars through
! input variable sp of type(spec_vars)
! input arguments:
! wave - real (2*mx) wave field if idir>0
! where mx=(jcap+1)*((iromb+1)*jcap+2)/2
! grid - real (imax,jmax) grid field (e->w,n->s) if idir<0
! idir - integer transform flag
! (idir>0 for wave to grid, idir<0 for grid to wave)
!
! output arguments:
! wave - real (2*mx) wave field if idir<0
! where mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
! grid - real (imax,jmax) grid field (e->w,n->s) if idir>0
!
! subprograms called:
! sptranf_s - perform a scalar spherical transform
!
! remarks: minimum grid dimensions for unaliased transforms to spectral:
! dimension linear quadratic
! ----------------------- --------- -------------
! imax 2*maxwv+2 3*maxwv/2*2+2
! jmax (idrt=4,iromb=0) 1*maxwv+1 3*maxwv/2+1
! jmax (idrt=4,iromb=1) 2*maxwv+1 5*maxwv/2+1
! jmax (idrt=0,iromb=0) 2*maxwv+3 3*maxwv/2*2+3
! jmax (idrt=0,iromb=1) 4*maxwv+3 5*maxwv/2*2+3
! jmax (idrt=256,iromb=0) 2*maxwv+1 3*maxwv/2*2+1
! jmax (idrt=256,iromb=1) 4*maxwv+1 5*maxwv/2*2+1
! ----------------------- --------- -------------
!
! attributes:
! language: fortran 77
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: zero
use general_specmod, only: spec_vars
implicit none
! Declare passed variables
type(spec_vars) ,intent(in ) :: sp_a,sp_b
integer(i_kind) ,intent(in ) :: idir
real(r_kind),dimension(sp_b%nc) ,intent(inout) :: wave
real(r_kind),dimension(sp_a%ijmax),intent(inout) :: grid
! Declare local variables
integer(i_kind) i
! Zero appropriate output array based on direction of transform
if (idir<0) then
do i=1,sp_b%nc
wave(i)=zero
end do
elseif (idir>0) then
do i=1,sp_a%ijmax
grid(i)=zero
end do
endif
! Call spectral <--> grid transform
call general_sptranf_s_b(sp_a,sp_b,wave,grid,idir)
return
end subroutine general_sptez_s_b
subroutine general_sptez_v_b(sp_a,sp_b,waved,wavez,gridu,gridv,idir,iuvflag)
!$$$ subprogram documentation block
! . . . .
! subprogram: sptez_v perform a simple vector spherical transform
! prgmmr: iredell org: np23 date: 1996-02-29
!
! abstract: this subprogram performs a spherical transform
! between spectral coefficients of divergence and curl
! and a vector field on a global cylindrical grid.
! the wave-space can be either triangular or rhomboidal.
! the grid-space can be either an equally-spaced grid
! (with or without pole points) or a gaussian grid.
! the wave field is in sequential 'ibm order'.
! the grid fiels is indexed east to west, then north to south.
! for more flexibility and efficiency, call sptran.
! subprogram can be called from a multiprocessing environment.
!
! This routine differs from splib routine sptezv in that
! 1) the calling list only contains the in/out arrays and
! flag for the direction in which to transform
! 2) it calls a version of sptranfv that does not invoke
! initialization routines on each entry
! 3) some generality built into the splib version is
! removed in the code below
!
! program history log:
! 1996-02-29 iredell
! 2004-08-23 treadon - adapt splib routine sptezv for gsi use
! 2007-04-25 errico - replace use of duplicate arguments in sptranf_v
! 2008-04-03 safford - rm unused vars
! 2010-02-18 parrish - copy to general_sptez_v, and pass specmod vars through
! input variable sp of type(spec_vars)
!
! input arguments:
! waved - real (2*mx) wave divergence field if idir>0
! where mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
! wavez - real (2*mx) wave vorticity field if idir>0
! where mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
! gridu - real (imax,jmax) grid u-wind (e->w,n->s) if idir<0
! gridv - real (imax,jmax) grid v-wind (e->w,n->s) if idir<0
! idir - integer transform flag
! (idir>0 for wave to grid, idir<0 for grid to wave)
! iuvflag - == 0 do both u and v
! > 0 do u (and polar rows for v)
! < 0 do v (and polar rows for u)
!
! output arguments:
! waved - real (2*mx) wave divergence field if idir<0
! where mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
! wavez - real (2*mx) wave vorticity field if idir>0
! where mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
! gridu - real (imax,jmax) grid u-wind (e->w,n->s) if idir>0
! gridv - real (imax,jmax) grid v-wind (e->w,n->s) if idir>0
!
! subprograms called:
! sptranf_v - perform a vector spherical transform
!
! remarks: minimum grid dimensions for unaliased transforms to spectral:
! dimension linear quadratic
! ----------------------- --------- -------------
! imax 2*maxwv+2 3*maxwv/2*2+2
! jmax (idrt=4,iromb=0) 1*maxwv+1 3*maxwv/2+1
! jmax (idrt=4,iromb=1) 2*maxwv+1 5*maxwv/2+1
! jmax (idrt=0,iromb=0) 2*maxwv+3 3*maxwv/2*2+3
! jmax (idrt=0,iromb=1) 4*maxwv+3 5*maxwv/2*2+3
! jmax (idrt=256,iromb=0) 2*maxwv+1 3*maxwv/2*2+1
! jmax (idrt=256,iromb=1) 4*maxwv+1 5*maxwv/2*2+1
! ----------------------- --------- -------------
!
! attributes:
! language: fortran 77
!
!$$$
use kinds, only: r_kind,i_kind
use general_specmod, only: spec_vars
use constants, only: zero
implicit none
! Declare passed variables
type(spec_vars) ,intent(in ) :: sp_a,sp_b
integer(i_kind) ,intent(in ) :: idir,iuvflag
real(r_kind),dimension(sp_b%nc) ,intent(inout) :: waved,wavez
real(r_kind),dimension(sp_a%ijmax),intent(inout) :: gridu,gridv
! Declare local variables
integer(i_kind) i
! Zero appropriate output array based on direction of transform
if (idir<0) then
do i=1,sp_b%nc
waved(i)=zero
wavez(i)=zero
end do
elseif (idir>0) then
do i=1,sp_a%ijmax
gridu(i)=zero
gridv(i)=zero
end do
endif
if(iuvflag == 0)then
! Call spectral <--> grid transform
call general_sptranf_v(sp_a,sp_b,waved,wavez,gridu,gridv,idir)
else if(idir > 0) then
if(iuvflag > 0)then
! Call spectral <--> grid transform u only (and polar v)
call general_sptranf_v_u(sp_a,sp_b,waved,wavez,gridu,gridv)
else
! Call spectral <--> grid transform v only (and polar u)
call general_sptranf_v_v(sp_a,sp_b,waved,wavez,gridu,gridv)
end if
else
if(idir < 0)then
write(6,*)'GENERAL_SPTRANF_V_B ***ERROR*** grid --> spectral transform NOT SAFE'
call stop2(330)
end if
end if
end subroutine general_sptez_v_b