
    module spectral
!  SPECTRAL_MOD - module containing the interface to GFS spectral
!  transform routines.
!
!  SYNOPSIS
!     use spectral
!
!  DESCRIPTION
!  The routines here handle the conversion between
!  the model observation grid and the scaled physical variables used in
!  our data assimilation routines.  In particular, this version allows us
!  to assimilate wind and, if desired, relative humidity, ozone, and
!  other tracers.
!
!  The GFS model uses observations grids with leading dimension 1:LONFD,
!  where LONFD > LONF, the number of actual longitudes.  We therefore
!  only create and use arrays of leading dimension 1:LONF for the data
!  assimilation.  We pad with zeros on output and strip the extra rows
!  on input as needed.
!
!  ROUTINES IMPLEMENTED
!   subroutine spectral_code_init(fn,header,coeff)
!   subroutine to_physical_space(spect,physgrid)
!   pure subroutine to_standard_units(gfs,phys)
!   subroutine from_physical_space(physgrid,spect)
!   generic elemental function pres_to_energy(lndpa)
!   generic elemental function temp_to_energy(lndpa,lev or z)
!   generic elemental function wind_to_energy(lndpa,lev or z)
!   generic elemental function energy_to_pres(p)
!   generic elemental function energy_to_temp(t,lev or z)
!   generic elemental function energy_to_wind(w,lev or z)
!   pure subroutine from_standard_units(phys,gfs)
!   subroutine read_sigma_file(fn,header,coeff)
!   subroutine write_sigma_file(fn,header,coeff)
!   pure function rmsof(f)
!   subroutine output_sinlat(fn)
!   subroutine dump_header(fn,iunit,header)
!   subroutine report_rms(coeff)
!
!  MACHINE DEPENDENCIES
!  On Linux clusters, we may need to set the environment variable
!  FORT90L="-Wl,-T1" and/or FORT90L="-Wl,-T2" to effect byte-swapping
!  on logical units 1 (INSWAP) and 2 (OUTSWAP), which we use for binary data.
!  The reason is that the model is compiled with the big-endian library
!  for creating sigma files and such, and Intel-based clusters are
!  little-endian machines.
!
!  MODIFICATION HISTORY
!  09/23/05 - Added variable_to_energy conversion routines for consistent
!             application throughout the code.  Removed the
!             GET_CONVERSION_COEFFS interface in view of this.
!  08/27/05 - added compatibility to 2004 GFS as well as 2001 version;
!             added optional arg to DUMP_HEADER for benefit of pre/postanalysis.
!  06/29/05 - added routines to implement energy units and corrected comments
!             in a few locations.
!  06/19/05 - changes to header and to WRITE_SIGMA_FILE to be compatible
!             with the 2001 GFS.
!  06/13/05 - fixed errors in calls to inverse transforms; eliminated
!             tread/twrite in favor of read/write_sigma_file.
!  05/23/05 - eliminated global variables previously used by preanalysis.
!  04/11/05 - corrections to NH/SH indexing in the MRF grid.  (I think this
!             is finally correct and consistent with the 2001 version.)
!  03/04/05 - revised version by Eric.
!
!  REQUIRED PACKAGES
!
    use grid_indexing
    use gfs_init
    use gfs_transforms
    use constant  ! the rest of these from the support directory
    use resolution
    use gfs_ids
    use io
    use ferror
    implicit none
!
!  Single-precision reals in sigma files
!
    integer,parameter::IO4=R
!
! Northern and Southern Hemispheres
!
    integer,parameter::NH=1,SH=2
!
!  Parameters needed to convert to and from physical units.
!  Initialized in spectral_code_init.
!
    real(D),save,private::sinlat(LATG2),wgt(LATG2)
    real(D),save,private::wslev(LEVS)
    real(D),save,private::wq,wte  ! for energy unit conversion
    real(D),save,private::pln(LNTD,LATG2),plntop(JCAP1,LATG2)
    logical,save,private::initialized=.false.
!
!  Name of the sigma file, if any, from which the above coefficients
!  were initialized.
!
    character(MAX_FILE_NAME_LEN),save,private::init_sigma_fn=' '
!
!  This header is valid for the 2001 version of the GFS code that we use.
!  We supply default values for most components, as the only ones that
!  are likely to change in our simulations are the forecast hour
!  and the date.  These default values simplify the construction of the
!  data assimilation scripts because we don't need to pass around a reference
!  sigma file for the XI and XL coefficients.
!
    type sigma_file_header
       sequence   ! important for working around Lahey byte-swapping deficiency
       real(IO4)::fhour4=0.0  ! forecast hour
       integer::idate(4)=0  ! hh dd mm yyyy
       !
       !  XI defines values that are used to compute WSLEV and energy units
       !
       real(IO4)::xi(LEVS+1)=(/1.0_IO4,0.99000001_IO4,0.974179983_IO4, &
          0.954590023_IO4,0.930540025_IO4,0.901350021_IO4,0.866419971_IO4, &
          0.825269997_IO4,0.777729988_IO4,0.724009991_IO4,0.664820015_IO4, &
          0.601350009_IO4,0.535290003_IO4,0.468600005_IO4,0.403340012_IO4, &
          0.341369987_IO4,0.284209996_IO4,0.232859999_IO4,0.187830001_IO4, &
          0.149159998_IO4,0.116540000_IO4,0.0894500017_IO4,0.0672300011_IO4, &
          0.0491999984_IO4,0.0346900001_IO4,0.0230899993_IO4,0.0138600003_IO4, &
          0.0065700002_IO4,0.0_IO4/)
       !
       !  XL defines the sigma levels
       !
       real(IO4)::xl(LEVS)=(/0.994996905_IO4,0.982082129_IO4,0.96437335_IO4, &
          0.942546666_IO4,0.915917397_IO4,0.88384378_IO4,0.84578532_IO4, &
          0.801416039_IO4,0.750755608_IO4,0.694264889_IO4,0.632895589_IO4, &
          0.568091393_IO4,0.501681089_IO4,0.435679018_IO4,0.372047603_IO4, &
          0.312478721_IO4,0.258230925_IO4,0.210057467_IO4,0.168230236_IO4, &
          0.132610902_IO4,0.102782242_IO4,0.0781517029_IO4,0.0580480471_IO4, &
          0.041794762_IO4,0.0287503004_IO4,0.0183360856_IO4,0.0100562088_IO4, &
          0.00272623356_IO4/)
       real(IO4)::dummy(201-LEVP1-LEVS)=0.0
       real(IO4)::waves=JCAP
       real(IO4)::xlayers=LEVS
       real(IO4)::trun=1.0
       real(IO4)::order=2.0
       real(IO4)::realform=1.0
       real(IO4)::gencode=78.0
       real(IO4)::rlond=LONF
       real(IO4)::rlatd=LATG
       real(IO4)::rlonp=LONF
       real(IO4)::rlatp=LATG
       real(IO4)::rlonr=LONF
       real(IO4)::rlatr=LATG
       real(IO4)::tracers=2.0  ! relative humidity and ozone, presumably
       real(IO4)::subcen=0.0
       real(IO4)::ensemble(2)=0.0
       real(IO4)::ppid=0.0
       real(IO4)::slid=0.0
       real(IO4)::vcid=0.0
       real(IO4)::vmid=0.0
       real(IO4)::vtid=0.0
       real(IO4)::runid4=0.0
       real(IO4)::usrid4=0.0
       !
       !  The next three fields are broken out in the 2004 GFS; the 2001 GFS
       !  has DUMMY2 as length 21.
       !
       real(IO4)::pdryini4=0.0
       real(IO4)::xncld=0.0
       real(IO4)::xgf=0.0
       real(IO4)::dummy2(18)=0.0
    end type sigma_file_header
