program format93
*_______________________________________________________________________

* date:  2/17/93

* author : Bill Ridgway (ridgway@climate.gsfc.nasa.gov, 301-286-9138)

* purpose: to covert gfdl, gsfc, or kiae data formats to fm93

* routines :

*    rgfdl : read data in gfdl format from unit 5
*    rgsfc : read data in gsfc format from unit 5
*    rkiae : read data in kiae format from unit 5
*    rfm93 : read data in fm93 format from unit 5
*    wfm93 : write data in fm93 format to unit 6

*_______________________________________________________________________
* proceedure: 1. ingest data at 5-10 cm-1 resolution into program memory
*             2. write data in newly specified format to output file
*_______________________________________________________________________
* parameters:
*    mlayer: maximum number of model layers
*    mlevel: maximum number of model levels
*    mbands: maximum number of spectral bands
*_______________________________________________________________________
* internal variables containing all model data
*    header(10): 10-line file header containing 80 characters per line
*    atmlbl    : 3-character atmospheric descriptor
*               'trp'  tropical
*               'mls'  mid-latitude summer
*               'mlw'  mid-latitude winter
*               'sas'  sub-arctic summer
*               'saw'  sub-arctic winter
*               '200'  200 (K) isothermal
*               '250'  250 (K) isothermal
*               '300'  300 (K) isothermal
*    iphase    : icrccm phase 1 or 2 (see report)
*    icase     : icrccm case number (see report)
*    nlayer    : actual number of model layers
*    nlevel    : actual number of model levels = nlayer+1
*    nbands    : actual number of spectral bands
*    psurf     : surface pressure (mb)
*    tsurf     : surface temperature (K)
*    ptrop     : tropopause pressure (mb)
*    flgctm    : flag to indicate h2o continuum is/is not (1/0) included
*    ppmco2    : co2 mixing ratio (ppm)
*    ppmch4    : ch4 mixing ratio (ppm)
*    ppmn2o    : n2o mixing ratio (ppm)
*    pthair    : total gas column for dry air (#air molecules/cm**2)
*    pthh2o    : total h2o column (#h2o molecules/cm**2)
*    pthco2    : co2 path (#co2 molecules/cm**2)
*    ptho3     : ozone path (#ozone molecules/cm**2)
*    pthch4    : ch4 path (#ch4 molecules/cm**2)
*    pthn2o    : n2o path (#n2o molecules/cm**2)
*    bandv1    : wavenumber lower limit of each band
*    bandv2    : wavenumber upper limit of each band
*    plevel    : pressure levels starting at top (mb)
*    tlevel    : temperature at pressure levels (K)
*    player    : pressure at mid-point of each layer (mb)
*    tlayer    : temperature characteristic of each layer (K)
*    alayer    : dry air mass of layer (#air molecules/cm**2)
*    wlayer    : water vapor volume mixing ratio of each layer (ppm)
*    olayer    : ozone mixing ratio of each layer (ppm)
*    fluxup    : radiative flux upward by level and band (w/m**2)
*    fluxdn    : radiative flux downward by level and band (w/m**2)
*    fluxnt    : net radiative flux upward by level and band (w/m**2)
*    heatrt    : layer heating rates indexed by layer and band (K/day)
*    tropup    : tropopause radiative flux upward by band (w/m**2)
*    tropdn    : tropopause radiative flux downward by band (w/m**2)
*    tropnt    : tropopause net radiative flux upward by band (w/m**2)
*_______________________________________________________________________

*  units and conversions:
*_______________________

*   dry air mass per layer =            [1/G]        * [1/Mair]  *  dP
*   (molecules/cm**2)                   [g/cm**2/mb]   (#air/g)     (mb)

*   H2O absorber amount    =  RH2O   *  [1/G]        * [1/Mh2o]  *  dP
*   (molecules/cm**2)         (g/g)     [g/cm**2/mb]   (#h2o/g)     (mb)

*   CO2 amount    =  10.e-6 * FCO2   *  [1/G]        * [1/Mair]  *  dP
*   (molecules/cm**2)         (ppmv)    [g/cm**2/mb]   (#air/g)     (mb)

*   O3  absorber amount    =  RO3    *  [1/G]        * [1/Mo3]   *  dP
*   (molecules/cm**2)         (g/g)     [g/cm**2/mb]   (#o3/g)      (mb)

* water:  RH2O is mass mixing ratio in (g/g)
* co2:    FCO2 is molecular mixing ratio in (ppmv)
* ozone:  RO3  is mass mixing ratio in (g/g)

*         [1/G]    = 1.02                g/cm**2/mb
*         [1/Mh2o] = 6.0256e+23/18.      #h2o/g
*         [1/Mo3]  = 6.0256e+23/48.      #o3/g
*         [1/Mair] = 6.0256e+23/28.97    #air/g   (dry air)
*         dP       = layer thickness     mb

* tropopause pressures for standard atmospheres:

*     trp=tropical                17 km or  93.7 mb
*     mls=mid-latitude summer     13 km or 179.0 mb
*     mlw=mid-latitude winter     10 km or 256.8 mb
*     sas=subarctic summer        10 km or 267.7 mb
*     saw=subarctic winter         9 km or 282.9 mb

*_______________________________________________________________________
* parameters defined in driver and passed to all subroutines:
      integer iounit,mlayer,mlevel,mbands
      parameter (mlayer=150,mlevel=151,mbands=600)

* variables passed to or created by subroutines:
      character*80 header(10)
      character*3  atmlbl
      integer      iphase,icase,nlayer,nlevel,nbands
      real         psurf,tsurf,ptrop,flgctm
      real         ppmco2,ppmch4,ppmn2o
      real         pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
      real         bandv1(0:mbands),bandv2(0:mbands)
      real         plevel(mlevel),tlevel(mlevel)
      real         player(mlayer),tlayer(mlayer)
      real         alayer(mlayer),wlayer(mlayer),olayer(mlayer)
      real         fluxup(mlevel,0:mbands)
      real         fluxdn(mlevel,0:mbands)
      real         fluxnt(mlevel,0:mbands)
      real         heatrt(mlayer,0:mbands)
      real         tropup(       0:mbands)
      real         tropdn(       0:mbands)
      real         tropnt(       0:mbands)
*_______________________________________________________________________

* identify file format as gsfc, gfdl, etc from first 4 characters
* ingest data by calling appropriate subroutine
* re-write data in 93 format using subroutine wfm93 (write format-93)

      character*4 first4

      iounit=5
      read(iounit,'(a)') first4
      rewind iounit

      if(first4.eq.'GFDL') then
         call rgfdl
     +        (iounit,mlayer,mlevel,mbands
     +        ,header,atmlbl,iphase,icase,nlayer,nlevel,nbands
     +        ,psurf,tsurf,ptrop,flgctm,ppmco2,ppmch4,ppmn2o
     +        ,pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
     +        ,bandv1,bandv2,plevel,tlevel
     +        ,player,tlayer,alayer,wlayer,olayer
     +        ,fluxup,fluxdn,fluxnt,heatrt
     +        ,tropup,tropdn,tropnt)

      else if(first4.eq.' ...') then
        call rgsfc
     +        (iounit,mlayer,mlevel,mbands
     +        ,header,atmlbl,iphase,icase,nlayer,nlevel,nbands
     +        ,psurf,tsurf,ptrop,flgctm,ppmco2,ppmch4,ppmn2o
     +        ,pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
     +        ,bandv1,bandv2,plevel,tlevel
     +        ,player,tlayer,alayer,wlayer,olayer
     +        ,fluxup,fluxdn,fluxnt,heatrt
     +        ,tropup,tropdn,tropnt)

      else if(first4.eq.'JUNE') then
        call rkiae
     +        (iounit,mlayer,mlevel,mbands
     +        ,header,atmlbl,iphase,icase,nlayer,nlevel,nbands
     +        ,psurf,tsurf,ptrop,flgctm,ppmco2,ppmch4,ppmn2o
     +        ,pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
     +        ,bandv1,bandv2,plevel,tlevel
     +        ,player,tlayer,alayer,wlayer,olayer
     +        ,fluxup,fluxdn,fluxnt,heatrt
     +        ,tropup,tropdn,tropnt)

      else if(first4.eq.'NASA') then
        call rfm93
     +        (iounit,mlayer,mlevel,mbands
     +        ,header,atmlbl,iphase,icase,nlayer,nlevel,nbands
     +        ,psurf,tsurf,ptrop,flgctm,ppmco2,ppmch4,ppmn2o
     +        ,pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
     +        ,bandv1,bandv2,plevel,tlevel
     +        ,player,tlayer,alayer,wlayer,olayer
     +        ,fluxup,fluxdn,fluxnt,heatrt
     +        ,tropup,tropdn,tropnt)

      else
        write(6,'(a)') 'data file format not recognized.'
        stop
      endif

      iounit=6
      call wfm93
     +        (iounit,mlayer,mlevel,mbands
     +        ,header,atmlbl,iphase,icase,nlayer,nlevel,nbands
     +        ,psurf,tsurf,ptrop,flgctm,ppmco2,ppmch4,ppmn2o
     +        ,pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
     +        ,bandv1,bandv2,plevel,tlevel
     +        ,player,tlayer,alayer,wlayer,olayer
     +        ,fluxup,fluxdn,fluxnt,heatrt
     +        ,tropup,tropdn,tropnt)

      stop
      end

      subroutine rgfdl
     +        (iounit,mlayer,mlevel,mbands
     +        ,header,atmlbl,iphase,icase,nlayer,nlevel,nbands
     +        ,psurf,tsurf,ptrop,flgctm,ppmco2,ppmch4,ppmn2o
     +        ,pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
     +        ,bandv1,bandv2,plevel,tlevel
     +        ,player,tlayer,alayer,wlayer,olayer
     +        ,fluxup,fluxdn,fluxnt,heatrt
     +        ,tropup,tropdn,tropnt)

*_______________________________________________________________________
* purpose: rgfdl - read gfdl data in earlier 1988 data format
*_______________________________________________________________________
* parameters defined in driver and passed to all subroutines:
      integer iounit,mlayer,mlevel,mbands

* variables passed to or created by subroutines:
      character*80 header(10)
      character*3  atmlbl
      integer      iphase,icase,nlayer,nlevel,nbands
      real         psurf,tsurf,ptrop,flgctm
      real         ppmco2,ppmch4,ppmn2o
      real         pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
      real         bandv1(0:mbands),bandv2(0:mbands)
      real         plevel(mlevel),tlevel(mlevel)
      real         player(mlayer),tlayer(mlayer)
      real         alayer(mlayer),wlayer(mlayer),olayer(mlayer)
      real         fluxup(mlevel,0:mbands)
      real         fluxdn(mlevel,0:mbands)
      real         fluxnt(mlevel,0:mbands)
      real         heatrt(mlayer,0:mbands)
      real         tropup(       0:mbands)
      real         tropdn(       0:mbands)
      real         tropnt(       0:mbands)
*_______________________________________________________________________

* icrccm phase, case #, atmosphere, etc. entered from separate file

      open(1,file='temp_header')
      read(1,*) iphase
      read(1,*) icase
      read(1,'(a3)') atmlbl
      read(1,*) flgctm
      read(1,*) ppmco2
      read(1,*) ppmch4
      read(1,*) ppmn2o
      close(1)

* blank out header initially

      do k=1,10
         write(header(k),'(80a1)') (' ',j=1,80)
      end do

* read 3 lines of header info from gfdl file, number of levels and bands

      read( iounit,'(a80)') header(1)
      read( iounit,'(a80)') header(3)
      read( iounit,'(a80)') header(5)
      read( iounit,*) nlevel,nbands
      nlayer=nlevel-1

* read spectral sums (k=0th elements)

      do k=0,0
         read( iounit,*) bandv1(k),bandv2(k)
         read( iounit,*) (fluxnt(l,k), l = 1,nlevel )
         read( iounit,*) (fluxup(l,k), l = 1,nlevel )
         read( iounit,*) (fluxdn(l,k), l = 1,nlevel )
         read( iounit,*) (heatrt(l,k), l = 1,nlayer )
         read( iounit,*) ptrop,tropnt(k),tropup(k),tropdn(k)
      end do

* read atmospheric profile information

      read( iounit,*) ( plevel(l), l = 1,nlevel )
      read( iounit,*) ( tlevel(l), l = 1,nlevel )
      read( iounit,*) ( player(l), l = 1,nlayer )
      read( iounit,*) ( tlayer(l), l = 1,nlayer )
      read( iounit,*) ( wlayer(l), l = 1,nlayer )
      read( iounit,*) ( olayer(l), l = 1,nlayer )

* ingest flux and heating rate data band-by-band

      do k=1,nbands
         read( iounit,*) bandv1(k),bandv2(k)
         read( iounit,*) (fluxnt(l,k), l = 1,nlevel )
         read( iounit,*) (fluxup(l,k), l = 1,nlevel )
         read( iounit,*) (fluxdn(l,k), l = 1,nlevel )
         read( iounit,*) (heatrt(l,k), l = 1,nlayer )
         read( iounit,*) ptrop,tropnt(k),tropup(k),tropdn(k)
      end do

* set up gfdl scaling factors for h2o and o3

      if(iphase.eq.1.and.icase.eq.17) then
        wfactor=0.75
      else if(iphase.eq.1.and.icase.eq.18) then
        wfactor=0.75
      else if(iphase.eq.1.and.icase.eq.21) then
        wfactor=1.25
      else if(iphase.eq.1.and.icase.eq.22) then
        wfactor=1.25
      else if(iphase.eq.2.and.icase.eq.15) then
        wfactor=0.75
      else if(iphase.eq.2.and.icase.eq.16) then
        wfactor=0.75
      else if(iphase.eq.2.and.icase.eq.17) then
        wfactor=1.25
      else if(iphase.eq.2.and.icase.eq.18) then
        wfactor=1.25
      else
        wfactor=1.
      end if

      if(icase.eq.24.and.header(3)(30:40).eq.'0.75 X CONC') then
        ofactor=0.75
      else if(icase.eq.24.and.header(3)(30:40).eq.'1.25 X CONC') then
        ofactor=1.25
      else
        ofactor=1.
      end if

* re-define units as ppm, molecules/cm**2

      do j=1,nlayer
         alayer(j) = 1.02 * (6.0256e+23/28.97) * (plevel(j+1)-plevel(j))
         wlayer(j) = wfactor * 1.e+6 * wlayer(j) * 28.97/18.00
         olayer(j) = ofactor * 1.e+6 * olayer(j) * 28.97/48.00
      end do

* sum column amounts for dry air and 5 gases

      pthair = 0.
      pthh2o = 0.
      ptho3  = 0.

      do j=1,nlayer
         pthair    = pthair + alayer(j)
         pthh2o    = pthh2o + 1.e-6 * wlayer(j) * alayer(j)
         ptho3     = ptho3  + 1.e-6 * olayer(j) * alayer(j)
      end do

      pthco2 = 1.e-6 * ppmco2 * pthair
      pthch4 = 1.e-6 * ppmch4 * pthair
      pthn2o = 1.e-6 * ppmn2o * pthair

* miscellaneous

      psurf = plevel(nlevel)
      tsurf = tlevel(nlevel)

      return

      end

      subroutine rkiae
     +        (iounit,mlayer,mlevel,mbands
     +        ,header,atmlbl,iphase,icase,nlayer,nlevel,nbands
     +        ,psurf,tsurf,ptrop,flgctm,ppmco2,ppmch4,ppmn2o
     +        ,pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
     +        ,bandv1,bandv2,plevel,tlevel
     +        ,player,tlayer,alayer,wlayer,olayer
     +        ,fluxup,fluxdn,fluxnt,heatrt
     +        ,tropup,tropdn,tropnt)

*_______________________________________________________________________
* purpose: rkiae - read kiae icrccm data in earlier 1988 data format
*_______________________________________________________________________
* parameters defined in driver and passed to all subroutines:
      integer iounit,mlayer,mlevel,mbands

* variables passed to or created by subroutines:
      character*80 header(10)
      character*3  atmlbl
      integer      iphase,icase,nlayer,nlevel,nbands
      real         psurf,tsurf,ptrop,flgctm
      real         ppmco2,ppmch4,ppmn2o
      real         pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
      real         bandv1(0:mbands),bandv2(0:mbands)
      real         plevel(mlevel),tlevel(mlevel)
      real         player(mlayer),tlayer(mlayer)
      real         alayer(mlayer),wlayer(mlayer),olayer(mlayer)
      real         fluxup(mlevel,0:mbands)
      real         fluxdn(mlevel,0:mbands)
      real         fluxnt(mlevel,0:mbands)
      real         heatrt(mlayer,0:mbands)
      real         tropup(       0:mbands)
      real         tropdn(       0:mbands)
      real         tropnt(       0:mbands)
*_______________________________________________________________________

* icrccm phase, case #, atmosphere, etc. entered from separate file

      open(1,file='temp_header')
      read(1,*) iphase
      read(1,*) icase
      read(1,'(a3)') atmlbl
      read(1,*) flgctm
      read(1,*) ppmco2
      read(1,*) ppmch4
      read(1,*) ppmn2o
      close(1)

* blank out header initially

      do k=1,10
         write(header(k),'(80a1)') (' ',j=1,80)
      end do

* read 3 lines of header info from gfdl file, number of levels and bands

      read( iounit,'(a80)') header(1)
      read( iounit,'(a80)') header(3)
      write(header(6),'(a)')
     + '  *************************************************************'
      write(header(7),'(a)')
     + '  PLEASE NOTE: WATER VAPOR AND OZONE PROFILES ARE NOT PROVIDED.'
      write(header(8),'(a)')
     + '               HEATING RATES ARE REPORTED AT FLUX LEVELS.      '
      write(header(9),'(a)')
     + '  *************************************************************'
      read( iounit,*) nlevel,nbands
      nlayer=nlevel-1

* read spectral sums (k=0th elements)

      do k=0,0
         read( iounit,*) bandv1(k),bandv2(k)
         read( iounit,*) (fluxnt(l,k), l = 1,nlevel )
         read( iounit,*) (fluxup(l,k), l = 1,nlevel )
         read( iounit,*) (fluxdn(l,k), l = 1,nlevel )
         read( iounit,*) (heatrt(l,k), l = 1,nlayer )
         read( iounit,*) tropnt(k),tropup(k),tropdn(k)
      end do

* read atmospheric profile information

      read( iounit,'(1x,7e11.5)') ( plevel(l), l = 1,nlevel )
      read( iounit,'(1x,7e11.5)') ( tlevel(l), l = 1,nlevel )
      read( iounit,'(1x,7e11.5)') ( player(l), l = 1,nlayer )
      read( iounit,'(1x,7e11.5)') ( tlayer(l), l = 1,nlayer )
      read( iounit,'(1x,7e11.5)') ( tlayer(l), l = 1,nlayer )
      read( iounit,'(1x,7e11.5)') ( tlayer(l), l = 1,nlayer )

* find two levels adjoining tropopause

      if(atmlbl.eq.'trp') ptrop= 93.7
      if(atmlbl.eq.'mls') ptrop=179.0
      if(atmlbl.eq.'mlw') ptrop=256.8
      if(atmlbl.eq.'sas') ptrop=267.7
      if(atmlbl.eq.'saw') ptrop=282.9

      do j=1,nlevel
         if(plevel(j).lt.ptrop.and.ptrop.le.plevel(j+1)) then
            jtrop1=j
            jtrop2=j+1
         endif
      end do

* interpolate fluxes to tropopause

      do k=0,0
         tropnt(k)=( (plevel(jtrop2)-ptrop) * fluxnt(jtrop1,k)
     +             + (ptrop-plevel(jtrop1)) * fluxnt(jtrop2,k))
     +             / (plevel(jtrop2)-plevel(jtrop1))
         tropup(k)=( (plevel(jtrop2)-ptrop) * fluxup(jtrop1,k)
     +             + (ptrop-plevel(jtrop1)) * fluxup(jtrop2,k))
     +             / (plevel(jtrop2)-plevel(jtrop1))
         tropdn(k)=( (plevel(jtrop2)-ptrop) * fluxdn(jtrop1,k)
     +             + (ptrop-plevel(jtrop1)) * fluxdn(jtrop2,k))
     +             / (plevel(jtrop2)-plevel(jtrop1))
      end do

* ingest flux and heating rate data band-by-band

      do k=1,nbands
         read( iounit,*) bandv1(k),bandv2(k)
         read( iounit,*) (fluxnt(l,k), l = 1,nlevel )
         read( iounit,*) (fluxup(l,k), l = 1,nlevel )
         read( iounit,*) (fluxdn(l,k), l = 1,nlevel )
         read( iounit,*) (heatrt(l,k), l = 1,nlayer )
         read( iounit,*) tropnt(k),tropup(k),tropdn(k)

         tropnt(k)=( (plevel(jtrop2)-ptrop) * fluxnt(jtrop1,k)
     +             + (ptrop-plevel(jtrop1)) * fluxnt(jtrop2,k))
     +             / (plevel(jtrop2)-plevel(jtrop1))
         tropup(k)=( (plevel(jtrop2)-ptrop) * fluxup(jtrop1,k)
     +             + (ptrop-plevel(jtrop1)) * fluxup(jtrop2,k))
     +             / (plevel(jtrop2)-plevel(jtrop1))
         tropdn(k)=( (plevel(jtrop2)-ptrop) * fluxdn(jtrop1,k)
     +             + (ptrop-plevel(jtrop1)) * fluxdn(jtrop2,k))
     +             / (plevel(jtrop2)-plevel(jtrop1))

      end do

* re-define units as ppm, molecules/cm**2

      do j=1,nlayer
         alayer(j) = 1.02 * (6.0256e+23/28.97) * (plevel(j+1)-plevel(j))
         wlayer(j) =  1.e+6 * wlayer(j) * 28.97/18.00
         olayer(j) =  1.e+6 * olayer(j) * 28.97/48.00
      end do

* sum column amounts for dry air and 5 gases

      pthair = 0.
      pthh2o = 0.
      ptho3  = 0.

      do j=1,nlayer
         pthair    = pthair + alayer(j)
         pthh2o    = pthh2o + 1.e-6 * wlayer(j) * alayer(j)
         ptho3     = ptho3  + 1.e-6 * olayer(j) * alayer(j)
      end do

      pthco2 = 1.e-6 * ppmco2 * pthair
      pthch4 = 1.e-6 * ppmch4 * pthair
      pthn2o = 1.e-6 * ppmn2o * pthair

* miscellaneous

      psurf = plevel(nlevel)
      tsurf = tlevel(nlevel)

      return

      end

      subroutine rgsfc
     +        (iounit,mlayer,mlevel,mbands
     +        ,header,atmlbl,iphase,icase,nlayer,nlevel,nbands
     +        ,psurf,tsurf,ptrop,flgctm,ppmco2,ppmch4,ppmn2o
     +        ,pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
     +        ,bandv1,bandv2,plevel,tlevel
     +        ,player,tlayer,alayer,wlayer,olayer
     +        ,fluxup,fluxdn,fluxnt,heatrt
     +        ,tropup,tropdn,tropnt)

*_______________________________________________________________________
* purpose: rgsfc - read gsfc icrccm data in earlier 1988 data format
*_______________________________________________________________________
* parameters defined in driver and passed to all subroutines:
      integer iounit,mlayer,mlevel,mbands

* variables passed to or created by subroutines:
      character*80 header(10)
      character*3  atmlbl
      integer      iphase,icase,nlayer,nlevel,nbands
      real         psurf,tsurf,ptrop,flgctm
      real         ppmco2,ppmch4,ppmn2o
      real         pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
      real         bandv1(0:mbands),bandv2(0:mbands)
      real         plevel(mlevel),tlevel(mlevel)
      real         player(mlayer),tlayer(mlayer)
      real         alayer(mlayer),wlayer(mlayer),olayer(mlayer)
      real         fluxup(mlevel,0:mbands)
      real         fluxdn(mlevel,0:mbands)
      real         fluxnt(mlevel,0:mbands)
      real         heatrt(mlayer,0:mbands)
      real         tropup(       0:mbands)
      real         tropdn(       0:mbands)
      real         tropnt(       0:mbands)
*_______________________________________________________________________

      character*80 csave

* blank out header initially

      do k=1,10
         write(header(k),'(80a1)') (' ',j=1,80)
      end do

* read 10 lines of header info from gsfc file

      read (iounit,'(16x,a80)') header

      do k=1,10
         csave=header(k)
         csave(1:77)=csave(4:80)
         header(k)=csave
      end do

* read key parameters which identify this case

      iphase=1
      read (iounit,910) icase,atmlbl,nlayer,nlevel,nbands,
     +                  psurf,tsurf,ptrop,
     +                  pthair,pthh2o,ptho3,flgctm,ppmco2,ppmch4,ppmn2o
 910  format(16x,i16,
     +     /,16x,13x,a3,
     +   3(/,16x,i16),
     +  10(/,16x,e16.6))

* ingest spectrally summed results (k=0) and all spectral band data

      do k=0,nbands
         read (iounit,'(16x,f16.2)') bandv1(k),bandv2(k)
         read (iounit,930) fluxup(nlevel,k),fluxdn(nlevel,k),
     +                     tropup(       k),tropdn(       k),
     +                     fluxup(     1,k)
 930  format(16x,e16.6,4(/,16x,e16.6))

         tropnt(k)=tropup(k)-tropdn(k)

         read (iounit,'(1x)')
         read (iounit,'(1x)')

         do j=1,nlayer
            read (iounit,'(1p,e16.6,3p,e16.4,1p,3e16.6)')
     +        player(j),tlayer(j),wlayer(j),olayer(j),heatrt(j,k)
         end do

         read (iounit,'(1x)')
         read (iounit,'(1x)')

         do j=1,nlevel
            read (iounit,'(1p,e16.6,3p,e16.4,1p,3e16.6)')
     +        plevel(j),tlevel(j),fluxup(j,k),fluxdn(j,k),fluxnt(j,k)
         end do

      end do

* re-define units as ppm, molecules/cm**2

      do j=1,nlayer
         alayer(j) = 1.02 * (6.0256e+23/28.97) * (plevel(j+1)-plevel(j))
         wlayer(j) = 1.e+6 * wlayer(j) * 28.97/18.00
         olayer(j) = 1.e+6 * olayer(j) * 28.97/48.00
      end do

* sum column amounts for dry air and 5 gases

      pthair = 0.
      pthh2o = 0.
      ptho3  = 0.

      do j=1,nlayer
         pthair    = pthair + alayer(j)
         pthh2o    = pthh2o + 1.e-6 * wlayer(j) * alayer(j)
         ptho3     = ptho3  + 1.e-6 * olayer(j) * alayer(j)
      end do

      pthco2 = 1.e-6 * ppmco2 * pthair
      pthch4 = 1.e-6 * ppmch4 * pthair
      pthn2o = 1.e-6 * ppmn2o * pthair

* icrccm phase, case # over-ridden using separate file

      open(1,file='temp_header')
      read(1,*) iphase
      read(1,*) icase
      read(1,'(a3)') atmlbl
      close(1)

      return
      end

      subroutine rfm93
     +        (iounit,mlayer,mlevel,mbands
     +        ,header,atmlbl,iphase,icase,nlayer,nlevel,nbands
     +        ,psurf,tsurf,ptrop,flgctm,ppmco2,ppmch4,ppmn2o
     +        ,pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
     +        ,bandv1,bandv2,plevel,tlevel
     +        ,player,tlayer,alayer,wlayer,olayer
     +        ,fluxup,fluxdn,fluxnt,heatrt
     +        ,tropup,tropdn,tropnt)

*_______________________________________________________________________
* purpose: rfm93 - read 1993-format icrccm data file, return data arrays
*_______________________________________________________________________
* parameters defined in driver and passed to all subroutines:
      integer iounit,mlayer,mlevel,mbands

* variables passed to or created by subroutines:
      character*80 header(10)
      character*3  atmlbl
      integer      iphase,icase,nlayer,nlevel,nbands
      real         psurf,tsurf,ptrop,flgctm
      real         ppmco2,ppmch4,ppmn2o
      real         pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
      real         bandv1(0:mbands),bandv2(0:mbands)
      real         plevel(mlevel),tlevel(mlevel)
      real         player(mlayer),tlayer(mlayer)
      real         alayer(mlayer),wlayer(mlayer),olayer(mlayer)
      real         fluxup(mlevel,0:mbands)
      real         fluxdn(mlevel,0:mbands)
      real         fluxnt(mlevel,0:mbands)
      real         heatrt(mlayer,0:mbands)
      real         tropup(       0:mbands)
      real         tropdn(       0:mbands)
      real         tropnt(       0:mbands)
*_______________________________________________________________________

* blank out header initially

      do k=1,10
         write(header(k),'(80a1)') (' ',j=1,80)
      end do

* read 10 lines of header info from fm93 file

      read (iounit,'(a80)') header

* read key parameters which identify this case

      read (iounit,910) iphase,icase,atmlbl,nlayer,nlevel,nbands,
     +                  psurf,tsurf,ptrop,flgctm,
     +                  ppmco2,ppmch4,ppmn2o,
     +                  pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
 910  format(24x,i16,
     +     /,24x,i16,
     +     /,24x,13x,a3,
     +   3(/,24x,i16),
     +   3(/,24x,f16.2),
     +   1(/,24x,f16.0),
     +   3(/,24x,f16.2),
     +   6(/,24x,f16.4))

* ingest spectrally summed results (k=0 only)

      k=0

      read (iounit,'(24x,f16.2)') bandv1(k),bandv2(k)
      read (iounit,'(24x,e16.4)') fluxup(nlevel,k),fluxdn(nlevel,k),
     +                           tropup(       k),tropdn(       k),
     +                           fluxup(     1,k)

      tropnt(k)=tropup(k)-tropdn(k)

      read (iounit,'(1x)')
      read (iounit,'(1x)')

      do j=1,nlayer
        read (iounit,'(6e12.4)')
     +    player(j),tlayer(j),alayer(j),wlayer(j),olayer(j),heatrt(j,k)
      end do

      read (iounit,'(1x)')
      read (iounit,'(1x)')

      do j=1,nlevel
         read (iounit,'(5e12.4)')
     +     plevel(j),tlevel(j),fluxup(j,k),fluxdn(j,k),fluxnt(j,k)
      end do

* read band-by-band flux summaries, fluxes, heating rates

      do k=1,nbands
         read(iounit,'(2f15.2)') bandv1(k),bandv2(k)
         read(iounit,'(5e15.6)') fluxup(nlevel,k),fluxdn(nlevel,k),
     +                           tropup(       k),tropdn(       k),
     +                           fluxup(     1,k)
         read(iounit,'(5e15.6)') (fluxup(j,k),j=1,nlevel)
         read(iounit,'(5e15.6)') (fluxdn(j,k),j=1,nlevel)
         read(iounit,'(5e15.6)') (fluxnt(j,k),j=1,nlevel)
         read(iounit,'(5e15.6)') (heatrt(j,k),j=1,nlayer)
      end do

      return
      end

      subroutine wfm93
     +        (iounit,mlayer,mlevel,mbands
     +        ,header,atmlbl,iphase,icase,nlayer,nlevel,nbands
     +        ,psurf,tsurf,ptrop,flgctm,ppmco2,ppmch4,ppmn2o
     +        ,pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
     +        ,bandv1,bandv2,plevel,tlevel
     +        ,player,tlayer,alayer,wlayer,olayer
     +        ,fluxup,fluxdn,fluxnt,heatrt
     +        ,tropup,tropdn,tropnt)

*_______________________________________________________________________
* purpose: wfm93 - write icrccm data in new standard (1993) format
*_______________________________________________________________________
* parameters defined in driver and passed to all subroutines:
      integer iounit,mlayer,mlevel,mbands

* variables passed to or created by subroutines:
      character*80 header(10)
      character*3  atmlbl
      integer      iphase,icase,nlayer,nlevel,nbands
      real         psurf,tsurf,ptrop,flgctm
      real         ppmco2,ppmch4,ppmn2o
      real         pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
      real         bandv1(0:mbands),bandv2(0:mbands)
      real         plevel(mlevel),tlevel(mlevel)
      real         player(mlayer),tlayer(mlayer)
      real         alayer(mlayer),wlayer(mlayer),olayer(mlayer)
      real         fluxup(mlevel,0:mbands)
      real         fluxdn(mlevel,0:mbands)
      real         fluxnt(mlevel,0:mbands)
      real         heatrt(mlayer,0:mbands)
      real         tropup(       0:mbands)
      real         tropdn(       0:mbands)
      real         tropnt(       0:mbands)
*_______________________________________________________________________

* print file header with 10 lines of info in any format

      write(iounit,'(a80)') header

* print descriptive parameters

      write(iounit,910) iphase,icase,atmlbl,nlayer,nlevel,nbands,
     +                  psurf,tsurf,ptrop,flgctm,
     +                  ppmco2,ppmch4,ppmn2o,
     +                  pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o
 910  format('.....icrccm.phase.number',i16
     +    ,/,'......icrccm.case.number',i16
     +    ,/,'.......atmospheric.label',a16
     +    ,/,'.....number.model.layers',i16
     +    ,/,'.....number.model.levels',i16
     +    ,/,'...number.spectral.bands',i16
     +    ,/,'........surface.pressure',0p,f16.2,'   mb'
     +    ,/,'.....surface.temperature',0p,f16.2,'   k'
     +    ,/,'.....tropopause.pressure',0p,f16.2,'   mb'
     +    ,/,'......h2o.continuum.flag',0p,f16.0,'   (0=off,1=included)'
     +    ,/,'........co2.mixing.ratio',0p,f16.2,'   ppm'
     +    ,/,'........ch4.mixing.ratio',0p,f16.2,'   ppm'
     +    ,/,'........n2o.mixing.ratio',0p,f16.2,'   ppm'
     +    ,/,'....total.dry.air.column',1p,e16.4,'   molecules/cm**2'
     +    ,/,'........total.h2o.column',1p,e16.4,'   molecules/cm**2'
     +    ,/,'........total.co2.column',1p,e16.4,'   molecules/cm**2'
     +    ,/,'.........total.o3.column',1p,e16.4,'   molecules/cm**2'
     +    ,/,'........total.ch4.column',1p,e16.4,'   molecules/cm**2'
     +    ,/,'........total.n2o.column',1p,e16.4,'   molecules/cm**2')

* print spectrally itegrated summary of icrccm fluxes

      k=0
      write(iounit,920) bandv1(k),bandv2(k)
 920  format('....spectral.region.from',0p,f16.2,'   cm-1  '
     +    ,/,'......................to',0p,f16.2,'   cm-1  ')
      write(iounit,930) fluxup(nlevel,k),fluxdn(nlevel,k),
     +                  tropup(       k),tropdn(       k),
     +                  fluxup(     1,k)
 930  format('.......flux.up.@.surface',1p,e16.4,'   watts/m**2'
     +    ,/,'.....flux.down.@.surface',1p,e16.4,'   watts/m**2'
     +    ,/,'....flux.up.@.tropopause',1p,e16.4,'   watts/m**2'
     +    ,/,'..flux.down.@.tropopause',1p,e16.4,'   watts/m**2'
     +    ,/,'..flux.up.@.top.of.atmos',1p,e16.4,'   watts/m**2')

* print table of atmospheric profile data and layer heating rates

      write(iounit,'(2a)')
     + ' layer >>  p           t    air mass'
     +,'         h2o          o3     heating'
     +,'        (mb)         (K)   (#/cm**2)'
     +,'       (ppm)       (ppm)   (deg/day)'

      do j=1,nlayer
         write(iounit,'(1p,e12.4,3p,e12.4,1p,4e12.4)')
     +     player(j),tlayer(j),alayer(j),wlayer(j),olayer(j),heatrt(j,k)
      end do

* print table of computed fluxes

      write(iounit,'(a)')
     + ' level >>  p           t     flux up   flux down    net flux'
     +,'        (mb)         (k)    (w/m**2)    (w/m**2)    (w/m**2)'

      do j=1,nlevel
         write(iounit,'(1p,e12.4,3p,e12.4,1p,3e12.4)')
     +     plevel(j),tlevel(j),fluxup(j,k),fluxdn(j,k),fluxnt(j,k)
      end do

* print band-by-band flux summaries, fluxes, heating rates

      do k=1,nbands
         write(iounit,'(0p,2f15.2)') bandv1(k),bandv2(k)
         write(iounit,'(1p,5e15.6)') fluxup(nlevel,k),fluxdn(nlevel,k),
     +                               tropup(       k),tropdn(       k),
     +                               fluxup(     1,k)
         write(iounit,'(1p,5e15.6)') (fluxup(j,k),j=1,nlevel)
         write(iounit,'(1p,5e15.6)') (fluxdn(j,k),j=1,nlevel)
         write(iounit,'(1p,5e15.6)') (fluxnt(j,k),j=1,nlevel)
         write(iounit,'(1p,5e15.6)') (heatrt(j,k),j=1,nlayer)
      end do

      return
      end