! F90 program TABCMP_X ! ! VERSION (update VERSID) ! 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 ! 16DEC14 AD Original. ! 01NOV14 AD Original. ! ! DESCRIPTION ! Subsample RFM .tab file and create irreg.grd file ! MODULE TABCMP_X_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 CONSTANTS CHARACTER(7), PARAMETER :: VERSID = '2.0' ! Program version identifier INTEGER(I4), PARAMETER :: LUNOLD = 1 ! LUN for reading .tab file INTEGER(I4), PARAMETER :: LUNNEW = 2 ! LUN for new .tab file INTEGER(I4), PARAMETER :: ND = 3 ! No of tabulation axes (3=q,T,p) ! ! GLOBAL VARIABLES LOGICAL :: BINNEW ! T=binary new .tab file, F=ASCII LOGICAL :: BINOLD ! T=binary old .tab file, F=ASCII INTEGER(I4) :: NMAX ! Max of NQ,NT or NP INTEGER(I4) :: NP ! No. original pressure axis points INTEGER(I4) :: NQ ! No. original VMR axis points INTEGER(I4) :: NT ! No. original temperature axis points INTEGER(I4) :: NTQ ! No.T axis pts * No.q axis points ! ! GLOBAL TYPES TYPE :: GRDTYP INTEGER(I4) :: FLG ! Flag 0=remove, 1=keep, 2=fixed pt INTEGER(I4) :: IAN ! Index of next retainted axis point INTEGER(I4) :: IAP ! Index of curr/prev retained axis point REAL(R4) :: AXV ! Tabulation axis values END TYPE GRDTYP TYPE :: TABTYP REAL(R4) :: FMT ! Tab. file format identifier CHARACTER(5) :: MOL ! Molecule/isotope ID INTEGER(I4) :: NV ! No. spectral points in file REAL(R8) :: V1 ! First spectral point in file REAL(R8) :: V2 ! Last spectral point in file REAL(R8) :: DV ! Nominal spectral grid interval INTEGER(I4) :: NX ! Total number of p*T*q grid points INTEGER(I4) :: NP ! Number of p-axis points INTEGER(I4) :: NT ! Number of T-axis points INTEGER(I4) :: NQ ! Number of q-axis points REAL(R4), POINTER :: PAX(:) ! [NP] Pressure axis [mb] REAL(R4), POINTER :: TEM(:) ! [NP] Temperature profile [K] REAL(R4), POINTER :: VMR(:) ! [NP] VMR profile [ppmv] REAL(R4), POINTER :: QAX(:) ! [NQ] VMR scale factor axis [%] REAL(R4), POINTER :: TAX(:) ! [NT] Temperature axis [K] END TYPE TABTYP ! CONTAINS ! LOGICAL FUNCTION BINLUT ( LUN, FILNAM ) ! ! VERSION ! 21MAR19 AD Original. ! ! DESCRIPTION ! Detect whether .tab files is binary (T) or ASCII (F) ! IMPLICIT NONE ! ! ARGUMENTS INTEGER(I4), INTENT(IN) :: LUN ! Logical Unit Number CHARACTER(*), INTENT(IN) :: FILNAM ! Filename ! ! LOCAL VARIABLES LOGICAL :: BINFIL ! Current assumption of binary (T) or ASCII (F) INTEGER(I4) :: IOS ! Saved value of IOSTAT for error messages INTEGER(I4) :: ITRY ! Counter for attempt to open existing file REAL(R4) :: RFMT ! Format identifier CHARACTER(80) :: RECORD ! Text header record ! ! EXECUTABLE CODE -------------------------------------------------------------- ! ! Initially assume ASCII unless filename contains string '.bin' BINFIL = INDEX ( FILNAM, '.bin' ) .GT. 0 DO ITRY = 1, 2 RECORD(1:1) = '!' IF ( BINFIL ) THEN OPEN ( UNIT=LUN, FILE=FILNAM, STATUS='OLD', ACTION='READ', & FORM='UNFORMATTED', IOSTAT=IOS, ERR=100 ) DO WHILE ( RECORD(1:1) .EQ. '!' ) READ ( LUN, IOSTAT=IOS, ERR=100 ) RECORD END DO ELSE OPEN ( UNIT=LUN, FILE=FILNAM, STATUS='OLD', ACTION='READ', & FORM='FORMATTED', IOSTAT=IOS, ERR=100 ) DO WHILE ( RECORD(1:1) .EQ. '!' ) READ ( LUN, '(A)', IOSTAT=IOS, ERR=100 ) RECORD END DO END IF READ ( RECORD, *, IOSTAT=IOS ) RFMT 100 CONTINUE CLOSE ( LUN ) IF ( IOS .EQ. 0 ) THEN BINLUT = BINFIL RETURN END IF BINFIL = .NOT. BINFIL ! first attempt failed, so try again END DO ! STOP 'F-BINLUT: Unable to determine if binary or ASCII LUT file' ! END FUNCTION BINLUT FUNCTION CLCLDU ( TAB ) ! ! VERSION ! 21MAR19 AD F90 original. ! ! DESCRIPTION ! Calculate ln(absorber amount [kmol/m2]) in each 'cell' of .tab file ! IMPLICIT NONE ! ! ARGUMENTS TYPE(TABTYP), INTENT(IN) :: TAB ! TAB file header ! ! FUNCTION TYPE REAL(R4) :: CLCLDU(TAB%NX) ! Function returns array size of NX=NP*NT*NQ ! ! LOCAL CONSTANTS REAL(R4), PARAMETER :: GRAVTY = 9.81 ! Acceleration due to gravity [m/s^2] REAL(R4), PARAMETER :: WGTAIR = 28.96 ! Weight of air [kg/kmole] ! ! LOCAL VARIABLES LOGICAL :: PHI2LO ! T=pressure axis starts with max p, F=min p INTEGER(I4) :: IP ! Counter for pressure axis points INTEGER(I4) :: IQ ! Counter for VMR axis points INTEGER(I4) :: IT ! Counter for temperature axis points INTEGER(I4) :: IX ! Counter for p,T,q axis points REAL(R4) :: DP ! Difference between pressure axis values ! ! EXECUTABLE CODE ------------------------------------------------------------- ! ! write (*,*) 'clcldu: changing definition of DP' PHI2LO = TAB%PAX(1) .GT. TAB%PAX(TAB%NP) IX = 0 DO IP = 1, TAB%NP DP = TAB%PAX(IP) ! pressure difference between grid points IF ( PHI2LO ) THEN IF ( IP .LT. TAB%NP ) DP = DP - TAB%PAX(IP+1) ELSE IF ( IP .GT. 1 ) DP = DP - TAB%PAX(IP-1) END IF dp = tab%pax(ip) DO IT = 1, TAB%NT DO IQ = 1, TAB%NQ IX = IX + 1 ! QAX is a percentage, DP is in [mb] so factors 100 cancel for DP in [Pa] ! 1.0e-6 to convert VMR [ppmv] to [ppv], DU in [kmole/m^2] CLCLDU(IX) = LOG ( MAX ( TINY(1.0), & DP * TAB%VMR(IP)*1.0E-6 * TAB%QAX(IQ) / GRAVTY / WGTAIR )) END DO ! IQ loop END DO ! IT loop END DO ! IP loop ! END FUNCTION CLCLDU SUBROUTINE COPHDR ( RECNEW ) ! ! VERSION ! 21MAR19 AD F90 conversion ! 21OCT14 AD Original ! ! DESCRIPTION ! Copy comment records from existing .tab file to new file. ! Also optionally inserts a new first record (if RECNEW starts with '!') ! IMPLICIT NONE ! ! ARGUMENTS CHARACTER(80), INTENT(IN) :: RECNEW ! New 1st record for new file ! ! LOCAL VARIABLES INTEGER(I4) :: IOS ! Saved value of IOSTAT for error message CHARACTER(80) :: RECORD ! Header records transferred between files ! ! EXECTUBLE CODE --------------------------------------------------------------- ! REWIND ( LUNOLD, IOSTAT=IOS, ERR=900 ) ! IF ( RECNEW(1:1) .EQ. '!' ) THEN IF ( BINNEW ) THEN WRITE ( LUNNEW, IOSTAT=IOS, ERR=900 ) RECNEW ELSE WRITE ( LUNNEW, '(A)', IOSTAT=IOS, ERR=900 ) RECNEW END IF END IF ! RECORD(1:1) = '!' DO WHILE ( RECORD(1:1) .EQ. '!' ) IF ( BINOLD ) THEN READ ( LUNOLD, IOSTAT=IOS, ERR=900 ) RECORD ELSE READ ( LUNOLD, '(A)', IOSTAT=IOS, ERR=900 ) RECORD END IF IF ( RECORD(1:1) .EQ. '!' ) THEN IF ( BINNEW ) THEN WRITE ( LUNNEW, IOSTAT=IOS, ERR=900 ) RECORD ELSE WRITE ( LUNNEW, '(A)', IOSTAT=IOS, ERR=900 ) RECORD END IF END IF END DO BACKSPACE ( LUNOLD, IOSTAT=IOS, ERR=900 ) RETURN ! 900 CONTINUE WRITE ( *, * ) 'F-COPHDR: Error copying header records, IOSTAT=', IOS STOP ! END SUBROUTINE COPHDR SUBROUTINE EXCGRD ( IA, ID, GRD ) ! ! VERSION ! 29JUL19 AD Original. ! ! DESCRIPTION ! Exclude grid axis point (IA,ID) ! ! ARGUMENTS INTEGER(I4), INTENT(IN) :: IA ! Axis coord of excluded point INTEGER(I4), INTENT(IN) :: ID ! Dimension of excluded point TYPE(GRDTYP), INTENT(INOUT) :: GRD(:,:) ! Tabulation axis data ! ! LOCAL VARIABLES INTEGER(I4) :: IAP ! Axis coord. of previous/current retained point INTEGER(I4) :: IAN ! Axis coord. of next retained point ! ! EXECUTABLE CODE ----------------------------------------------------------------------- ! GRD(IA,ID)%FLG = 0 GRD(IA,ID)%IAP = GRD(IA-1,ID)%IAP ! replace current pt with prev IAP = GRD(IA,ID)%IAP ! Identify prev/next remaining points IAN = GRD(IA,ID)%IAN GRD(IAP:IA-1,ID)%IAN = IAN ! Set prev/next links to skip current pt GRD(IA+1:IAN,ID)%IAP = IAP ! END SUBROUTINE EXCGRD INTEGER(I4) FUNCTION IAXIS ( IX, ID ) ! ! VERSION ! 29JUL19 AD Original ! ! DESCRIPTION ! Return axis coordinate of tabulated point IX on ID axis ! ! ARGUMENTS INTEGER(I4), INTENT(IN) :: IX ! Index of point in tabulation INTEGER(I4), INTENT(IN) :: ID ! Dimension ! ! EXECUTABLE CODE ---------------------------------------------------------------------- ! SELECT CASE ( ID ) CASE ( 1 ) ; IAXIS = 1 + MOD ( IX-1, NQ ) ! q-axis coordinate CASE ( 2 ) ; IAXIS = 1 + MOD ( (IX-1)/NQ, NT ) ! T-axis coordinate CASE ( 3 ) ; IAXIS = 1 + (IX-1)/NTQ ! p-axis coordinate END SELECT ! END FUNCTION IAXIS SUBROUTINE INIGRD ( TAB, GRD, NAD ) ! ! VERSION ! 29JUL19 AD Original. ! ! DESCRIPTION ! Initialise GRD structure ! ! ARGUMENTS TYPE(TABTYP), INTENT(IN) :: TAB ! Tabulated data header TYPE(GRDTYP), INTENT(OUT) :: GRD(:,:) ! Tabulation axis data INTEGER(I4), INTENT(OUT) :: NAD(:) ! No.points in each dimension ! ! LOCAL VARIABLES INTEGER(I4) :: IA ! Counter for axis points INTEGER(I4) :: ID ! Counter for dimensions ! EXECUTABLE CODE ----------------------------------------------------------------------- ! NAD(1) = TAB%NQ NAD(2) = ABS ( TAB%NT ) NAD(3) = TAB%NP ! ! Convert axis values to lin.interp values GRD(1:NAD(1),1)%AXV = LOG ( TAB%QAX ) ! lnq interpolation GRD(1:NAD(2),2)%AXV = TAB%TAX ! T interpolation GRD(1:NAD(3),3)%AXV = LOG ( TAB%PAX ) ! lnp interpolation ! ! Initialise flags and weights for axes DO ID = 1, ND ! Loop over q, T, p axes DO IA = 1, NAD(ID) ! Loop over axis internal points GRD(IA,ID)%FLG = 1 ! flag removable point GRD(IA,ID)%IAP = IA ! location of previous/current point on axis GRD(IA,ID)%IAN = IA + 1 ! location of next retained point on axis END DO ! Special treatment of end points along each table axis GRD(1,ID)%FLG = 2 ! flag unremovable end-axis points GRD(NAD(ID),ID)%FLG = 2 GRD(NAD(ID),ID)%IAN = NAD(ID) END DO ! END SUBROUTINE INIGRD SUBROUTINE INTLKU ( IX, GRD, LKU, LKUINT ) ! ! VERSION ! 29JUL19 AD F90 original. ! ! DESCRIPTION ! Calculate max interpolation error for tabulation point. ! IMPLICIT NONE ! ! ARGUMENTS INTEGER(I4), INTENT(IN) :: IX ! Tab.pt to be interpolated TYPE(GRDTYP), INTENT(IN) :: GRD(:,:) ! Table axis data REAL(R4), INTENT(IN) :: LKU(:,:) ! Original ln(k*U) values for each spec.pt REAL(R4), INTENT(OUT) :: LKUINT(:) ! Interpolated values for each spec.pt ! ! LOCAL CONSTANTS INTEGER(I4), PARAMETER :: MW = 2**ND ! No. of interp. points in q,T,p ! ! LOCAL VARIABLES INTEGER(I4) :: IA ! Axis coordinate of IX for one table dimension INTEGER(I4) :: IA1(MW) ! q-axis coordinates for 3D interp. INTEGER(I4) :: IA2(MW) ! T-axis coordinates for 3D interp. INTEGER(I4) :: IA3(MW) ! p-axis coordinates for 3D interp. INTEGER(I4) :: IAN ! Next axis point for IX INTEGER(I4) :: IAP ! Prev axis point for IX INTEGER(I4) :: IW ! Counter for (MW=8) interpolation weights INTEGER(I4) :: IXP ! Tab.points used for 3D interp. REAL(R4) :: WP ! Interpolation weight for IAP to IX along one axis REAL(R4) :: WGT(MW) ! Interpolation weights for up to 8 points ! ! EXECUTABLE CODE ---------------------------------------------------------------- ! ! Establish set of points required to interpolate this grid point ! q-axis: interp = PPPPNNNN IA = IAXIS(IX,1) IAP = GRD(IA,1)%IAP IAN = GRD(IA,1)%IAN IF ( IAP .EQ. IAN ) THEN ! Only one point on q-axis WP = 1.0 ELSE WP = ( GRD(IAN,1)%AXV - GRD(IA,1)%AXV ) / & ( GRD(IAN,1)%AXV - GRD(IAP,1)%AXV ) END IF IA1(1:4) = IAP WGT(1:4) = WP IA1(5:8) = IAN WGT(5:8) = 1.0 - WP ! T-axis: interp = PPNNPPNN IA = IAXIS(IX,2) IAP = GRD(IA,2)%IAP IAN = GRD(IA,2)%IAN IF ( IAP .EQ. IAN ) THEN ! Only one point on T-axis WP = 1 ELSE WP = ( GRD(IAN,2)%AXV - GRD(IA,2)%AXV ) / & ( GRD(IAN,2)%AXV - GRD(IAP,2)%AXV ) END IF IA2(1:2) = IAP IA2(5:6) = IAP IA2(3:4) = IAN IA2(7:8) = IAN WGT(1:2) = WGT(1:2) * WP WGT(5:6) = WGT(5:6) * WP WGT(3:4) = WGT(3:4) * ( 1.0 - WP ) WGT(7:8) = WGT(7:8) * ( 1.0 - WP ) ! p-axis: interp = PNPNPNPN IA = IAXIS(IX,3) IAP = GRD(IA,3)%IAP IAN = GRD(IA,3)%IAN IF ( IAP .EQ. IAN ) THEN ! Only one point on p-axis WP = 1 ELSE WP = ( GRD(IAN,3)%AXV - GRD(IA,3)%AXV ) / & ( GRD(IAN,3)%AXV - GRD(IAP,3)%AXV ) END IF IA3(1:7:2) = IAP IA3(2:8:2) = IAN WGT(1:7:2) = WGT(1:7:2) * WP WGT(2:8:2) = WGT(2:8:2) * ( 1.0 - WP ) ! LKUINT = 0.0 DO IW = 1, MW IXP = IA1(IW) + (IA2(IW)-1)*NQ + (IA3(IW)-1)*NTQ LKUINT(:) = LKUINT(:) + WGT(IW) * LKU(:,IXP) END DO ! END SUBROUTINE INTLKU SUBROUTINE OPNTAB ( LREAD, LUN, BINFIL ) ! ! VERSION ! 08JUL19 AD Original ! ! DESCRIPTION ! Open .tab file ! IMPLICIT NONE ! ! ARGUMENTS LOGICAL, INTENT(IN) :: LREAD ! T=read/old file, F=write/new file INTEGER(I4), INTENT(IN) :: LUN ! Logical Unit Number LOGICAL, INTENT(INOUT) :: BINFIL ! T=binary file, F=ASCII file ! ! LOCAL VARIABLES INTEGER(I4) :: IOS ! Saved value of IOSTAT CHARACTER(200) :: FILNAM ! Name of .tab file ! ! EXECUTABLE CODE -------------------------------------------------------------- ! IF ( LREAD ) THEN WRITE ( *, '(A)', ADVANCE='NO' ) 'Input file: ' ELSE WRITE ( *, '(A)', ADVANCE='NO' ) 'Output file: ' END IF READ ( *, '(A)' ) FILNAM ! IF ( LREAD ) THEN ! Open existing file for read BINFIL = BINLUT ( LUN, FILNAM ) ! determine if binary or ASCII IF ( BINFIL ) THEN OPEN ( UNIT=LUN, FILE=FILNAM, STATUS='OLD', ACTION='READ', & FORM='UNFORMATTED', IOSTAT=IOS ) ELSE OPEN ( UNIT=LUN, FILE=FILNAM, STATUS='OLD', ACTION='READ', & FORM='FORMATTED', IOSTAT=IOS ) END IF ! ELSE ! Open new file for write IF ( BINFIL .AND. INDEX ( FILNAM, '.asc' ) .GT. 0 ) THEN WRITE ( *, '(A)' ) 'I-OPNTAB: Converting to ASCII file' BINFIL = .FALSE. ELSE IF ( .NOT. BINFIL .AND. INDEX ( FILNAM, '.bin' ) .GT. 0 ) THEN WRITE ( *, '(A)' ) 'I-OPNTAB: Converting to Binary file' BINFIL = .TRUE. ENDIF ! IF ( BINFIL ) THEN OPEN ( UNIT=LUN, FILE=FILNAM, STATUS='UNKNOWN', ACTION='WRITE', & FORM='UNFORMATTED', IOSTAT=IOS ) ELSE OPEN ( UNIT=LUN, FILE=FILNAM, STATUS='UNKNOWN', ACTION='WRITE', & FORM='FORMATTED', IOSTAT=IOS ) ! END IF END IF ! IF ( IOS .NE. 0 ) THEN WRITE (*,*) 'F-OPNTAB: Failed to open file, IOSTAT=', IOS STOP END IF ! END SUBROUTINE OPNTAB SUBROUTINE RWAXIS ( LREAD, LUN, BINFIL, TAB ) ! ! VERSION ! 21MAR19 AD Rewritten in F90 ! 15DEC14 AD Add format identifier record ! 04NOV14 AD Add DV ! 27OCT14 AD Original. ! ! DESCRIPTION ! Read/write axis information from/to .tab file ! IMPLICIT NONE ! ! ARGUMENTS LOGICAL, INTENT(IN) :: LREAD ! T=read header, F=write header INTEGER(I4), INTENT(IN) :: LUN ! Logical Unit Number LOGICAL, INTENT(IN) :: BINFIL ! T=binary file, F=ASCII file TYPE(TABTYP), INTENT(INOUT) :: TAB ! Header structure ! ! LOCAL VARIABLES INTEGER(I4) :: IOS ! Saved value of IOSTAT for error messages REAL(R4) :: RFMT ! File format identifier CHARACTER(80) :: FMTREC ! Text header record containing format CHARACTER(80) :: RECORD ! Text header record ! ! EXECUTABLE CODE -------------------------------------------------------------- ! ! Read in header data from file IF ( LREAD ) THEN RECORD = '!' DO WHILE ( RECORD(1:1) .EQ. '!' ) IF ( BINFIL ) THEN READ ( LUN, IOSTAT=IOS, ERR=900 ) RECORD ELSE READ ( LUN, '(A)', IOSTAT=IOS, ERR=900 ) RECORD END IF END DO ! ! First record contains format identifier variable READ ( RECORD, * ) RFMT IF ( RFMT .GT. 1.0 ) THEN ! Currently only format 1.0 recognised WRITE ( *, * ) & 'F-RWAXIS: unrecognised LUT file format, version=', RFMT STOP END IF TAB%FMT = RFMT ! IF ( BINFIL ) THEN READ ( LUN, IOSTAT=IOS, ERR=900 ) RECORD ELSE READ ( LUN, '(A)', IOSTAT=IOS, ERR=900 ) RECORD END IF ! TAB%MOL = RECORD(1:5) READ ( RECORD(6:), * ) & TAB%NV, TAB%V1, TAB%V2, TAB%DV, TAB%NX, TAB%NP, TAB%NT, TAB%NQ ! ALLOCATE ( TAB%PAX(TAB%NP), TAB%TEM(TAB%NP), TAB%VMR(TAB%NP), & TAB%TAX(ABS(TAB%NT)), TAB%QAX(TAB%NQ) ) IF ( BINFIL ) THEN READ ( LUN, IOSTAT=IOS, ERR=900 ) TAB%PAX READ ( LUN, IOSTAT=IOS, ERR=900 ) TAB%TEM READ ( LUN, IOSTAT=IOS, ERR=900 ) TAB%VMR READ ( LUN, IOSTAT=IOS, ERR=900 ) TAB%TAX READ ( LUN, IOSTAT=IOS, ERR=900 ) TAB%QAX ELSE READ ( LUN, *, IOSTAT=IOS, ERR=900 ) TAB%PAX READ ( LUN, *, IOSTAT=IOS, ERR=900 ) TAB%TEM READ ( LUN, *, IOSTAT=IOS, ERR=900 ) TAB%VMR READ ( LUN, *, IOSTAT=IOS, ERR=900 ) TAB%TAX READ ( LUN, *, IOSTAT=IOS, ERR=900 ) TAB%QAX END IF ! ! Write header data to file ELSE WRITE ( FMTREC, '(F10.2)' ) TAB%FMT WRITE ( RECORD, '(A5,I10,2F12.4,F8.4,4I8)' ) & TAB%MOL, TAB%NV, TAB%V1, TAB%V2, TAB%DV, TAB%NX, TAB%NP, TAB%NT, TAB%NQ IF ( BINFIL ) THEN WRITE ( LUN, IOSTAT=IOS, ERR=900 ) FMTREC WRITE ( LUN, IOSTAT=IOS, ERR=900 ) RECORD WRITE ( LUN, IOSTAT=IOS, ERR=900 ) TAB%PAX WRITE ( LUN, IOSTAT=IOS, ERR=900 ) TAB%TEM WRITE ( LUN, IOSTAT=IOS, ERR=900 ) TAB%VMR WRITE ( LUN, IOSTAT=IOS, ERR=900 ) TAB%TAX WRITE ( LUN, IOSTAT=IOS, ERR=900 ) TAB%QAX ELSE WRITE ( LUN, '(A)', IOSTAT=IOS, ERR=900 ) FMTREC WRITE ( LUN, '(A)', IOSTAT=IOS, ERR=900 ) RECORD WRITE ( LUN, *, IOSTAT=IOS, ERR=900 ) TAB%PAX WRITE ( LUN, *, IOSTAT=IOS, ERR=900 ) TAB%TEM WRITE ( LUN, *, IOSTAT=IOS, ERR=900 ) TAB%VMR WRITE ( LUN, *, IOSTAT=IOS, ERR=900 ) TAB%TAX WRITE ( LUN, *, IOSTAT=IOS, ERR=900 ) TAB%QAX END IF END IF ! Normal exit RETURN ! 900 CONTINUE WRITE ( *, * ) 'F-RWAXIS: I/O error, IOSTAT=', IOS STOP ! END SUBROUTINE RWAXIS SUBROUTINE RWTAB ( LREAD, LUN, BINFIL, V, K ) ! ! VERSION ! 21MAR19 AD Orginal. ! ! DESCRIPTION ! Read/write data record from/to .tab file ! IMPLICIT NONE ! ! ARGUMENTS LOGICAL, INTENT(IN) :: LREAD ! T=read header, F=write header INTEGER(I4), INTENT(IN) :: LUN ! Logical Unit Number LOGICAL, INTENT(IN) :: BINFIL ! T=binary file, F=ASCII file REAL(R8), INTENT(INOUT) :: V ! Wavenumber [cm-1] REAL(R4), INTENT(INOUT) :: K(:) ! ln(k) absorp.coeff. ! ! LOCAL VARIABLES INTEGER(I4) :: IOS ! Saved value of IOSTAT for error messages ! ! EXECUTABLE CODE -------------------------------------------------------------- ! IF ( LREAD ) THEN IF ( BINFIL ) THEN READ ( LUN, IOSTAT=IOS, ERR=900 ) V, K ELSE READ ( LUN, *, IOSTAT=IOS, ERR=900 ) V, K END IF ELSE IF ( BINFIL ) THEN WRITE ( LUN, IOSTAT=IOS, ERR=900 ) V, K ELSE WRITE ( LUN, *, IOSTAT=IOS, ERR=900 ) V, K END IF END IF RETURN 900 CONTINUE WRITE ( *, * ) 'F-RWTAB: Error reading/writing .tab data record, IOSTAT=',IOS STOP END SUBROUTINE RWTAB END MODULE TABCMP_X_MOD PROGRAM TABCMP_X ! ! DESCRIPTION ! Main program. ! USE TABCMP_X_MOD ! IMPLICIT NONE ! ! LOCAL VARIABLES LOGICAL :: OFFTAB ! T=temperature axis relative to TEMTAB,F=absolute LOGICAL :: RETAIN ! T=tabulated grid point retained, F=removed INTEGER(I4) :: IA ! Axis coordinate of IX for one table dimension INTEGER(I4) :: IAMIN ! Axis point removed with minimum interp.error INTEGER(I4) :: ID ! Dimension counter INTEGER(I4) :: IDMIN ! Dimension with pt removed with min. interp.error INTEGER(I4) :: IOS ! Saved value of IOSTAT for error messages INTEGER(I4) :: IP ! Counter for pressure axis points INTEGER(I4) :: IPN ! Counter for new pressure axis points INTEGER(I4) :: IQ ! Counter for VMR axis points INTEGER(I4) :: IQN ! Counter for new VMR axis points INTEGER(i4) :: IT ! Counter for temperature axis points INTEGER(i4) :: ITN ! Counter for new temperature axis points INTEGER(I4) :: IV ! Counter for spectral points INTEGER(I4) :: IX ! Counter for tabulated points INTEGER(I4) :: NAD(ND) ! No. axis points in each table dimensions q,T,p INTEGER(I4) :: NPNEW ! No. modified pressure axis points INTEGER(I4) :: NQNEW ! No. modified VMR axis points INTEGER(I4) :: NTNEW ! No. modified temperature axis points INTEGER(I4) :: NV ! No. spectral points in original file INTEGER(I4) :: NX ! No. original tabulated k(p,T,q) axis coord points INTEGER(I4) :: NXNEW ! No. modified, tabulated k(p,T,q) points REAL(R4) :: ERRLIM ! User-defined maximum absorption error REAL(R4) :: ERRMIN ! Minimum error for removing any point REAL(R4) :: ERRX ! Maximum interpolation error if axis point removed REAL(R8) :: V ! Individual spectral axis value CHARACTER(80) :: RECORD ! Record read from file TYPE(TABTYP) :: TAB ! Tabulated data header INTEGER(I4), ALLOCATABLE :: IXLIST(:) ! Original indices of ret. pts REAL(R4), ALLOCATABLE :: DA(:) ! [NV] Difference in absorption 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 :: KU(:) ! [NV] (k*U) values at grid pt. REAL(R4), ALLOCATABLE :: KUINT(:) ! [NV] interp.(k*U) values at grid pt. REAL(R4), ALLOCATABLE :: LKU(:,:) ! [NV,NX] ln(k*U) values for spc.pt REAL(R4), ALLOCATABLE :: LKUINT(:) ! [NV] Interp. ln(k) at grid point REAL(R4), ALLOCATABLE :: LNK(:,:) ! [NV,NX] ln(k) values REAL(R4), ALLOCATABLE :: PSAV(:) ! Saved pressure axis REAL(R4), ALLOCATABLE :: TSAV(:) ! Saved temperature profile or axis REAL(R4), ALLOCATABLE :: VSAV(:) ! Saved VMR profile or axis REAL(R8), ALLOCATABLE :: VTAB(:) ! [NV] List of wno values [cm-1] TYPE(GRDTYP), ALLOCATABLE :: GRD(:,:) ! [NMAX,ND] Grid axis data TYPE(GRDTYP), ALLOCATABLE :: GRDSAV(:,:) ! [NMAX,ND] Saved GRD ! ! EXECUTABLE CODE -------------------------------------------------------------- ! WRITE ( *, '(A)' ) 'R-TABCMP_X: Running TABCMP_X v' // VERSID ! ! Open original .tab file CALL OPNTAB ( .TRUE., LUNOLD, BINOLD ) ! ! Read header data into TAB CALL RWAXIS ( .TRUE., LUNOLD, BINOLD, TAB ) ! ! Calculate ln(abs.amt) in each 'cell' of tabulation NX = TAB%NX ALLOCATE ( LDU(NX) ) LDU = CLCLDU ( TAB ) ! WRITE ( *, '(A)', ADVANCE='NO' ) 'Max absorption error (eg 0.001): ' READ ( *, * ) ERRLIM ! BINNEW = BINOLD CALL OPNTAB ( .FALSE., LUNNEW, BINNEW ) ! ! Read in data and convert to cell optical thickness NV = TAB%NV ALLOCATE ( VTAB(NV), KU(NV), KUINT(NV), LKUINT(NV), DA(NV), & LKU(NV,NX), LNK(NV,NX), KTAB(NX) ) WRITE ( *, '(A)' ) 'I-TABCMP_X: Reading in LUT ...' DO IV = 1, NV CALL RWTAB ( .TRUE., LUNOLD, BINOLD, VTAB(IV), KTAB ) LKU(IV,:) = KTAB + LDU END DO ! OFFTAB = TAB%NT .LT. 0 NT = ABS ( TAB%NT ) NP = TAB%NP NQ = TAB%NQ NTQ = NT * NQ ! NMAX = MAX ( NP, NT, NQ ) ALLOCATE ( GRD(NMAX,ND), GRDSAV(NMAX,ND) ) CALL INIGRD ( TAB, GRD, NAD ) ! WRITE ( *, '(A,I5)' ) 'I-TABCMP_X: Starting to remove points. ' & // 'Initial No.=', NX ! ! Loop back to here to find next axis point for removal DO ERRMIN = 2*ERRLIM ! Calculate interpolation error for each axis point DO ID = 1, ND ! Loop over table dimensions q,T,p DO IA = 2, NAD(ID)-1 ! Loop over internal axis points IF ( GRD(IA,ID)%FLG .EQ. 1 ) THEN ! Try removing axis point ! Temporarily remove axis point GRDSAV = GRD CALL EXCGRD ( IA, ID, GRD ) ERRX = 0.0 DO IX = 1, NX ! Loop over all grid points CALL INTLKU ( IX, GRD, LKU, LKUINT ) KU = EXP ( LKU(:,IX) ) KUINT = EXP ( LKUINT ) DA = ABS ( EXP ( -KU ) - EXP ( -KUINT ) ) ERRX = MAX ( ERRX, MAXVAL(DA) ) END DO ! IX loop IF ( ERRX .LT. ERRMIN ) THEN IAMIN = IA IDMIN = ID ERRMIN = ERRX END IF ! Restore tabulation axis point GRD = GRDSAV END IF ! If IA is remaining axis point END DO ! IA loop END DO ! ID loop IF ( ERRMIN .GT. ERRLIM ) EXIT SELECT CASE ( IDMIN ) CASE ( 1 ) WRITE ( *, '(A,F10.2,A,F10.7)' ) 'Removing q-axis value=', & TAB%QAX(IAMIN), ' [%], ERR=', ERRMIN CASE ( 2 ) WRITE ( *, '(A,F10.2,A,F10.7)' ) 'Removing T-axis value=', & TAB%TAX(IAMIN), ' [K], ERR=', ERRMIN CASE ( 3 ) WRITE ( *, '(A,ES10.3,A,F10.7)' ) 'Removing p-axis value=', & TAB%PAX(IAMIN), ' [mb], ERR=', ERRMIN END SELECT CALL EXCGRD ( IAMIN, IDMIN, GRD ) END DO ! NXNEW = 0 ALLOCATE ( IXLIST(NX) ) DO IX = 1, NX RETAIN = .TRUE. DO ID = 1, ND IA = IAXIS(IX,ID) IF ( GRD(IA,ID)%FLG .EQ. 0 ) RETAIN = .FALSE. END DO IF ( RETAIN ) THEN NXNEW = NXNEW + 1 IXLIST(NXNEW) = IX END IF END DO ! WRITE ( *, '(A,I9,A,I9,A,F5.1,A)' ) 'I-TABCMP_X: Summary: Orig Npts=', NX, & ' New Npts=', NXNEW, ' Red.Factor=', & 100 * FLOAT ( NXNEW ) / FLOAT ( NX ), '%:' WRITE ( *, '(A)' ) 'I-TABCMP_X: Writing new .tab file' ! RECORD = '! Reduced grid using program TABCMP_X v.'//VERSID CALL COPHDR ( RECORD ) ! CALL RWAXIS ( .TRUE., LUNOLD, BINOLD, TAB ) ! NPNEW = COUNT ( GRD(:,3)%FLG .NE. 0 ) IF ( NPNEW .NE. NP ) THEN ALLOCATE ( PSAV(NP), TSAV(NP), VSAV(NP) ) PSAV = TAB%PAX TSAV = TAB%TEM VSAV = TAB%VMR DEALLOCATE ( TAB%PAX, TAB%TEM, TAB%VMR ) ALLOCATE ( TAB%PAX(NPNEW), TAB%TEM(NPNEW), TAB%VMR(NPNEW) ) IPN = 0 DO IP = 1, NP IF ( GRD(IP,3)%FLG .NE. 0 ) THEN IPN = IPN + 1 TAB%PAX(IPN) = PSAV(IP) TAB%TEM(IPN) = TSAV(IP) TAB%VMR(IPN) = VSAV(IP) END IF END DO TAB%NP = NPNEW DEALLOCATE ( PSAV, TSAV, VSAV ) END IF ! NTNEW = COUNT ( GRD(:,2)%FLG .NE. 0 ) IF ( NTNEW .NE. NT ) THEN ALLOCATE ( TSAV(NT) ) TSAV = TAB%TAX DEALLOCATE ( TAB%TAX ) ALLOCATE ( TAB%TAX(NTNEW) ) ITN = 0 DO IT = 1, NT IF ( GRD(IT,2)%FLG .NE. 0 ) THEN ITN = ITN + 1 TAB%TAX(ITN) = TSAV(IT) END IF END DO TAB%NT = NTNEW IF ( OFFTAB ) TAB%NT = -NTNEW DEALLOCATE ( TSAV ) END IF ! NQNEW = COUNT ( GRD(:,1)%FLG .NE. 0 ) IF ( NQNEW .NE. NQ ) THEN ALLOCATE ( VSAV(NQ) ) VSAV = TAB%QAX DEALLOCATE ( TAB%QAX ) ALLOCATE ( TAB%QAX(NQNEW) ) IQN = 0 DO IQ = 1, NQ IF ( GRD(IQ,1)%FLG .NE. 0 ) THEN IQN = IQN + 1 TAB%QAX(IQN) = VSAV(IQ) END IF END DO TAB%NQ = NQNEW DEALLOCATE ( VSAV ) END IF ! IF ( NXNEW .NE. NPNEW * NTNEW * NQNEW ) STOP 'F-TABCMP_X: Logical error' TAB%NX = NXNEW ! CALL RWAXIS ( .FALSE., LUNNEW, BINNEW, TAB ) ! DO IV = 1, NV CALL RWTAB ( .TRUE., LUNOLD, BINOLD, V, KTAB ) KTAB(1:NXNEW) = KTAB(IXLIST(1:NXNEW)) CALL RWTAB ( .FALSE., LUNNEW, BINNEW, V, KTAB(1:NXNEW) ) END DO ! CLOSE ( LUNOLD, IOSTAT=IOS, ERR=900 ) CLOSE ( LUNNEW, IOSTAT=IOS, ERR=900 ) ! 900 CONTINUE IF ( IOS .EQ. 0 ) THEN WRITE ( *, '(A)' ) 'R-TABCMP_X: Successful completion' ELSE WRITE ( *, * ) 'F-TABCMP_X: I/O error, IOSTAT=', IOS STOP END IF ! END PROGRAM TABCMP_X