!
!  The following is needed for the messy workaround to Lahey runtime 
!  byte-swapping bugs.  SIGMA_HEADER_LENGTH is the number of default real
!  and integer elements in the file header.
!
    integer,parameter::SIGMA_HEADER_LENGTH=250
!
!  The naming conventions below are the same as those in the 2001 GFS code
!
    type spectral_coeffs
       sequence
       real(D)::gz(LNTD,2)   ! geopotential height
       real(D)::q(LNTD,2)  ! surface pressure
       real(D)::te(LNTD,2,LEVS)  ! temperature
       real(D)::di(LNTD,2,LEVS)  ! divergence
       real(D)::ze(LNTD,2,LEVS)  ! vorticity
       real(D)::rq(LNTD,2,LEVS)  ! relative humidity
       real(D)::oz(LNTD,2,LEVS,NTRAC)  ! ozone plus other tracers
    end type spectral_coeffs
!
!  Physical ("standard") unit to energy unit conversions
!
    interface pres_to_energy
       module procedure pres_to_energy_r,pres_to_energy_d
    end interface
    interface temp_to_energy
       module procedure temp_to_energy_lev,temp_to_energy_z
    end interface
    interface wind_to_energy
       module procedure wind_to_energy_lev,wind_to_energy_z
    end interface
!
!  Energy to physical ("standard") unit conversions
!
    interface energy_to_pres
       module procedure energy_to_pres_r,energy_to_pres_d
    end interface
    interface energy_to_temp
       module procedure energy_to_temp_lev,energy_to_temp_z
    end interface
    interface energy_to_wind
       module procedure energy_to_wind_lev,energy_to_wind_z
    end interface
!
    private::report_rms
    contains
!---------------------------------------------------------------------------
    subroutine spectral_code_init(file,header,coeff)
!  SPECTRAL_CODE_INIT - initialize the GFS spectral transforms.
!  FILE  : if present, then we read in the header from FILE; otherwise,
!     we use the default values supplied above.
!  HEADER := if present, we fill in the values using the sigma file header
!     from FILE or the default header values otherwise.
!  COEFF := if present, the spectral coefficients from the sigma file if present
!     or from the default header values otherwise.
!
    character(*),intent(in),optional::file
    type(sigma_file_header),intent(out),optional::header
    type(spectral_coeffs),intent(out),optional::coeff
!
!  Local variables
!
    type(sigma_file_header)::hdr  ! inherits default values from above
    real(D)::rmu(LEVS),rnu(LEVS),tov(LEVS)
    real(D)::colat1,coslat(LATG2),rcs2(LATG2)
    real(D)::rdel2(LEVS),slk(LEVS),rk,sl1,colrad(LATG2)
    real(D)::sl(LEVS),si(LEVP1),cl(LEVS),del(LEVS),qln(LNTD)
    integer::levshc,j
!
    if(present(file)) call read_sigma_file(file,hdr,coeff)
    if(present(header)) header=hdr
!
!  Initialize module variables
!
    wte=sqrt(CP/TR)  ! these values in module constant
    wq=SQRT(RD*TR)
    wslev=real(sqrt(hdr%xi(1:LEVS)-hdr%xi(2:LEVS+1)),kind(wslev))
!
    si=real(hdr%xi,kind(si))  ! any type conversion is deliberate
    sl=real(hdr%xl,kind(sl))
    call init_vert(cl,sl,del,si,tov,rmu,rnu,rdel2,rk,slk,sl1,levshc)
    call glats(colrad,wgt,rcs2,sinlat,coslat,colat1)  ! SINLAT is module var
    call gpln2i
    call ggozrm
    do j=1,LATG2  !  PLN and PLNTOP are module vars
       call pln2i(pln(1,j),qln,plntop(1,j),rcs2(j),colrad(j))
    enddo
    call toevenodd(LATG2,pln)
    initialized=.true.  ! module var
    if(present(file)) init_sigma_fn=file  ! for record-keeping purposes
    return
    end subroutine spectral_code_init
!---------------------------------------------------------------------------
    subroutine to_physical_space(spect,physgrid)
