SVD-DECOMPRESSION FOR TABULATED K-COEFFICIENTS =============(PO-TN-OXF-GS-0011)============== A. Dudhia 08-DEC-04 08-DEC-04 Change MWCODE from C*6 to C*8 13-NOV-97 Add note on format of MWCODE+ID+TAB for binary files Renumber notes 01-NOV-97 Change (p,T) interpolation: interpolate in ln(k) for all types of tabulation (LIN,LOG,4RT) 31-OCT-97 Document changes from Progress Meeting (24-OCT-97) Define ordering of NX dimension of K Correct indices in U(NV,*) record 17-JUL-97 Correction: define k in units of (m^2/mole). 17-JUL-97 Add Tabulation code to header: C*3 TAB Define units of k (moles/m^2) 16-JUL-97 Add MWCODE and GAS ID record to LUT file header Absorption coefficient tabulated as ln(k) instead of k directly. Remove comments on using relative T(p) profile (not applicable) Add comments on handling p,T at table edges 10-MAR-97 1st Draft 1. FILE FORMAT ============== ASCII or Binary files, one file per microwindow per absorber. --- ! An arbitrary number of initial comments records, starting with `!' MWCODE ID TAB [ MW id, Gas id, Tabulation id - Note 4] NL NV V1 DV NP P1 DP NT T1 DT [ table dimensions ] U(1,1) U(1,2) ... U(1,NL) [ NL elements in each record ] U(2,1) U(2,2) ... U(2,NL) ... ... U(NV,1) U(NV,2) ... U(NV,NL) [total of NV records for U-matrix] K(1,1) K(2,1) ... K(NL,1) [ NL elements in each record ] K(1,2) K(2,2) ... K(NL,2) ... ... K(1,NX) K(2,NX) ... K(NL,NX) [total of NX=NP*NT records for K-matrix] --- where: MWCODE (C*8) is the microwindow identifier (eg 'O3__0001'), followed by 1 space ID (C*2) is the HITRAN code for the modeclue (eg ' 3' for ozone), then 1 space TAB (C*3) is the Tabulation code: 'LIN', 'LOG', '4RT' (see note 3) NL (I*4) Number of basis vectors (typically 10) NV (I*4) Number of wavenumber points in microwindow (typically 2000) V1 (R*4) Microwindow lower wavenumber boundary [cm-1] DV (R*4) Spacing [cm-1] for U-matrix tabulation (normally 0.0005 cm-1) NP (I*4) Number of -ln(pressure) tabulation points (typically 25) P1 (R*4) Lowest -ln(pressure) [ p in mb] (typically -6.0) DP (R*4) Spacing of -ln(pressure) tabulation (typically 1.0) NT (I*4) Number of temperature tabulation points (typically 10) T1 (R*4) Lowest tabulated temperature [K] (typically 180 K) DT (R*4) Spacing of temperature tabulation [K] (typically 15 K) U(i,j) (R*4) are the U-matrix elements, i=1,NV, j=1,NL K(j,k) (R*4) are the K-matrix elements, j=1,NL, k=1,NX where NX=NP*NT. K(j,1)=K(j,P1,T1); K(j,2)=K(j,P1+DP,T1) ... K(j,NP+1)=K(j,P1,T1+DT) --- Notes 1. Although it would be quicker to read in U if it were transposed (ie if each record contained U(i,j),i=1,NV instead of U(i,j),j=1,NL), it is assumed that this is outweighed by the convenience of having each data record in the file the same length (NL elements) for both the U and K matrices. 2. The (p,T) values are treated as a single dimension (1:NX) to allow minimum array storage within the code (ie so NP and NT do not have to be dimensioned separately in the array declarations, just a maximum value for NX). 3. The actual compressed/reconstructed tables can either represent the absorption coefficient directly (TAB='LIN') or some function: TAB='LOG' implies tabulation is of ln(k) TAB='4RT' implies tabulation is of SQRT(SQRT(k)) The reading program should check that tabulation ID corresponds to one of the functions coded in the decompression stage. 4. For Binary Files, the MWCODE+ID+TAB will be written as a single C*13 string equivalent to (A6,X,I2,X,A3) 2.DECOMPRESSION =============== Purpose: return the absorption coefficient vector KABS(1:NV) across the whole microwindow for path conditions (p,T), where p is -ln(pressure/mb), T is the temperature in K, and KABS is in units of m^2/mole. It is assumed that the compressed data are stored in arrays U(NV,NL), K(NL,NX), with additional indices (not shown here) to identify the appropriate microwindow/absorber combination required. Variables beginning I,J,N are integers, others are real. 2.1 Calculation of Interpolation weights for ln(k) in (p,T) domain ------------------------------------------------------------------ Interpolation points and weights in -ln(pressure) axis XP = ( P - P1 ) / DP + 1.0 XP = MIN ( MAX ( 1.0, XP ), FLOAT ( NP ) ) ! limit to 1:XP IP = MIN ( INT ( XP ), NP-1 ) ! limit to 1:NP-1 DP = XP - FLOAT ( IP ) Interpolation points and weights in Temperature axis XT = ( T - T1 ) / DT + 1.0 XT = MIN ( MAX ( 1.0, XT ), FLOAT ( NT ) ) ! limit to 1:XT IT = MIN ( INT ( XT ), NT-1 ) ! limit to 1:NT-1 DT = XT - FLOAT ( IT ) Indices and weights in 1:NX dimension for (p,T) interpolation II = IP + NP * ( IT - 1 ) ! low -lnp, low temperature JI = II + 1 ! high -lnp, low temperature IJ = II + NP ! low -lnp, high temperature JJ = IJ + 1 ! high -lnp, high temperature WII = ( 1.0 - DP ) * ( 1.0 - DT ) WJI = DP * ( 1.0 - DT ) WIJ = ( 1.0 - DP ) * DT WJJ = DP * DT Notes: 5. XP and XT are limited to the ranges 1:NP, 1:NT to ensure there is no extrapolation in cases where the required p,T are outside the tabulated range (in this case the edge values are used). 6. IP and IT are limited to the ranges 1:NP-1, 1:NT-1 to ensure that when XP=NP, XT=NT the interpolation does not attempt to access undefined memory elements. 2.2 Expansion from basis vectors -------------------------------- DO IV = 1, NV ! loop over all microwindow spectral grid pts KII = 0.0 ! Initialise to zero KIJ = 0.0 KJI = 0.0 KJJ = 0.0 DO IL = 1, NL ! loop over all basis vectors KII = KII + U(IV,IL) * K(IL,II) ! Reconstruct tabulated function at KIJ = KIJ + U(IV,IL) * K(IL,IJ) ! the four interpolation points KJI = KJI + U(IV,IL) * K(IL,JI) KJJ = KJJ + U(IV,IL) * K(IL,JJ) END DO IF ( USELOG ) THEN ! If tabulated ln(k) KABS(IV) = EXP ( WII*KII + WIJ*KIJ + WJI*KJI + WJJ*KJJ ) ELSE ! If tabulated k or k**0.25 KABS(IV) = EXP ( WII * LOG ( MAX ( KII, KMIN ) ) + & WIJ * LOG ( MAX ( KIJ, KMIN ) ) + & WJI * LOG ( MAX ( KJI, KMIN ) ) + & WJJ * LOG ( MAX ( KJJ, KMIN ) ) ) IF ( USE4RT ) KABS(IV) = KABS(IV)**4 ! Tabulated k**0.25 END IF END DO Notes 7. The above assumes that the absorption coefficients are required on the same wavenumber grid as the U tabulation, so no interpolation is performed in the wavenumber dimension. 8. The interpolation in (p,T) is always carried out in ln(k) since ln(k) is generally a linear function of -ln(p) (in the Lorentz limit). 9. KMIN (eg set equal to 1.0E-38) is necessary to ensure that ln(k) returns a reasonable value when the reconstructed absorption coefficients are close to zero or negative (possible with k or k**0.25 tabulations) 10.Since ln(k**0.25) = 0.25*ln(k), interpolation in (p,T) domain for 'LIN' and '4RT' reconstructions is the same, with the expansion from k**0.25 to k left until the last step when KABS is calculated. --------