PROGRAM TABCMP_V ! ! VERSION (update VERSID) ! 04JUN20 AD v2.01 Use TABAUX for various auxiliary subroutines ! 27AUG19 AD v2.0 Converted from F77 to F90 ! 05AUG15 AD v1.21 Change '$' in WRITE statements to avoid gfortran warnings ! 24FEB15 AD v1.2: modified OPNTAB to improve detection of file-type ! 29DEC14 AD v1.1: modified RWAXIS to allow for format identifier ! Assume tabulation is ln(k) rather than k ! 08DEC14 AD Check for zero-absorption blocks ! 04NOV14 AD Original. ! ! DESCRIPTION ! Compress RFM .tab file in spectral domain ! NB if using binary files, use same compiler as the RFM. ! Compile with eg: gfortran tabaux.f90 tabcmp_v.f90 -o tabcmp_v ! ! USE TABAUX ! General modules for manipulating RFM .tab files ! IMPLICIT NONE ! ! GLOBAL CONSTANTS CHARACTER(7), PARAMETER :: VERSID = '2.01' ! Program version identifier INTEGER(I4), PARAMETER :: LUNOLD = 1 ! LUN for reading .tab file INTEGER(I4), PARAMETER :: LUNNEW = 2 ! LUN for new .tab file ! ! LOCAL VARIABLES LOGICAL :: BINNEW ! T=output to binary file, F=ASCII LOGICAL :: BINOLD ! T=input from binary file, F=ASCII INTEGER(I4) :: IDIV ! Counter for fine pts within convol.interval INTEGER(I4) :: IDUM(1) ! Dummy array for MINLOC result INTEGER(I4) :: IV ! Counter for orig.grd points INTEGER(I4) :: IV1, IV2 ! Range of orig.grd pts for current block INTEGER(I4) :: IVBEST ! Local index of grd pt with smallest error INTEGER(I4) :: IW ! Counter for convolved points INTEGER(I4) :: IW1, IW2 ! Limits for convolved spectral points INTEGER(I4) :: NDIV ! No. fine pts within each convol.pt interval INTEGER(I4) :: NV ! Original no. of spectral points INTEGER(I4) :: NVLEFT ! No. original pts left to read in file INTEGER(I4) :: NVNEW ! Counter for no.spectral points in new grid INTEGER(I4) :: NW ! Total no. of convolved points INTEGER(I4) :: NX ! Total number of p*T*q grid points REAL(R4) :: ERRLIM ! User-defined maximum absorption error REAL(R4) :: ERRMIN ! Minimum error for removing any point REAL(R8) :: DW ! Convolved resolution REAL(R8) :: V ! Individual spectral axis value REAL(R4) :: WGTTOT ! Total weight of all spectral pts CHARACTER(80) :: DWSTR ! User response to convolved resolution CHARACTER(80) :: RECORD ! Record read from file TYPE(TABTYP) :: TAB ! Tabulated data header INTEGER(I4), ALLOCATABLE :: IFLAG(:) ! Flag=1 or 2 for retained pts INTEGER(I4), ALLOCATABLE :: IVNEXT(:) ! [NV] Index of next kept pt INTEGER(I4), ALLOCATABLE :: IVPREV(:) ! [NV] Index of prev kept pt INTEGER(I4), ALLOCATABLE :: IWNEXT(:) ! [NV] Index of next cnv. pt INTEGER(I4), ALLOCATABLE :: IWPREV(:) ! [NV] Inex of prev cnv. pt INTEGER(I4), ALLOCATABLE :: IWV(:) ! [NV] Lower cnv.pt affctd by spc pt REAL(R4), ALLOCATABLE :: ERRV(:) ! [NV] Max.Error for each pt. REAL(R4), ALLOCATABLE :: ERRW(:,:) ! [NW,NX] Error on each cnv pt REAL(R4), ALLOCATABLE :: KTAB(:) ! [NX] ln(k) values for each pt REAL(R4), ALLOCATABLE :: LDU(:) ! [NX] Ln(Abs.Amt) in each cell REAL(R4), ALLOCATABLE :: LKU(:,:) ! [NV,NX] ln(k*U) values for spc.pt REAL(R8), ALLOCATABLE :: VTAB(:) ! [NV] List of wno values [cm-1] REAL(R4), ALLOCATABLE :: WGT(:) ! [NV] Weight of each spc.pt ! ! EXECUTABLE CODE ------------------------------------------------------------- ! WRITE ( *, '(A)' ) 'R-TABCMP_V: Running TABCMP_V v' // VERSID ! ! Open original .tab file CALL OPNTAB ( .TRUE., LUNOLD, BINOLD ) ! ! Read header data into TAB CALL RWAXIS ( .TRUE., LUNOLD, BINOLD, TAB ) ! ! Check original file is on a regular spectral grid IF ( ABS ( TAB%DV - ( TAB%V2 - TAB%V1 ) / DBLE ( TAB%NV - 1 ) ) .GE. & TAB%DV * 0.001 ) STOP 'F-TABCMP_V: Input file not on reg.grid' ! ! Calculate ln(abs.amt) in each 'cell' of tabulation NX = TAB%NX ALLOCATE ( LDU(NX) ) LDU = CLCLDU ( TAB ) ! ! Determine convolved resolution (if any) WRITE ( *, '(A)', ADVANCE='NO' ) 'Convolved resolution [cm-1]: ' READ ( *, '(A)' ) DWSTR ! IF ( DWSTR .EQ. '' ) THEN ! assume no convolution DW = TAB%DV ELSE READ ( DWSTR, * ) DW ! convolved resolution END IF ! ! Set weight for error at each point (assume triangular filter function) NDIV = NINT ( DW / TAB%DV ) ! No. spc. points to each convolved pt. WGTTOT = 1.0 / FLOAT(NDIV) ! Tot weight of fine grid to all o/p grid pts IW = 0 IDIV = NDIV NV = TAB%NV ALLOCATE ( IWV(NV), WGT(NV) ) DO IV = 1, NV IF ( IDIV .EQ. NDIV ) THEN IDIV = 0 IW = IW + 1 END IF IWV(IV) = IW ! Lower cnv.pt using IV WGT(IV) = ( NDIV - IDIV ) * WGTTOT**2 ! Wgt of contrib of IV to IW IDIV = IDIV + 1 END DO NW = IW ! WRITE ( *, '(A)', ADVANCE='NO' ) 'Max absorption error (eg 0.001): ' READ ( *, * ) ERRLIM ! CALL OPNTAB ( .FALSE., LUNNEW, BINNEW ) ! ! Read in data and convert to cell optical thickness ALLOCATE ( VTAB(NV), LKU(NV,NX), KTAB(NX) ) WRITE ( *, '(A)' ) 'I-TABCMP_V: Reading in LUT ...' DO IV = 1, NV CALL RWTAB ( .TRUE., LUNOLD, BINOLD, VTAB(IV), KTAB ) LKU(IV,:) = KTAB + LDU END DO ! IF ( MAXVAL ( LKU ) .LT. LOG ( ERRLIM ) ) & STOP 'W-TABCMP_V: Entire table has negligible absorption.' ! ALLOCATE ( IFLAG(NV), IVPREV(NV), IVNEXT(NV), IWPREV(NV), IWNEXT(NV) ) IFLAG = 1 IFLAG(1) = 2 IFLAG(NV) = 2 IVPREV = (/ 1, (IV-1, IV=2,NV) /) ! 1,1,2,3 ... NV-1 IVNEXT = (/ (IV+1,IV=1,NV-1), NV /) ! 2,3, ... NV, NV IWPREV = IWV IWNEXT = IWV + 1 ! ! Set up array containing current convolved pt errors ALLOCATE ( ERRW(NW,NX) ) ERRW = 0.0 ! ! Evaluate initial interpolation errors for each spectral point ALLOCATE ( ERRV(NV) ) CALL INIERR ( NX, NV, WGTTOT, WGT, LKU, ERRV ) ! WRITE ( *, '(A,I8)' ) 'I-TABCMP_V: Starting to remove ' // & 'points. Initial No.=', NV ! NVLEFT = NV ! Repeat from here after each point is removed DO ! Find point with the smallest error ERRMIN = MINVAL ( ERRV, MASK=IFLAG .EQ. 1 ) IF ( ERRMIN .GT. ERRLIM ) THEN WRITE ( *, '(4X, I10, A10, ES13.4)' ) NVLEFT, ' Err.lim.=', ERRMIN EXIT ! no more points END IF IDUM = MINLOC ( ERRV, MASK=IFLAG .EQ. 1 ) IVBEST = IDUM(1) ! ! Remove point with the smallest error and write progress every 100 pts NVLEFT = NVLEFT - 1 IF ( MOD ( NVLEFT, 10000 ) .EQ. 0 .OR. NVLEFT .EQ. NV-1 ) & WRITE ( *, '(4X,I10, F10.4, ES13.4, F6.1, A)' ) & NVLEFT, VTAB(IVBEST), ERRMIN, FLOAT(NVLEFT)/FLOAT(NV)*100.0, '%' ! ! Remove point IFLAG(IVBEST) = 0 IV1 = IVPREV(IVBEST) IVNEXT(IV1) = IVNEXT(IVBEST) IWNEXT(IV1) = IWNEXT(IVBEST) IV2 = IVNEXT(IVBEST) IVPREV(IV2) = IVPREV(IVBEST) IWPREV(IV2) = IWPREV(IVBEST) ! ! Update ERRW with current interpolation errors on each convolved point ! And replace LKU with interpolated values CALL UPDATE ( NX, WGTTOT, IWV, WGT, LKU, ERRW, IV1, IV2 ) ! ! Fine grid range IV1:IV2 is now interpolated, so find range of convolved grid ! points for which error has changed IW1 = IWPREV(IVBEST) IW2 = IWNEXT(IVBEST) ! ! Now find range JV1:JV2 of fine pts whose removal influences changed cnv pts IV = IV1 DO IF ( IWNEXT(IV) .LT. IW1 ) EXIT IV = IVPREV(IV) if ( iflag(iv) .eq. 0 ) stop 'logical error#0' IF ( IFLAG(IV) .EQ. 2 ) EXIT END DO ! ! Evalute new interpolation error for each spc.pt within interval DO IV1 = IVPREV(IV) IV2 = IVNEXT(IV) ERRV(IV) = ERRINT ( NX, WGTTOT, IWV, WGT, LKU, IV1, IV2, ERRW ) IV = IV2 IF ( IFLAG(IV) .EQ. 2 ) EXIT ! Reached upper end of spectral range IF ( IWPREV(IV) .GT. IW2 ) EXIT ! Reached upper end of affected range END DO END DO ! ! Continue from here when no more points can be removed NVNEW = COUNT ( IFLAG .NE. 0 ) WRITE ( *, '(A,I9,A,I9,A,F5.1,A)' ) 'I-TABCMP_V: Summary: Orig Npts=', & NV, ' New Npts=', NVNEW, ' Red.Factor=', & 100 * FLOAT ( NVNEW ) / FLOAT ( NV ), '%' WRITE ( *, '(A)' ) 'I-TABCMP_V: Writing new .tab file ...' ! RECORD = '! Reduced grid using program TABCMP_V v.' // VERSID CALL COPHDR ( LUNOLD, LUNNEW, BINOLD, BINNEW, RECORD ) ! also rewinds LUNOLD ! ! Re-read header data into TAB CALL RWAXIS ( .TRUE., LUNOLD, BINOLD, TAB ) TAB%NV = NVNEW CALL RWAXIS ( .FALSE., LUNNEW, BINNEW, TAB ) ! DO IV = 1, NV CALL RWTAB ( .TRUE., LUNOLD, BINOLD, V, KTAB ) IF ( IFLAG(IV) .NE. 0 ) & CALL RWTAB ( .FALSE., LUNNEW, BINNEW, V, KTAB ) END DO CLOSE ( LUNOLD ) CLOSE ( LUNNEW ) ! WRITE (*,'(A)') 'R-TABCMP_V: Successful completion' ! END PROGRAM TABCMP_V