!  TO_PHYSICAL_SPACE - given the spectral coefficients in SPECT (and
!  other coeffs computed from SPECTRAL_CODE_INIT), compute the physical
!  grid.  The variables in the physical grid are in "standard" units
!  viz., surface pressure in ln(dPa);
!        temperature in Kelvin;
!        wind in meters per second.
!  and are indexed using the same conventions as the GFS.
!  SPECT : input spectral coefficients.
!  PHYSGRID := the physical grid (in units as listed above).  At any given
!     atmospheric level, the physical grid is indexed just as the GFS.
!     (However, the longitudinal dimension of the physical grid is LONF,
!     not LONFD, as in the GFS.)
!
    type(spectral_coeffs),intent(in)::spect
    real(GRIDKIND),intent(out)::physgrid(LONF,NAVARS,LATG,LEVS)
!
!  Local variables
!
    real(D)::coeffs(LNTD,2,NAVARS)
    real(D)::ftop(0:JCAP,2,2)
    real(D)::gfsgrid(LONFD,2,NAVARS,LATG2)  ! transform routines are d.p.
    integer::k,ivect(2)
!
    if(.not.initialized) &
       call fexit('to_physical_space: spectral_code_init must be called first')
!
!  Proceed by level to convert from spectral to physical units.
!  The surface pressure is transformed at each atmospheric level.
!
    do k=1,LEVS
       ftop=0.0  ! just in case; output arg in DZ2UV
       gfsgrid=0.0  ! just in case; output arg in FLN2FOUR
       call dz2uv(spect%di(:,:,k),spect%ze(:,:,k),coeffs(:,:,UWIND_TYPE), &
          coeffs(:,:,VWIND_TYPE),ftop(:,:,1),ftop(:,:,2))
       coeffs(:,:,PRES_TYPE)=spect%q
       coeffs(:,:,TEMP_TYPE)=spect%te(:,:,k)
       call toevenodd(2*NAVARS,coeffs)
       ivect=(/UWIND_TYPE,VWIND_TYPE/)  ! indexes GFSGRID in FLN2FOUR
       call fln2four(NAVARS,LATG2,NAVARS,coeffs,gfsgrid, &
          pln,2,ivect,ftop,plntop,LATG2,LONSPERLAT)
       call four2grid(NAVARS,LATG2,NAVARS,gfsgrid,LATG2,LONSPERLAT)
       call to_standard_units(gfsgrid,physgrid(:,:,:,k))
    enddo
    call to_energy_units(physgrid,1,LATG)
    return
    end subroutine to_physical_space
!---------------------------------------------------------------------------
    pure subroutine to_standard_units(gfs,phys)
!  TO_STANDARD_UNITS  - Utility routine to convert the transformed
!  coefficients in GFS to "standard" units in PHYS, as follows:
!     surface pressure in ln(dPa);
!     temperature in Kelvin;
!     wind components in meters per second (from pseudo-wind).
!  This routine does not alter the units of temperature or surface pressure.
!  It also handles Northern/Southern Hemisphere conversions and so must be
!  called with the entire physical grid at a given atmospheric level.
!  GFS : intermediate coefficients from GFS transform routines for
!    one atmospheric level.  This array is assumed to contain pseudo-wind,
!    and we convert to wind.
!  PHYS := converted values for one atmospheric leve.
!
    real(D),intent(in)::gfs(LONFD,2,NAVARS,LATG2)
    real(GRIDKIND),intent(out)::phys(LONF,NAVARS,LATG)
!
!  Local variables
!
    integer::j
    real(kind(gfs))::rcl
!
!  The data assimilation code uses the convention that latitude 1 is
!  the latitude closest to the South Pole and latitude LATG is the latitude
!  closest to the North Pole.  The GFS code uses two subscripts: (lat,k)
!  for this purpose.  K=1 corresponds to the Northern Hemisphere and K=2
!  to the Southern Hemisphere.  The latitude index runs from 1 (pole) to
!  LATG/2 (equator).  We make these adjustments here.
!  We do not use the extra LONF+1:LONFD rows of data in the GFS grid
!  in the data assimilation.
!
    do j=1,LATG2
       !
       !  South Pole to Equator
       !
       phys(:,PRES_TYPE,j)=real(gfs(1:LONF,SH,PRES_TYPE,j),kind(phys))
       phys(:,TEMP_TYPE,j)=real(gfs(1:LONF,SH,TEMP_TYPE,j),kind(phys))
       rcl=1.0/sqrt(1.0-sinlat(j)**2)  ! pseudo-wind -> wind
       phys(:,UWIND_TYPE,j)=real(gfs(1:LONF,SH,UWIND_TYPE,j)*rcl,kind(phys))
       phys(:,VWIND_TYPE,j)=real(gfs(1:LONF,SH,VWIND_TYPE,j)*rcl,kind(phys))
       !
       !  Equator to North Pole; otherwise identical conversions
       !
       phys(:,PRES_TYPE,j+LATG2)= &
          real(gfs(1:LONF,NH,PRES_TYPE,LATG2-j+1),kind(phys))
       phys(:,TEMP_TYPE,j+LATG2)= &
          real(gfs(1:LONF,NH,TEMP_TYPE,LATG2-j+1),kind(phys))
       rcl=1.0/sqrt(1.0-sinlat(LATG2-j+1)**2)
       phys(:,UWIND_TYPE,j+LATG2)= &
          real(gfs(1:LONF,NH,UWIND_TYPE,LATG2-j+1)*rcl,kind(phys))
       phys(:,VWIND_TYPE,j+LATG2)= &
          real(gfs(1:LONF,NH,VWIND_TYPE,LATG2-j+1)*rcl,kind(phys))
    enddo
    return
    end subroutine to_standard_units
!---------------------------------------------------------------------------
    elemental function pres_to_energy_r(lndpa)
!  PRES_TO_ENERGY_R - return the surface pressure (in ln(dPa)) in energy units.
!
    real(R)::pres_to_energy_r
    real(R),intent(in)::lndpa
!
    pres_to_energy_r=real(lndpa*wq,R)
    return
    end function pres_to_energy_r
!---------------------------------------------------------------------------
    elemental function pres_to_energy_d(lndpa)
