! F90 Program GEIXSC ! ! VERSION (update VERSID) ! 16NOV18 AD v2.00 F90 version ! 05AUG15 AD v1.01 Change '$' in WRITE statements to avoid gfortran warnings ! 25SEP13 AD v1.00 Original. ! ! DESCRIPTION ! Convert GEISA cross-section data to HITRAN .xsc format. ! See http://eodg.atm.ox.ac.uk/RFM/geixsc.html for documentation ! User supplies a file containing list of GEISA filenames for conversion. ! GEISA files may contain more than one spectral band and non-uniform ! spectral grids - this program generates separate files for each spectral ! band and puts these on to a regular grid required for HITRAN format. ! ! Fatal errors result in a message 'F-GEIXSC ...' printed to terminal. ! Information messages begin 'I-GEIXSC ...' ! ! The structure of this file has the GEIXSC MODULE before the main program in ! order to simplify the compilation (eg use: gfortran geixsc.f90 -o geixsc ) ! MODULE GEIXSC ! VARIABLE KINDS ! GLOBAL DATA ! GLOBAL TYPES ! CONTAINS ! SUBROUTINE REAGEI ! SUBROUTINE CHKGRD ! SUBROUTINE SPCGEI ! SUBROUTINE REGRID ! ! REFERENCES ! GEISA format: http://cds-espri.ipsl.upmc.fr/etherTypo/?id=1778 ! HITRAN format: Table 1 of http://hitran.org/media/refs/HITRAN-2004.pdf ! MODULE GEIXSC_MOD ! ! VARIABLE KINDS INTEGER, PARAMETER :: I4 = SELECTED_INT_KIND(9) INTEGER, PARAMETER :: R4 = SELECTED_REAL_KIND(6) INTEGER, PARAMETER :: R8 = SELECTED_REAL_KIND(15,200) ! ! GLOBAL DATA CHARACTER(4) :: VERSID = '2.00' ! ! GLOBAL TYPES TYPE :: GEITYP ! Structure for GEISA .xsc header info INTEGER(I4) :: NPT ! Total No. data points REAL(R4) :: PRE ! Pressure [Pa] REAL(R4) :: RES ! Resolution [cm-1] REAL(R4) :: TEM ! Temperature [K] REAL(R8) :: WN1 ! Lower Wno [cm-1] of GEISA file REAL(R8) :: WN2 ! Upper Wno [cm-1] of GEISA file CHARACTER(20) :: MOL ! Molecule name from GEISA file CHARACTER(15) :: REF ! Reference REAL(R8), POINTER :: WNO(:) ! GEISA wavenumber values [cm-1] REAL(R4), POINTER :: DAT(:) ! GEISA absorption cross-sections [cm2/molec] END TYPE GEITYP ! CONTAINS SUBROUTINE REAGEI ( LUN, FIL, GEI ) ! ! VERSION ! 09NOV18 AD F90 conversion ! ! DESCRIPTION ! Read GEISA x/s file ! IMPLICIT NONE ! ! ARGUMENTS INTEGER(I4), INTENT(IN) :: LUN ! LUN for reading data CHARACTER(*), INTENT(IN) :: FIL ! Name of GEISA .xsc file TYPE(GEITYP), INTENT(OUT) :: GEI ! GEISA data ! ! LOCAL VARIABLES INTEGER(I4) :: IHDR ! Counter for header records INTEGER(I4) :: IOS ! Saved value of IOSTAT INTEGER(I4) :: IPT ! Counter for data points/records INTEGER(I4) :: NPT ! Number data points in file CHARACTER(80) :: REC ! Record read from file ! ! EXECUTABLE CODE ------------------------------------------------------------- ! OPEN ( UNIT=LUN, FILE=FIL, STATUS='OLD', ACTION='READ', ERR=900, IOSTAT=IOS ) ! ! First 9 records are header records REC = '' DO IHDR = 1, 9 READ ( LUN, '(A)', ERR=900, IOSTAT=IOS ) REC IF ( REC(1:1) .NE. '#' ) THEN WRITE ( *, * ) 'F-REAGEI: unexpected header format in file: ' // FIL STOP END IF IF ( IHDR .EQ. 6 ) THEN ! Rec#6 contains header data READ ( REC, '(X,A20,X,F6.2,X,F9.2,F8.4,A15,I9,2F11.4)', & ERR=900, IOSTAT=IOS ) & GEI%MOL, GEI%TEM, GEI%PRE, GEI%RES, GEI%REF, GEI%NPT, GEI%WN1, GEI%WN2 IF ( IOS .NE. 0 ) THEN WRITE ( *, * ) 'F-REAGEI: Unexpected header contents in file: ' // FIL WRITE ( *, * ) 'F-REAGEI: Error reading header data, IOSTAT=', IOS STOP END IF END IF END DO ! NPT = GEI%NPT ALLOCATE ( GEI%WNO(NPT), GEI%DAT(NPT) ) READ ( LUN, *, ERR=900, IOSTAT=IOS ) & ( GEI%WNO(IPT), GEI%DAT(IPT), IPT = 1, NPT ) CLOSE ( LUN, IOSTAT=IOS, ERR=900 ) ! RETURN ! Normal exit ! 900 CONTINUE WRITE ( *, * ) 'IOSTAT=', IOS, ' error reading file: ', FIL STOP 'F-REAGEI: fatal error' ! END SUBROUTINE REAGEI SUBROUTINE CHKGRD ( WNOGRD, DATGRD ) ! ! VERSION ! 09NOV18 AD Original. ! ! DESCRIPTION ! Check GEISA spectral grid is monotonic and remove any duplicate points ! IMPLICIT NONE ! ! ARGUMENTS REAL(R8), INTENT(INOUT) :: WNOGRD(:) ! Original wavenumbe grid REAL(R4), INTENT(INOUT) :: DATGRD(:) ! Data on original grid ! ! LOCAL VARIABLES INTEGER(I4) :: IGRD ! Counter for input grid points INTEGER(I4) :: JGRD ! Secondary counter for input grid points INTEGER(I4) :: NGRD ! No. points in grid ! ! EXECUTABLE CODE -------------------------------------------------------------- ! NGRD = SIZE ( WNOGRD ) JGRD = 1 DO IGRD = 2, NGRD ! Check for duplicate wavenumber grid points IF ( WNOGRD(IGRD) .GT. WNOGRD(IGRD-1) ) THEN JGRD = JGRD + 1 ELSE IF ( WNOGRD(IGRD) .LT. WNOGRD(IGRD-1) ) THEN WRITE ( *, '(A,2F10.4,A)' ) 'F-REGRID: Wavenumber grid for range ', & WNOGRD(1), WNOGRD(NGRD), ' [cm-1] does not increase monotonically' STOP 'F-REGRID: Fatal error' END IF WNOGRD(JGRD) = WNOGRD(IGRD) DATGRD(JGRD) = DATGRD(IGRD) END DO ! ! Warn if any duplicate spectral grid points have been removed IF ( JGRD .NE. NGRD ) THEN WRITE ( *, '(A,I6,A,I6)' ) 'W-GEIXSC: Removing ', NGRD-JGRD, & ' duplicate spectral points from original ', NGRD NGRD = JGRD STOP 'haven''t yet handled changing original grid' END IF ! END SUBROUTINE CHKGRD SUBROUTINE SPCGEI ( GEI, NSPC, SPC ) ! ! VERSION ! 09NOV18 AD F90 conversion ! ! DESCRIPTION ! Identify spectral bands within GEISA data ! IMPLICIT NONE ! ! ARGUMENTS TYPE(GEITYP), INTENT(IN) :: GEI ! GEISA data INTEGER(I4), INTENT(OUT) :: NSPC ! No.different spectral bands INTEGER(I4), ALLOCATABLE, INTENT(OUT) :: SPC(:,:) ! Start/End pts in each spectral band ! ! LOCAL CONSTANTS INTEGER(I4), PARAMETER :: MINPTS = 10 ! Min no. points in spectral subset REAL(R4), PARAMETER :: WNOGAP = 5.0 ! Min separation between spc.subsets ! ! LOCAL VARIABLES INTEGER(I4) :: IPT ! Counter for spectral points/records INTEGER(I4) :: IPT1 ! First point in spectral band INTEGER(I4) :: IPT2 ! Last point in spectral band INTEGER(I4), ALLOCATABLE :: SPCSAV(:,:) ! Saved value of SPC ! ! EXECUTABLE CODE ------------------------------------------------------------- ! ! Read through spectral data defining new spectral range by gaps in waveno. NSPC = 0 IPT2 = 1 IPT1 = 2 DO DO IPT = IPT1+1, GEI%NPT IF ( ( GEI%WNO(IPT) - GEI%WNO(IPT2) ) .GT. WNOGAP ) EXIT IPT2 = IPT END DO IF ( IPT2 - IPT1 + 1 .GE. MINPTS ) THEN ! new spectral band IF ( NSPC .GT. 0 ) CALL MOVE_ALLOC ( SPC, SPCSAV ) NSPC = NSPC + 1 ALLOCATE ( SPC(2,NSPC) ) IF ( NSPC .GT. 1 ) SPC(:,1:NSPC-1) = SPCSAV SPC(:,NSPC) = (/ IPT1, IPT2 /) END IF IF ( IPT2 .EQ. GEI%NPT ) EXIT IPT2 = IPT2 + 1 IPT1 = IPT2 END DO ! END SUBROUTINE SPCGEI SUBROUTINE REGRID ( WNOGRD, DATGRD, NREG, WN1REG, DWNREG, DATREG ) ! ! VERSION ! 01NOV18 AD Original. ! ! DESCRIPTION ! Reinterpolate spectral data to regular grid. ! IMPLICIT NONE ! ! ARGUMENTS REAL(R8), INTENT(IN) :: WNOGRD(:) ! Original wavenumbe grid REAL(R4), INTENT(IN) :: DATGRD(:) ! Data on original grid INTEGER(I4), INTENT(OUT) :: NREG ! No. regular data points REAL(R8), INTENT(OUT) :: WN1REG ! 1st Wno.pt [cm-1] on regular grid REAL(R8), INTENT(OUT) :: DWNREG ! Wno increment [cm-1] on regular grid REAL(R4), ALLOCATABLE, & INTENT(OUT) :: DATREG(:) ! Data interpolated to regular grid ! ! LOCAL CONSTANTS REAL(R4), PARAMETER :: DPREC = 0.08 ! Precision of wno grid ! ! LOCAL VARIABLES INTEGER(I4) :: IGRD ! Counter for input grid points INTEGER(I4) :: IMAX ! Saved location of best grid point INTEGER(I4) :: IREG ! Counter for points on new regular grid INTEGER(I4) :: JGRD ! Secondary counter for input grid points INTEGER(I4) :: NFIT ! No. of data pts fitted without interpolation INTEGER(I4) :: NGRD ! Size of original grid INTEGER(I4) :: NMAX ! Max. value of NFIT found so far REAL(R4) :: D ! Fractional grid position REAL(R8) :: DGRD ! Difference between adjacent pts on original grid REAL(R8) :: DMAX ! Best estimate of grid spacing REAL(R8) :: DMIN ! Minimum spacing [cm-1] for new grid REAL(R8) :: WNO ! Wavenumber [cm-1] from new regular grid ! ! EXECUTABLE CODE -------------------------------------------------------------- ! NGRD = SIZE ( WNOGRD ) ! ! Start by looking at end points and see if all wavenumber grid points fit ! a linear interpolation WN1REG = WNOGRD(1) DWNREG = ( WNOGRD(NGRD) - WN1REG ) / ( NGRD - 1 ) NFIT = 0 WNO = WN1REG DO IGRD = 1, NGRD D = SNGL ( ( WNO - WNOGRD(IGRD) ) / DWNREG ) IF ( ABS ( D - NINT(D) ) .LT. DPREC ) NFIT = NFIT + 1 WNO = WNO + DWNREG END DO ! ! Normal exit for a regular grid IF ( NFIT .EQ. NGRD ) THEN NREG = NGRD ALLOCATE ( DATREG(NREG) ) DATREG = DATGRD RETURN END IF ! ! Next looking at grid spacing of adjacent data points to see if these can ! be used to construct a regular grid which contains most of the other data ! points without interpolation. ! ! It could be that two grid points are really closely spaced, so also restrict ! the constructed grid resolution to be not much smaller than average spacing DMIN = 0.9 * DWNREG ! NMAX = 0 IMAX = 1 DO IGRD = 1, NGRD-1 DGRD = WNOGRD(IGRD+1) - WNOGRD(IGRD) IF ( DGRD .GT. DMIN .AND. DGRD .LT. 2*DMIN ) THEN NFIT = 0 DO JGRD = 1, NGRD D = SNGL ( ( WNOGRD(JGRD) - WNOGRD(IGRD) ) / DGRD ) IF ( ABS ( D - NINT(D) ) .LT. DPREC ) THEN NFIT = NFIT + 1 IF ( NINT(D) .NE. 0 ) & DGRD = ( WNOGRD(JGRD) - WNOGRD(IGRD) ) / NINT(D) END IF END DO IF ( NFIT .GT. NMAX ) THEN NMAX = NFIT IMAX = IGRD DMAX = DGRD END IF ! If 90% of data points lie on this grid, don't look any further IF ( NMAX .GT. 0.9 * NGRD ) EXIT END IF END DO DWNREG = DMAX ! ! If 10% or more of data points lie on this new grid, then use it .. IF ( NMAX .GT. 0.1 * NGRD ) THEN ! find lowest wno point within data IGRD = NINT ( ( WNOGRD(1) - WNOGRD(IMAX) ) / DWNREG ) ! could be 0 WN1REG = WNOGRD(IMAX) + IGRD * DWNREG IF ( WN1REG .LT. WNOGRD(1) - 0.01 * DWNREG ) & WN1REG = WN1REG + DWNREG ! next pt NREG = INT ( ( WNOGRD(NGRD) - WN1REG ) / DWNREG ) + 1 ! round down ! ... otherwise just reinterpolate original grid to regular spacing IF ( NGRD .EQ. NREG ) THEN ! only minor adjustment to wno axis WRITE ( *, '(A,2F10.4,A,F10.6,A)' ) 'I-REGRID: Adjusting ', WNOGRD(1), & WNOGRD(NGRD), ' spectral range to uniform ', DWNREG, ' [cm-1] spacing' ELSE WRITE ( *, '(A,2F10.4,A,I6,A,I6,A,F10.6,A)' ) 'W-REGRID: Interpolating ',& WNOGRD(1), WNOGRD(NGRD), ' from', NGRD, ' to', NREG, ' points at ', & DWNREG, ' [cm-1] spacing' END IF ELSE NREG = NGRD WN1REG = WNOGRD(1) DWNREG = ( WNOGRD(NGRD) - WNOGRD(1) ) / ( NREG - 1 ) WRITE ( *, '(A,2F10.4,A,F10.6,A)' ) 'W-REGRID: Re-interpolating ', & WNOGRD(1), WNOGRD(NGRD), ' to uniform grid at ', DWNREG, ' [cm-1] spacing' END IF ! ! Interpolate original data to new regular grid ALLOCATE ( DATREG(NREG) ) WNO = WN1REG IGRD = 1 DO IREG = 1, NREG DO WHILE ( WNOGRD(IGRD+1) .LT. WNO .AND. IGRD .LT. NGRD-1 ) IGRD = IGRD + 1 END DO D = SNGL ( ( WNO - WNOGRD(IGRD) ) / ( WNOGRD(IGRD+1) - WNOGRD(IGRD) ) ) ! Don't alter original data points unless significant interpolation required IF ( ABS ( D ) .LT. DPREC ) D = 0.0 IF ( ABS ( D - 1.0 ) .LT. DPREC ) D = 1.0 DATREG(IREG) = D * DATGRD(IGRD+1) + (1.0-D)*DATGRD(IGRD) WNO = WNO + DWNREG END DO ! ! Normal exit with new grid ! END SUBROUTINE REGRID END MODULE GEIXSC_MOD PROGRAM GEIXSC ! ! DESCRIPTION ! Main program ! USE GEIXSC_MOD ! IMPLICIT NONE ! ! LOCAL CONSTANTS INTEGER(I4), PARAMETER :: LUNGEI = 1 ! LUN for GEISA input files INTEGER(I4), PARAMETER :: LUNLST = 2 ! LUN for list of GEISA filenames INTEGER(I4), PARAMETER :: LUNXSC = 3 ! LUN for output .xsc files REAL(R4), PARAMETER :: PA2TOR = 0.0075 ! Conversion of Pa to Torr ! approx 760.0/101325.0 ! LOCAL VARIABLES INTEGER(I4) :: IOS ! Saved value of IOSTAT INTEGER(I4) :: IPT1 ! Start index of spectral band within GEISA data INTEGER(I4) :: IPT2 ! End index of spectral band within GEISA data INTEGER(I4) :: ISPC ! Counter for spc.ranges within file INTEGER(I4) :: NGRD ! Number of abs.coeff. data points INTEGER(I4) :: NSPC ! No. spcranges within file for current molec. REAL(R8) :: DWGRD ! Resolution [cm-1] of output grid REAL(R8) :: W1GRD ! Lower Wno [cm-1] of HITRAN output grid REAL(R8) :: W2GRD ! Upper Wno [cm-1] of HITRAN output grid CHARACTER(200) :: FILGEI ! Name of next GEISA file to be read CHARACTER(200) :: FILLST ! File containing list of GEISA filenames CHARACTER(50) :: FILXSC ! Name of output [molec].xsc file CHARACTER(5) :: PSTR ! Pressure [Torr] inserted into HITRAN filename CHARACTER(5) :: TSTR ! Temperature [K] inserted into HITRAN filename CHARACTER(6) :: W1STR ! Lower Wno [cm-1] inserted into HITRAN filename CHARACTER(6) :: W2STR ! Upper Wno [cm-1] inserted into HITRAN filename TYPE(GEITYP) :: GEI ! GEISA .xsc file contents INTEGER(I4), ALLOCATABLE :: SPC(:,:) ! Indices of spectral bands REAL(R4), ALLOCATABLE :: ACSGRD(:) ! Abs.x/s data on output grid ! ! EXECUTABLE CODE ------------------------------------------------------------- ! WRITE ( *, '(A)' ) 'R-GEIXSC: Running program GEIXSC v' // VERSID ! WRITE ( *, '(A)', ADVANCE='NO' ) 'File containing list of GEISA filenames: ' READ ( *, '(A)' ) FILLST ! OPEN ( UNIT=LUNLST, FILE=FILLST, STATUS='OLD', ACTION='READ', IOSTAT=IOS ) IF ( IOS .NE. 0 ) STOP 'F-GEIXSC: Error opening file' ! ! Read through all GEISA files DO READ ( LUNLST, '(A)', IOSTAT=IOS ) FILGEI IF ( IOS .GT. 0 ) THEN WRITE ( *, * ) & 'F-GEIXSC: Error reading list of GEISA files, IOSTAT=', IOS STOP END IF IF ( IOS .LT. 0 ) EXIT ! assume end-of-file reached IF ( FILGEI .EQ. '' ) CYCLE ! Skip any blank records WRITE (*,*) 'I-GEIXSC: Reading file: ' // TRIM ( FILGEI ) CALL REAGEI ( LUNGEI, FILGEI, GEI ) ! Read GEISA file CALL CHKGRD ( GEI%WNO, GEI%DAT ) CALL SPCGEI ( GEI, NSPC, SPC ) ! Split into spectral ranges ! DO ISPC = 1, NSPC ! IPT1 = SPC(1,ISPC) IPT2 = SPC(2,ISPC) CALL REGRID ( GEI%WNO(IPT1:IPT2), GEI%DAT(IPT1:IPT2), & NGRD, W1GRD, DWGRD, ACSGRD ) W2GRD = W1GRD + ( NGRD - 1 ) * DWGRD ! ! Create output file using molecule name ! 1234567890123456789012345678901234567890 ! eg CCl3F_233.1_198.0_810.0-880.0_00.xsc WRITE ( TSTR, '(F5.1)' ) GEI%TEM WRITE ( PSTR, '(F5.1)' ) MAX ( 0.0, GEI%PRE*PA2TOR ) ! Convert Pa to Torr WRITE ( W1STR, '(F6.1)' ) W1GRD WRITE ( W2STR, '(F6.1)' ) W2GRD FILXSC = TRIM ( ADJUSTL ( GEI%MOL ) ) // '_' // TSTR // '_' // & TRIM ( ADJUSTL ( PSTR ) ) // '_' // & TRIM ( ADJUSTL ( W1STR ) ) // '-' // TRIM ( ADJUSTL ( W2STR ) ) & // '.xsc' WRITE ( *, * ) 'I-GEIXSC: Writing file: ' // FILXSC OPEN ( UNIT=LUNXSC, FILE=FILXSC, STATUS='UNKNOWN', ACTION='WRITE', & ERR=900, IOSTAT=IOS ) ! WRITE ( LUNXSC, '(A20,2F10.4,I7,F7.2,F6.1,ES10.3,F5.2,A15,7X,A3)',& IOSTAT=IOS, ERR=900 ) & GEI%MOL, W1GRD, W2GRD, NGRD, GEI%TEM, MAX(0.0,GEI%PRE*PA2TOR), & MAXVAL(ACSGRD), DWGRD, GEI%MOL(6:20), GEI%REF(13:15) WRITE ( LUNXSC, '(10ES10.3)' ) ACSGRD CLOSE ( LUNXSC, ERR=900, IOSTAT=IOS ) END DO END DO ! IFIL loop ! STOP 'R-GEIXSC: Successful completion' ! 900 CONTINUE WRITE ( *, * ) 'F-GEIXSC: I/O error on output file, IOSTAT=', IOS STOP ! END PROGRAM GEIXSC