c
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c  .                                                             .
c  .                  copyright (c) 1998 by UCAR                 .
c  .                                                             .
c  .       University Corporation for Atmospheric Research       .
c  .                                                             .
c  .                      all rights reserved                    .
c  .                                                             .
c  .                                                             .
c  .                         SPHEREPACK                          .
c  .                                                             .
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c
c
c ... file trvsph.f
c
c     this file contains documentation and code for subroutine trvsph
c
c ... required files
c
c     sphcom.f, hrfft.f, gaqd.f, vhaec.f, vhsec.f, vhagc.f, vhsgc.f
c
c     subroutine trvsph (intl,igrida,nlona,nlata,iveca,ua,va,
c    +igridb,nlonb,nlatb,ivecb,ub,vb,wsave,lsave,lsvmin,work,
c    +lwork,lwkmin,dwork,ldwork,ier)
c
c *** author
c
c     John C. Adams (NCAR 1997), email: johnad@ncar.ucar.edu
c
c *** purpose
c
c     subroutine trvsph transfers vector data given in (ua,va) on a grid on
c     the full sphere to vector data in (ub,vb) on a grid on the full sphere.
c     the grids on which (ua,va) is given and (ub,vb) is generated can be
c     specified independently of each other (see the input arguments igrida,
c     igridb,iveca,ivecb).  ua and ub are the east longitudinal components of
c     the given and transformed vector fields.  va is either the latitudinal
c     or colatitudinal component of the given vector field (see iveca).
c     vb is either the latitudinal or colatitudinal component of the
c     transformed vector field (see ivecb).  for transferring scalar data
c     on the sphere, use subroutine trssph.
c
c *   notice that scalar and vector quantities are fundamentally different
c     on the sphere.  for example, vectors are discontinuous and multiple
c     valued at the poles.  scalars are continuous and single valued at the
c     poles. erroneous results would be produced if one attempted to transfer
c     vector fields between grids with subroutine trssph applied to each
c     component of the vector.
c
c *** underlying grid assumptions and a description
c
c     discussions with the ncar scd data support group and others indicate
c     there is no standard grid for storing observational or model generated
c     data on the sphere.  subroutine trvsph was designed to handle most
c     cases likely to be encountered when moving data from one grid format
c     to another.
c
c     the grid on which (ua,va) is given must be equally spaced in longitude
c     and either equally spaced or gaussian in latitude (or colatitude).
c     longitude, which can be either the first or second dimension of ua,va
c     subdivides [0,2pi) excluding the periodic point 2pi.  (co)latitude,
c     which can be the second or first dimension of ua,va, has south
c     to north or north to south orientation with increasing subscript
c     value in ua,va (see the argument igrida).
c
c     the grid on which ub,vb is generated must be equally spaced in longitude
c     and either equally spaced or gaussian in latitude (or colatitude).
c     longitude, which can be either the first or second dimension of ub,vb
c     subdivides [0,2pi) excluding the periodic point 2pi.  (co)latitude,
c     which can be the second or first dimension of ub,vb, has south
c     to north or north to south orientation with increasing subscript
c     value in db (see the argument igridb).
c
c     let nlon be either nlona or nlonb (the number of grid points in
c     longitude.  the longitude grid subdivides [0,2pi) into nlon spaced
c     points
c
c          (j-1)*2.*pi/nlon  (j=1,...,nlon).
c
c     it is not necessary to communicate to subroutine trvsph whether the
c     underlying grids are in latitude or colatitude.  it is only necessary
c     to communicate whether they run south to north or north to south with
c     increasing subscripts.  a brief discussion of latitude and colatitude
c     follows.  equally spaced latitude grids are assumed to subdivide
c     [-pi/2,pi/2] with the south pole at -pi/2 and north pole at pi/2.
c     equally spaced colatitude grids subdivide [0,pi] with the north pole
c     at 0 and south pole at pi.  equally spaced partitions on the sphere
c     include both poles.  gaussian latitude grids subdivide (-pi/2,pi/2)
c     and gaussian colatitude grids subdivide (0,pi).  gaussian grids do not
c     include the poles.  the gaussian grid points are uniquely determined by
c     the size of the partition.  they can be computed in colatitude in
c     (0,pi) (north to south) in double precision by the spherepack subroutine
c     gaqd.  let nlat be nlata or nlatb if either the ua,va or ub,vb grid is
c     gaussian.  let
c
c        north pole                             south pole
c        ----------                             ----------
c           0.0    <  cth(1) < ... < cth(nlat)  <   pi
c
c
c     be nlat gaussian colatitude points in the interval (0,pi) and let
c
c        south pole                        north pole
c        ----------                        ----------
c           -pi/2  < th(1) < ... < th(nlat) < pi/2
c
c     be nlat gaussian latitude points in the open interval (-pi/2,pi/2).
c     these are related by
c
c          th(i) = -pi/2 + cth(i)  (i=1,...,nlat)
c
c     if the (ua,va) or (ub,vb) grid is equally spaced in (co)latitude then
c
c          ctht(i) = (i-1)*pi/(nlat-1)
c                                               (i=1,...,nlat)
c          tht(i) = -pi/2 + (i-1)*pi/(nlat-1)
c
c     define the equally spaced (north to south) colatitude and (south to
c     north) latitude grids.
c
c *** method (simplified description)
c
c    (1)
c
c     the vector field (ua,va) is reformated to a vector field in mathematical
c     spherical coordinates using array transpositions, subscript reordering
c     and negation of va as necessary (see arguments igrida,iveca).
c
c     (2)
c
c     a vector harmonic analysis is performed on the result from (1)
c
c     (3)
c
c     a vector harmonic synthesis is performed on the (ub,vb) grid
c     using as many coefficients from (2) as possible (i.e., as
c     as is consistent with the size of the ub,vb grid).
c
c     (4)
c
c     the vector field generated in (3) is transformed from mathematical
c     spherical coordinates to the form flagged by ivecb and igridb in
c     (ub,vb) using array transpositions, subscript reordering and negation
c     as necessary
c
c
c *** advantages
c
c     the use of vector spherical harmonics to transfer vector data is
c     highly accurate and preserves properties of vectors on the sphere.
c     the method produces a weighted least squares fit to vector data in
c     which waves are resolved uniformly on the full sphere.  high frequencies
c     induced by closeness of grid points near the poles (due to computational
c     or observational errors) are smoothed.  the method is consistent with
c     methods used to generate vector data in numerical spectral models based
c     on spherical harmonics.  for more discussion of these and related issues,
c     see "on the spectral approximation of discrete scalar and vector
c     functions on the sphere," siam j. numer. anal., vol. 16, december 1979,
c     pp. 934-949, by paul swarztrauber.
c
c
c *** comment
c
c     on a nlon by nlat or nlat by nlon grid (gaussian or equally spaced)
c     spherical harmonic analysis generates and synthesis utilizes
c     min0(nlat,(nlon+2)/2)) by nlat coefficients.  consequently, for
c     ua,va and ub,vb,  if either
c
c             min0(nlatb,(nlonb+2)/2) < min0(nlata,(nlona+2)/2)
c
c     or if
c
c             nlatb < nlata
c
c     then all the coefficients generated by an analysis of ua,va cannot be
c     used in the synthesis which generates ub,vb.  in this case "information"
c     can be lost in generating ub,vb.  more precisely, information will be
c     lost if the analysis of ua,va yields nonzero coefficients which are
c     outside the coefficient bounds determined by the ub,vb grid. still
c     transference with vector spherical harmonics will yield results
c     consistent with grid resolution and is highly accurate.
c
c *** input arguments
c
c ... intl
c
c     an initialization argument which should be zero on an initial call to
c     trvsph.  intl should be one if trvsph is being recalled and
c
c          igrida,nlona,nlata,iveca,igridb,nlonb,nlatb,ivecb
c
c     have not changed from the previous call.  if any of these arguments have
c     changed intl=0 must be used to avoid undetectable errors.  when allowed,
c     calls with intl=1 bypass redundant computation and save time.  it can
c     be used when transferring multiple vector data sets with the same
c     underlying grids.
c
c ... igrida
c
c     an integer vector dimensioned two which identifies the underlying grid
c     on the full sphere for the given vector data (ua,va) as follows:
c
c     igrida(1)
c
c     = -1
c     if the latitude (or colatitude) grid for ua,va is an equally spaced
c     partition of [-pi/2,pi/2] ( or [0,pi] ) including the poles which
c     runs north to south with increasing subscript value
c
c     = +1
c     if the latitude (or colatitude) grid for ua,va is an equally spaced
c     partition of [-pi/2,pi/2] ( or [0,pi] ) including the poles which
c     runs south to north with increasing subscript value
c
c     = -2
c     if the latitude (or colatitude) grid for ua,va is a gaussian partition
c     of (-pi/2,pi/2) ( or (0,pi) ) excluding the poles which runs north
c     to south with increasing subscript value
c
c     = +2
c     if the latitude (or colatitude) grid for ua,va is a gaussian partition
c     of (-pi/2,pi/2) ( or (0,pi) ) excluding the poles which runs south
c     north with increasing subscript value
c
c     igrida(2)
c
c     = 0 if the underlying grid for ua,va is a nlona by nlata
c
c     = 1 if the underlying grid for ua,va is a nlata by nlona
c
c
c ... nlona
c
c     the number of longitude points on the uniform grid which partitions
c     [0,2pi) for the given vector (ua,va).  nlona is also the first or second
c     dimension of ua,va (see igrida(2)) in the program which calls trvsph.
c     nlona determines the grid increment in longitude as 2*pi/nlona. for
c     example nlona = 72 for a five degree grid.  nlona must be greater than
c     or equal to 4.  the efficiency of the computation is improved when
c     nlona is a product of small prime numbers
c
c ... nlata
c
c     the number of points in the latitude (or colatitude) grid for the
c     given vector (ua,va).  nlata is also the first or second dimension
c     of ua and va (see igrida(2)) in the program which calls trvsph.
c     if nlata is odd then the equator will be located at the (nlata+1)/2
c     gaussian grid point.  if nlata is even then the equator will be
c     located half way between the nlata/2 and nlata/2+1 grid points.
c
c ... iveca
c
c     if iveca=0 is input then va is the latitudinal component of the
c     given vector field. if iveca=1 then va is the colatitudinal
c     compoenent of the given vector field.  in either case, ua must
c     be the east longitudinal component of the given vector field.
c
c *** note:
c     igrida(1)=-1 or igrida(1)=-2, igrida(2)=1, and iveca=1 corresponds
c     to the "usual" mathematical spherical coordinate system required
c     by most of the drivers in spherepack2.  igrida(1)=1 or igrida(1)=2,
c     igrida(2)=0, and iveca=0 corresponds to the "usual" geophysical
c     spherical coordinate system.
c
c
c ... ua
c
c     ua is the east longitudinal component of the given vector field.
c     ua must be dimensioned nlona by nlata in the program calling trvsph if
c     igrida(2) = 0.  ua must be dimensioned nlata by nlona in the program
c     calling trvsph if igrida(2) = 1.  if ua is not properly dimensioned
c     and if the latitude (colatitude) values do not run south to north or
c     north to south as flagged by igrida(1) (this cannot be checked!) then
c     incorrect results will be produced.
c
c
c ... va
c
c     va is either the latitudinal or colatitudinal componenet of the
c     given vector field (see iveca).  va must be dimensioned nlona by
c     nlata in the program calling trvsph if igrida(2)=0.  va must be
c     dimensioned nlata by nlona in the program calling trvsph if
c     igrida(2)=1.  if va is not properly dimensioned or if the latitude
c     (colatitude) values do not run south to north or north to south
c     as flagged by igrida(1) (this cannot be checked!) then incorrect
c     results will be produced.
c
c ... igridb
c
c     an integer vector dimensioned two which identifies the underlying grid
c     on the full sphere for the transformed vector (ub,vb) as follows:
c
c     igridb(1)
c
c     = -1
c     if the latitude (or colatitude) grid for ub,vb is an equally spaced
c     partition of [-pi/2,pi/2] ( or [0,pi] ) including the poles which
c     north to south
c
c     = +1
c     if the latitude (or colatitude) grid for ub,vb is an equally spaced
c     partition of [-pi/2,pi/2] ( or [0,pi] ) including the poles which
c     south to north
c
c     = -2
c     if the latitude (or colatitude) grid for ub,vb is a gaussian partition
c     of (-pi/2,pi/2) ( or (0,pi) ) excluding the poles which runs north to
c     south
c
c     = +2
c     if the latitude (or colatitude) grid for ub,vb is a gaussian partition
c     of (-pi/2,pi/2) ( or (0,pi) ) excluding the poles which runs south to
c     north
c
c     igridb(2)
c
c     = 0 if the underlying grid for ub,vb is a nlonb by nlatb
c
c     = 1 if the underlying grid for ub,vb is a nlatb by nlonb
c
c
c ... nlonb
c
c     the number of longitude points on the uniform grid which partitions
c     [0,2pi) for the transformed vector (ub,vb).  nlonb is also the first or
c     second dimension of ub and vb (see igridb(2)) in the program which calls
c     trvsph.  nlonb determines the grid increment in longitude as 2*pi/nlonb.
c     for example nlonb = 72 for a five degree grid.  nlonb must be greater
c     than or equal to 4.  the efficiency of the computation is improved when
c     nlonb is a product of small prime numbers
c
c ... nlatb
c
c     the number of points in the latitude (or colatitude) grid for the
c     transformed vector (ub,vb).  nlatb is also the first or second dimension
c     of ub and vb (see igridb(2)) in the program which calls trvsph.
c     if nlatb is odd then the equator will be located at the (nlatb+1)/2
c     gaussian grid point.  if nlatb is even then the equator will be
c     located half way between the nlatb/2 and nlatb/2+1 grid points.
c
c ... ivecb
c
c     if ivecb=0 is input then vb is the latitudinal component of the
c     given vector field. if ivecb=1 then vb is the colatitudinal
c     compoenent of the given vector field.  in either case, ub must
c     be the east longitudinal component of the given vector field.
c
c *** note:
c     igridb(1)=-1 or igridb(1)=-2, igridb(2)=1, and ivecb=1 corresponds
c     to the "usual" mathematical spherical coordinate system required
c     by most of the drivers in spherepack2.  igridb(1)=1 or igridb(1)=2,
c     igridb(2)=0, and ivecb=0 corresponds to the "usual" geophysical
c     spherical coordinate system.
c
c ... wsave
c
c     a saved work space array that can be utilized repeatedly by trvsph
c     as long as the arguments nlata,nlona,nlatb,nlonb remain unchanged.
c     wsave is set by a intl=0 call to trvsph.  wsave must not be altered
c     when trvsph is being recalled with intl=1.
c
c ... lsave
c
c     the dimension of the work space wsave as it appears in the program
c     that calls trvsph.  the minimum required value of lsave for the
c     current set of input arguments is set in the output argument lsvmin.
c     it can be determined by calling trvsph with lsave=0 and printing lsvmin.
c
c          la1 = min0(nlata,(nlona+1)/2), la2 = (nlata+1)/2
c
c          lb1 = min0(nlatb,(nlonb+1)/2), lb2 = (nlatb+1)/2
c
c          lwa = 4*nlata*la2+3*max0(la1-2,0)*(2*nlata-la1-1)+la2+nlona+15
c
c          lwb = 4*nlatb*lb2+3*max0(lb1-2,0)*(2*nlatb-lb1-1)+nlonb+15
c
c      then
c
c          lsvmin = lwa + lwb
c
c      is the minimal required work space length of wsave
c
c
c ... work
c
c     a work array that does not have to be preserved
c
c ... lwork
c
c     the dimension of the array work as it appears in the program that
c     calls trvsph. the minimum required value of lwork for the current
c     set of input arguments is set in the output argument lwkmin.
c     it can be determined by calling trvsph with lwork=0 and printing
c     lwkmin.  an estimate for lwork follows.  let nlat = max0(nlata,nlatb),
c     nlon = max0(nlona,nlonb) and l1 = min0(nlat,(nlon+2)/2).  with these
c     these definitions, the quantity
c
c            2*nlat*(8*l1 + 4*nlon + 3)
c
c     will suffice as a length for the unsaved work space.  this formula
c     may overestimate the required minimum value for lwork.  the exact
c     minimum value can be predetermined by calling trvsph wtih lwork=0
c     and printout of lwkmin.
c
c ... dwork
c
c     a double precision work array that does not have to be preserved.
c
c ... ldwork
c
c     the length of dwork in the routine calling trvsph
c     Let
c
c       nlat = max0(nlata,nlatb)
c
c     ldwork must be at least 2*nlat*(nlat+1)+1
c
c
c *** output arguments
c
c
c ... ub
c
c     a two dimensional array that contains the east longitudinal component
c     of the transformed vector data.  ub
c     must be dimensioned nlonb by nlatb in the program calling trvsph if
c     igridb(2)=0.  ub must be dimensioned nlatb by nlonb in the program
c     calling trvsph if igridb(2)=1.  if ub is not properly dimensioned
c     and if the latitude (colatitude) values do not run south to north or
c     north to south as flagged by igrdb(1) (this cannot be checked!) then
c     incorrect results will be produced.
c
c
c ... vb
c
c     a two dimensional array that contains the latitudinal or colatitudinal
c     component of the transformed vector data (see ivecb).
c     vb must be dimensioned nlonb by nlatb in the program calling trvsph if
c     igridb(2)=0.  vb must be dimensioned nlatb by nlonb in the program
c     calling trvsph if igridb(2)=1.  if vb is not properly dimensioned
c     and if the latitude (colatitude) values do not run south to north or
c     north to south as flagged by igrdb(1) (this cannot be checked!) then
c     incorrect results will be produced.
c
c ... lsvmin
c
c     the minimum length of the saved work space in wsave.
c     lsvmin is computed even if lsave < lsvmin (ier = 10).
c
c ... lwkmin
c
c     the minimum length of the unsaved work space in work.
c     lwkmin is computed even if lwork < lwkmin (ier = 11).
c
c
c *** error argument
c
c ... ier = 0  if no errors are detected
c
c         = 1  if intl is not 0 or 1
c
c         = 2  if igrida(1) is not -1 or +1 or -2 or +2
c
c         = 3  if igrida(2) is not 0 or 1
c
c         = 4  if nlona is less than 4
c
c         = 5  if nlata is less than 3
c
c         = 6  if iveca is not 0 or 1
c
c         = 7  if igridb(1) is not -1 or +1 or -2 or +2
c
c         = 8  if igridb(2) is not 0 or 1
c
c         = 9  if nlonb is less than 4
c
c         =10  if nlatb is less than 3
c
c         =11  if ivecb is not 0 or 1
c
c         =12  if there is insufficient saved work space (lsave < lsvmin)
c
c         =13  if there is insufficient unsaved work space (lwork < lwkmin)
c
c         =14  indicates failure in an eigenvalue routine which computes
c              gaussian weights and points
c
c         =15  if ldwork is too small (insufficient double precision
c              unsaved work space)
c
c *****************************************************
c *****************************************************
c
c     end of argument description ... code follows
c
c *****************************************************
c *****************************************************
c