!  PRES_TO_ENERGY_D - return the surface pressure (in ln(dPa)) in energy units.
!  Like PRES_TO_ENERGY_R, but the PRES argument - and the value of
!  the function - is real(D).
!
    real(D)::pres_to_energy_d
    real(D),intent(in)::lndpa
!
    pres_to_energy_d=lndpa*wq
    return
    end function pres_to_energy_d
!---------------------------------------------------------------------------
    elemental function temp_to_energy_lev(tkelvin,lev)
!  TEMP_TO_ENERGY_LEV - return the temperature (in Kelvin) in energy units
!  assuming that one is at level LEV in the model grid.
!  If LEV < 1 or LEV > 28, then return 0.
!
    real(R)::temp_to_energy_lev
    real(R),intent(in)::tkelvin
    integer,intent(in)::lev
!
    if(lev.lt.1.or.lev.gt.LEVS) then
       temp_to_energy_lev=0.0
    else
       temp_to_energy_lev=real(tkelvin*wte*wslev(lev),R)
    endif
    return
    end function temp_to_energy_lev
!---------------------------------------------------------------------------
    elemental function temp_to_energy_z(tkelvin,z)
!  TEMP_TO_ENERGY_Z - return the temperature (in Kelvin) in energy units
!  assuming that one is at level Z in the model grid.  Do the necessary
!  linear interpolation by model level.
!  If Z < 1, then return 0.
!
    real(R)::temp_to_energy_z
    real(R),intent(in)::tkelvin,z
!
    integer::iw
    real(D)::w
!
    iw=int(z)
    if(iw.ge.LEVS) then
       w=wslev(LEVS)
    else if(iw.ge.1) then
       w=wslev(iw)+(z-iw)*(wslev(iw+1)-wslev(iw))
    else
       w=0.0
    endif
    temp_to_energy_z=real(tkelvin*wte*w,R)
    return
    end function temp_to_energy_z
!---------------------------------------------------------------------------
    elemental function wind_to_energy_lev(wms,lev)
!  WIND_TO_ENERGY_LEV - return the wind (in m/s) in energy units
!  assuming that one is at level LEV in the model grid.
!  If LEV < 1 or LEV > 28, then return 0.
!
    real(R)::wind_to_energy_lev
    real(R),intent(in)::wms
    integer,intent(in)::lev
!
    if(lev.lt.1.or.lev.gt.LEVS) then
       wind_to_energy_lev=0.0
    else
       wind_to_energy_lev=real(wms*wslev(lev),R)
    endif
    return
    end function wind_to_energy_lev
!---------------------------------------------------------------------------
    elemental function wind_to_energy_z(wms,z)
!  WIND_TO_ENERGY_Z - return the wind (in m/s) in energy units
!  assuming that one is at level Z in the model grid.  Do the necessary
!  linear interpolation by model level.
!  If Z < 1, then return 0.
!
    real(R)::wind_to_energy_z
    real(R),intent(in)::wms,z
!
    integer::iw
    real(D)::w
!
    iw=int(z)
    if(iw.ge.LEVS) then
       w=wslev(LEVS)
    else if(iw.ge.1) then
       w=wslev(iw)+(z-iw)*(wslev(iw+1)-wslev(iw))
    else
       w=0.0
    endif
    wind_to_energy_z=real(wms*w,R)
    return
    end function wind_to_energy_z
!---------------------------------------------------------------------------
    pure subroutine to_energy_units(phys,south,north)
!  TO_ENERGY_UNITS  - Utility routine to convert from standard units to
!  energy units (as in the original version of the data assimilation code).
!  PHYS :=: on entry, contains "standard" units for the physical grid,
!     as output by TO_STANDARD_UNITS.  Overwritten on return by the 
!     corresponding energy units.
!  This routine is the inverse of FROM_ENERGY_UNITS.
!
    integer,intent(in)::south,north
    real(GRIDKIND),intent(inout)::phys(LONF,NAVARS,south:north,LEVS)
!
!  Local variables
!
    integer::lev
!
    phys(:,PRES_TYPE,:,:)=pres_to_energy(phys(:,PRES_TYPE,:,:))
    do lev=1,LEVS
       phys(:,TEMP_TYPE,:,lev)=temp_to_energy(phys(:,TEMP_TYPE,:,lev),lev)
       phys(:,UWIND_TYPE,:,lev)=wind_to_energy(phys(:,UWIND_TYPE,:,lev),lev)
       phys(:,VWIND_TYPE,:,lev)=wind_to_energy(phys(:,VWIND_TYPE,:,lev),lev)
    enddo
    return
    end subroutine to_energy_units
!---------------------------------------------------------------------------
    elemental function energy_to_pres_r(p)
!  PRES_TO_ENERGY_R - return the surface pressure (in energy units) in terms
!  of ln(dPa).
!
    real(R)::energy_to_pres_r
    real(R),intent(in)::p
!
    energy_to_pres_r=real(p/wq,R)
    return
    end function energy_to_pres_r
!---------------------------------------------------------------------------
    elemental function energy_to_pres_d(p)
!  PRES_TO_ENERGY_D - return the surface pressure (in energy units) in terms
!  of ln(dPa).  Double-precision argument and return value.
!
    real(D)::energy_to_pres_d
    real(D),intent(in)::p
!
    energy_to_pres_d=p/wq
    return
    end function energy_to_pres_d
!---------------------------------------------------------------------------
    elemental function energy_to_temp_lev(t,lev)
!  TEMP_TO_ENERGY_LEV - return the temperature (in energy units) as Kelvin
!  assuming that one is at level LEV in the model grid.
!  If LEV < 1 or LEV > 28, then return 0.
!
    real(R)::energy_to_temp_lev
    real(R),intent(in)::t
    integer,intent(in)::lev
!
    if(lev.lt.1.or.lev.gt.LEVS) then
       energy_to_temp_lev=0.0
    else
       energy_to_temp_lev=real(t/(wte*wslev(lev)),R)
    endif
    return
    end function energy_to_temp_lev
!---------------------------------------------------------------------------
    elemental function energy_to_temp_z(t,z)
!  TEMP_TO_ENERGY_Z - return the temperature (in Kelvin) in energy units
!  assuming that one is at level Z in the model grid.  Do the necessary
!  linear interpolation by model level.
!  If Z < 1, then return 0.
!
    real(R)::energy_to_temp_z
    real(R),intent(in)::t,z
!
    integer::iw
    real(D)::w
!
    iw=int(z)
    if(iw.ge.LEVS) then
       w=wslev(LEVS)
    else if(iw.ge.1) then
       w=wslev(iw)+(z-iw)*(wslev(iw+1)-wslev(iw))
    else
       w=0.0
    endif
    energy_to_temp_z=real(t/(wte*w),R)
    return
    end function energy_to_temp_z
!---------------------------------------------------------------------------
    elemental function energy_to_wind_lev(wind,lev)
!  WIND_TO_ENERGY_LEV - return the wind (in energy units) as m/s
!  assuming that one is at level LEV in the model grid.
!  If LEV < 1 or LEV > 28, then return 0.
!
    real(R)::energy_to_wind_lev
    real(R),intent(in)::wind
    integer,intent(in)::lev
!
    if(lev.lt.1.or.lev.gt.LEVS) then
       energy_to_wind_lev=0.0
    else
       energy_to_wind_lev=real(wind/wslev(lev),R)
    endif
    return
    end function energy_to_wind_lev
!---------------------------------------------------------------------------
    elemental function energy_to_wind_z(wind,z)
!  WIND_TO_ENERGY_Z - return the wind (in energy units) as m/s
!  assuming that one is at level Z in the model grid.  Do the necessary
!  linear interpolation by model level.
!  If Z < 1, then return 0.
!
    real(R)::energy_to_wind_z
    real(R),intent(in)::wind,z
!
    integer::iw
    real(D)::w
!
    iw=int(z)
    if(iw.ge.LEVS) then
       w=wslev(LEVS)
    else if(iw.ge.1) then
       w=wslev(iw)+(z-iw)*(wslev(iw+1)-wslev(iw))
    else
       w=0.0
    endif
    energy_to_wind_z=real(wind/w,R)
    return
    end function energy_to_wind_z
!---------------------------------------------------------------------------
    subroutine from_physical_space(physgrid,spect)
!  FROM_PHYSICAL_SPACE - get one ensemble solution from the array PHYS,
!  which defines the physical grid in energy units, to spectral
!  coordinates.
!  PHYSGRID :=: On input, physical grid coeffs as described above.
!    On output, the units are changed to "standard" units,
!    viz., surface pressure in ln(dPa);
!        temperature in Kelvin;
!        wind in meters/second.
!  SPECT  :=: The Q, TE, DI, and ZE components are set and the
!     other levels and components retain their input values (if any).
!
    real(GRIDKIND),intent(inout)::physgrid(LONF,NAVARS,LATG,LEVS)
    type(spectral_coeffs),intent(inout)::spect
!
!  Local variables
!
    real(D)::gfsgrid(LONFD,2,NAVARS,LATG2)  ! because spectral codes are d.p.
    real(D)::utop(0:JCAP,2),vtop(0:JCAP,2)
    real(D)::coeff(LNTD,2,NAVARS)
    real(D)::dudt(LONFD,2,LATG2),dvdt(LONFD,2,LATG2)
    integer::k
!
    if(.not.initialized) &
      call fexit('from_physical_space: spectral_code_init must be called first')
!
    call from_energy_units(physgrid,1,LATG)
!
!  Proceed one atmospheric level at a time.  The calling sequence here
!  is not exactly the inverse of that in TO_PHYSICAL_SPACE.
!  Among other differences, FOUR2FLN calls FROMEVENODD so that an
!  explicit call to the latter is unnecessary (contrast the call to
!  TOEVENODD in TO_PHYSICAL_SPACE).
!
    do k=1,LEVS
       call from_standard_units(physgrid(:,:,:,k),gfsgrid)
       call grid2four(NAVARS,LATG2,NAVARS,gfsgrid,LATG2,wgt,LONSPERLAT)
       dudt=gfsgrid(:,:,UWIND_TYPE,:)
       dvdt=gfsgrid(:,:,VWIND_TYPE,:)
       utop=0.0  ! essential; UTOP and VTOP are intent(inout) variables
       vtop=0.0
       !
       !  To turn on the bug that was present in the original version
       !  of the code, comment out the call to UVSUMS.
       !
       call uvsums(LATG2,dudt,dvdt,utop,vtop,plntop,LATG2,LONSPERLAT)
       coeff=0.0  ! just in case
       call four2fln(NAVARS,LATG2,NAVARS,gfsgrid,coeff,pln,LATG2,LONSPERLAT)
       call uv2dz(coeff(:,:,UWIND_TYPE),coeff(:,:,VWIND_TYPE), &
          utop,vtop,spect%di(:,:,k),spect%ze(:,:,k))  ! divergence and vorticity
       spect%q=coeff(:,:,PRES_TYPE)  ! surface pressure 
       spect%te(:,:,k)=coeff(:,:,TEMP_TYPE)   ! temperature
    enddo
    return
    end subroutine from_physical_space
!---------------------------------------------------------------------------
    pure subroutine from_energy_units(phys,south,north)
!  FROM_ENERGY_UNITS  - Utility routine to convert from energy units to
!  standard units.
!  PHYS :=: on entry, contains energy units for the physical grid,
!     as output by TO_ENERGY_UNITS.  Overwritten on return by the 
!     corresponding physical units.
!  This routine is the inverse of TO_ENERGY_UNITS.
!
    integer,intent(in)::south,north
    real(GRIDKIND),intent(inout)::phys(LONF,NAVARS,south:north,LEVS)
!
!  Local variables
!
    integer::lev
!
    phys(:,PRES_TYPE,:,:)=energy_to_pres(phys(:,PRES_TYPE,:,:))
    do lev=1,LEVS
       phys(:,TEMP_TYPE,:,lev)=energy_to_temp(phys(:,TEMP_TYPE,:,lev),lev)
       phys(:,UWIND_TYPE,:,lev)=energy_to_wind(phys(:,UWIND_TYPE,:,lev),lev)
       phys(:,VWIND_TYPE,:,lev)=energy_to_wind(phys(:,VWIND_TYPE,:,lev),lev)
    enddo
    return
    end subroutine from_energy_units
!---------------------------------------------------------------------------
    pure subroutine from_standard_units(phys,gfs)
!  FROM_STANDARD_UNITS - convert from standard units in PHYS,
!  viz., surface pressure in ln(hPa);
!        temperature in Kelvin;
!        wind in meters/second
!  to units that are compatible with the GFS transform routines,
!  namely, surface pressure in ln(dPa), temperature in Kelvin, and
!  pseudo-wind.
!
!  PHYS : one atmospheric level in standard physical units (as output
!     in TO_PHYSICAL_SPACE).
!  GFS := intermediate units suitable for direct computation with GFS
!     transform routines.
!
    real(GRIDKIND),intent(in)::phys(LONF,NAVARS,LATG)
    real(D),intent(out)::gfs(LONFD,2,NAVARS,LATG2)
!
!  Local variables
!
    real(kind(gfs))::rcl
    integer::j
!
!  WARNING: going from physical to spectral space is NOT the
!  inverse procedure as from spectral to physical space!  Note especially
!  the scaling from wind to pseudo-wind.
!  Alter the code below at your own risk; the whole data assimilation
!  scheme will break if this is messed up.
!
!  The data assimilation code uses the convention that latitude 1 is
!  the latitude closest to the North Pole and latitude LATG is the latitude
!  closest to the South Pole.  The GFS code uses two subscripts: (lat,k)
!  for this purpose.  K=1 corresponds to the Northern Hemisphere and K=2
!  to the Southern Hemisphere.  The latitude index runs from 1 (pole) to
!  LATG/2 (equator).  We make these adjustments here.
!  We do not use the extra LONF+1:LONFD rows of data in the GFS grid
!  in the data assimilation but simply zero them out for safety.
!
    gfs(LONF+1:LONFD,:,:,:)=0.0
    do j=1,LATG2
       !
       !  South Pole to Equator
       !
       gfs(1:LONF,SH,PRES_TYPE,j)=phys(:,PRES_TYPE,j)
       gfs(1:LONF,SH,TEMP_TYPE,j)=phys(:,TEMP_TYPE,j)
       rcl=real(1.0/sqrt(1.0-sinlat(j)**2),kind(gfs))  ! wind -> pseudo-wind
       gfs(1:LONF,SH,UWIND_TYPE,j)=phys(:,UWIND_TYPE,j)*rcl
       gfs(1:LONF,SH,VWIND_TYPE,j)=phys(:,VWIND_TYPE,j)*rcl
       !
       !  Equator to North Pole; otherwise identical conversions
       !
       gfs(1:LONF,NH,PRES_TYPE,LATG2-j+1)=phys(:,PRES_TYPE,j+LATG2)
       gfs(1:LONF,NH,TEMP_TYPE,LATG2-j+1)=phys(:,TEMP_TYPE,j+LATG2)
       rcl=real(1.0/sqrt(1.0-sinlat(LATG2-j+1)**2),kind(gfs))
       gfs(1:LONF,NH,UWIND_TYPE,LATG2-j+1)=phys(:,UWIND_TYPE,j+LATG2)*rcl
       gfs(1:LONF,NH,VWIND_TYPE,LATG2-j+1)=phys(:,VWIND_TYPE,j+LATG2)*rcl
    enddo
    return
    end subroutine from_standard_units
! -------------------------------------------------------------------------
    subroutine read_sigma_file(fn,header,coeff)
!  READ_SIGMA_FILE - read the sigma file, extracting header information from the
!  second record and, optionally, spectral model data from subsequent records.
!  FN : name of the sigma file to open.
!  IVERB : if positive, then display forecast hour and date on the log file.
!  HEADER := the sigma file header contents.
!  COEFF := if present, the specral coefficients in the file.  If absent,
!     then reading halts after the header.
!
!  Side effects:  opens and closes INSWAP.  Also prints RMS values to the
!  log file if IVERB > VERB_TERSE.
!
    character(*),intent(in)::fn
    type(sigma_file_header),intent(out)::header
    type(spectral_coeffs),intent(out),optional::coeff
!
!  Local variables
!
!---real(D),parameter::GA2=GRAV/(REARTH**2)
!---real(D)::snnp1(LNTD)
    real(IO4)::buff(2,LNTD)
    integer::j,k,n,ktr,istat
    integer::ntraci,ntozi,ntcwi,ncldi
!
!  On Lahey, the environment variable FORT90L does not enable byte-swapping
!  of components of derived types.  So, we create a dummy real(IO4) array
!  of length 250, which is the current size of the sigma file header record,
!  read this array (which the Lahey runtime byte swaps), then copy the
!  equivalenced header to the dummy argument (sigh).  This code will work
!  correctly on systems without this problem at the price of slight
!  inefficiency, provided that REAL(IO4) and INTEGER occupy the same
!  amount of storage (which they should, assuming IO4 is default real).
!
    real(IO4)::lahey_bug_workaround(SIGMA_HEADER_LENGTH)
    type(sigma_file_header)::byte_swapped_header
    equivalence(lahey_bug_workaround,byte_swapped_header)
!
!  Skip the first record in the sigma file and read the second, which contains
!  the sigma file header.
!
    open(INSWAP,file=fn,form='unformatted',status='old',iostat=istat)
    if(istat.ne.0) call fexit('read_sigma_file: cannot open '//trim(fn),2)
!
!  File header.  Lahey does not byte-swap derived types in binary i/o,
!  which is a PITA.  So we copy the header from an object that's equivalenced
!  to a sufficiently large array, which the Lahey runtime byte-swaps.
!  This code is portable, since the header is a sequence type, and will
!  work correctly (with a slight loss of efficiency) on implementations
!  that don't have Lahey's problem.
!
    read(INSWAP,end=911)
    read(INSWAP,end=911) lahey_bug_workaround  ! to get the byte-swapped header
    header=byte_swapped_header
!
!  Adjustments in inputs per 2004 GFS code (see PARA_TREAD in mpi_io.f)
!
    if(header%fhour4.lt.1.0) header%pdryini4=0.0
!
    call dump_header(fn,LOGFILE,header,verbosity())
!
    if(.not.present(coeff)) then  ! no sense in reading the remaining data
       close(INSWAP)
       return
    endif
!
!  Geopotential height
!
    read(INSWAP,end=911) buff
    coeff%gz=transpose(buff)
!
!  These adjustments are taken from the GFS code.
!  Note: If these adjustments are made, then they must be undone in
!  WRITE_SIGMA_FILE.
!
!---   k=0
!---   do j=0,JCAP
!---      do n=j,JCAP
!---         k=k+1
!---         snnp1(k)=n*(n+1)
!---      enddo
!---   enddo
!---   coeff%gz(:,1)=coeff%gz(:,1)*snnp1*GA2
!---   coeff%gz(:,2)=coeff%gz(:,2)*snnp1*GA2
!
!  Surface pressure
!
    read(INSWAP,end=911) buff
    coeff%q=transpose(buff)
!
!  Temperature
!
    do k=1,LEVS
       read(INSWAP,end=911) buff
       coeff%te(:,:,k)=transpose(buff)
    enddo
!
!  Divergence and vorticity
!
    do k=1,LEVS
       read(INSWAP,end=911) buff
       coeff%di(:,:,k)=transpose(buff)
       read(INSWAP,end=911) buff
       coeff%ze(:,:,k)=transpose(buff)
    enddo
!
!  Relative humidity
!
    do k=1,LEVS
       read(INSWAP,end=911) buff
       coeff%rq(:,:,k)=transpose(buff)
    enddo
!
!  Ozone and other tracers
!
    coeff%oz=0.0
    ntraci=nint(header%tracers-1.0)
    if(header%vtid.gt.0.0) then
       ntcwi=nint(header%vtid/10.0)
       ntozi=nint(header%vtid-ntcwi*10)
       ncldi=nint(header%xncld)
    else if(ntraci.eq.1) then
       ntozi=1
       ntcwi=0
       ncldi=0
    else if(ntraci.eq.2) then
       ntozi=1
       ntcwi=2
       ncldi=1
    else
       ntozi=1
       ntcwi=0
       ncldi=0
    endif
!
!  Determine which tracer we're reading; this logic from the 2004 GFS code
!
    do k=1,ntraci
       if(k.eq.ntozi) then
          ktr=NTOZ
       else if(k.ge.ntcwi.and.ncldi.gt.0) then
          do n=1,ncldi
             if(k.eq.ntcwi+n-1) ktr=NTCW+n-1
          enddo
       else
          ktr=0
       endif
       do j=1,LEVS
          read(INSWAP,end=911) buff
          if(ktr.gt.0.and.ktr.le.NTRAC) coeff%oz(:,:,j,ktr)=transpose(buff)
       enddo
    enddo
!
!  Compute and print (to log file) the global RMS of each dynamical
!  variable at level 5 if requested.
!
    if(verbosity().gt.VERB_TERSE) call report_rms(coeff,5,header%tracers)
    close(INSWAP)
    return
!
!  Handle errors
!
911 continue
    call finfo('read_sigma_file: unexpected EOF on read')
    call fexit('read_sigma_file: on Lahey, export FORT90L="-Wl,-T1"')
    end subroutine read_sigma_file
!---------------------------------------------------------------------------
    subroutine write_sigma_file(fn,header,coeff)
!  WRITE_SIGMA_FILE - write out the spectral coefficients in sigma file format.
!  Note: the sigma file format uses single-precision values in the 2001/2004
!  versions of the GFS; we denote the loss of precision with explicit
!  conversions.
!  FN  : name of the output file (any predecessor is overwritten).
!  HEADER : sigma file header
!  COEFF : spectral coefficients to write out.
!
!  Side effects:  opens and closes OUTSWAP.
!
    character(*),intent(in)::fn
    type(sigma_file_header),intent(in)::header
    type(spectral_coeffs),intent(in)::coeff
!
!  Local variables
!
    character(8),parameter::LAB(4)= &
       char(0)//char(0)//char(0)//char(0)//char(0)//char(0)//char(0)//char(0)
    real(IO4)::buff(2,LNTD)
!
!  On Lahey, the environment variable FORT90L does not enable byte-swapping
!  of components of derived types.  So, we create a dummy real(IO4) array
!  of length SIGMA_HEADER_LENGTH, which is the current size of the sigma file
!  header record, read this array (which the Lahey runtime byte swaps), then 
!  copy the bits to the dummy argument (sigh).  This code will work
!  correctly on systems without this problem at the price of slight
!  inefficiency, provided that REAL(IO4) and INTEGER occupy the same
!  amount of storage (which they should, assuming IO4 is default real).
!
    real(IO4)::lahey_bug_workaround(SIGMA_HEADER_LENGTH)
!
    integer::k,n
!
    call dump_header(fn,LOGFILE,header,verbosity())
    open(OUTSWAP,file=fn,form='unformatted',action='write',position='rewind')
!
!  Write a NUL lable at the beginning of the file
!
    write(OUTSWAP) LAB
!
!  Sigma file header.
!  Lahey does not byte-swap derived types in binary i/o,
!  which is a PITA.  So we copy the header to a sufficiently large array,
!  which the Lahey runtime byte-swaps.  This code is portable, since
!  the header is a sequence type, and will work correctly (with a slight
!  loss of efficiency) on implementations that don't have Lahey's problem.
!
    lahey_bug_workaround=transfer(header,lahey_bug_workaround)
    write(OUTSWAP) lahey_bug_workaround  ! to get the byte-swapped header
!
!  Geopotential height.  Note:  If any rescaling is done in READ_SIGMA_FILE,
!  then we must undo it here.
!
    buff=real(transpose(coeff%gz),kind(buff))  ! possible double -> single 
    write(OUTSWAP) buff
!
!  Surface pressure
!
    buff=real(transpose(coeff%q),kind(buff))
    write(OUTSWAP) buff
!
!  Temperature
!
    do k=1,LEVS
       buff=real(transpose(coeff%te(:,:,k)),kind(buff))
       write(OUTSWAP) buff
    enddo
!
!  Divergence and vorticity
!
    do k=1,LEVS
       buff=real(transpose(coeff%di(:,:,k)),kind(buff))
       write(OUTSWAP) buff
       buff=real(transpose(coeff%ze(:,:,k)),kind(buff))
       write(OUTSWAP) buff
    enddo
!
!  Relative humidity
!
    do k=1,LEVS
       buff=real(transpose(coeff%rq(:,:,k)),kind(buff))
       write(OUTSWAP) buff
    enddo
!
!  Ozone and other tracers
!
    do n=1,NTRAC
       do k=1,LEVS
          buff=real(transpose(coeff%oz(:,:,k,n)),kind(buff))
          write(OUTSWAP) buff
       enddo
    enddo
    close(OUTSWAP)
!
!  Verbose output about RMS at level 5 if requested.
!
    if(verbosity().gt.VERB_TERSE) call report_rms(coeff,5,header%tracers)
    return
    end subroutine write_sigma_file
!---------------------------------------------------------------------------
    pure function rmsof(f)
    real(D)::rmsof
    real(D),intent(in)::f(LNTD,2)
!
    rmsof=sqrt(0.5*sum(f**2))
    return
    end function rmsof
!---------------------------------------------------------------------------
    subroutine report_rms(coeff,lev,tracers)
    real(IO4),intent(in)::tracers
    integer,intent(in)::lev
    type(spectral_coeffs),intent(in)::coeff
!
    real(D)::rms(6)
    character(2)::label(6)
    integer::j,n
!
    if(lev.eq.1) then
       rms(1)=rmsof(coeff%q)  ! surf pres & geopotential height at level 1 only
       label(1)='q '
       rms(2)=rmsof(coeff%gz)
       label(2)='gz'
       write(LOGFILE,100) 1,(trim(label(j)),rms(j),j=1,2)
    endif
!
    rms(1)=rmsof(coeff%te(:,:,lev))
    rms(2)=rmsof(coeff%di(:,:,lev))
    rms(3)=rmsof(coeff%ze(:,:,lev))
    rms(4)=rmsof(coeff%rq(:,:,lev))
    label(1:4)=(/'te','di','ze','rq'/)
    n=4
    if(tracers.ge.2.0) then
       n=n+1
       rms(n)=rmsof(coeff%oz(:,:,lev,1))
       label(n)='oz'
    endif
    write(LOGFILE,100) lev,(trim(label(j)),rms(j),j=1,n)
100 format('lev',i3,' rms',6(1x,a,':',es10.3))
    return
    end subroutine report_rms
!---------------------------------------------------------------------------
    subroutine output_coeffs(fn)
!  OUTPUT_COEFFS - output the spectral arrays to the named file.
!
    character(*),intent(in)::fn
!
    integer::j
!
    wte=sqrt(CP/TR)  ! these values in module constant
    wq=SQRT(RD*TR)
    open(unit=OUTNOSWAP,file=fn,status='replace',action='write')
    write(OUTNOSWAP,100) (j,sinlat(j),j=1,LATG2)
100 format('sinlat',i4,es21.14)
    write(OUTNOSWAP,101) (j,wgt(j),j=1,LATG2)
101 format('wgt',i4,es21.14)
    close(OUTNOSWAP)
    return
    end subroutine output_coeffs
!---------------------------------------------------------------------------
    subroutine dump_header(fn,iunit,header,iverb)
!  DUMP_HEADER - dump file header info for sigma file FN to the indicated unit.
!
    character(*),intent(in)::fn
    integer,intent(in)::iunit
    type(sigma_file_header),intent(in)::header
    integer,intent(in),optional::iverb
!
    integer::ichat
!
    if(present(iverb)) then
       ichat=iverb
    else
       ichat=verbosity()
    endif
    if(ichat.eq.VERB_SILENT) return
!
!  Display forcast hour and date only unless greater verbosity is requested.
!
    if(iverb.gt.VERB_TERSE) then
       write(iunit,100) trim(fn)
       write(iunit,101) 'idate: ',header%idate,'fhour4: ',header%fhour4
    endif
100 format('header info for: ',a)
101 format(a,4i5,3x,a,f8.3)
!
    if(iverb.le.VERB_TERSE) return
!
    write(iunit,102) 'xi',header%xi(1:LEVS)
102 format(a,1x,7es11.3/(3x,7es11.3))
    write(iunit,102) 'xl',header%xl
!
    write(iunit,103) 'waves: ',header%waves,'xlayers: ',header%xlayers
103 format(2(a9,1x,es11.3,1x))
    write(iunit,103) 'trun: ',header%trun,'order: ',header%order
    write(iunit,103) 'gencode: ',header%gencode,'tracers: ',header%tracers
    write(iunit,103) 'rlond: ',header%rlond,'rlatd: ',header%rlatd
    write(iunit,103) 'rlonp: ',header%rlonp,'rlatp: ',header%rlatp
    write(iunit,103) 'rlonr: ',header%rlonr,'rlatr: ',header%rlatr
    write(iunit,103) 'subcen: ',header%subcen,'ppid: ',header%ppid
    write(iunit,103) 'slid: ',header%slid,'vcid: ',header%vcid
    write(iunit,103) 'vmid: ',header%vmid,'vtid: ',header%vtid
    write(iunit,103) 'runid4: ',header%runid4,'usrid4: ',header%usrid4
    write(iunit,103) 'pdryini4: ',header%pdryini4,'xncld: ',header%xncld
    write(iunit,103) 'xgf: ',header%xgf
    return
    end subroutine dump_header
!---------------------------------------------------------------------------
    end module spectral
