PROGRAM RTSPEC C Program for computing thermal atmospheric radiative transfer with C a multiple layer cloud in an atmosphere. Outputs a radiance spectrum, C suitably averaged. Three approximate solutions for the scattering C radiative transfer in the cloud layers are available: Single scattering C approximation, Eddington second approximation and a hybrid single C scattering/Eddington approximation. Comparison with an "exact" C doubling-adding model for ice clouds has shown that the hybrid model C is superior in the mid IR (300-3000 cm^-1) but not always in the C submm (0-50 cm^-1) (for particle sizes Dme=30 to 300 um). C Reference: Deeter, M. N. and K. F. Evans, 1998: C A hybrid Eddington-single scattering radiative transfer model C for computing radiances from thermally emitting atmospheres. C J. Quant. Spectosc. Radiat. Transfer, 60, 635-648. C The scattering approximations do radiative transfer for a single C homogeneous layer. The input to them is the incident radiation on C the layer. The incident radiance is computed by a nonscattering C radiative transfer integration of cloud and gaseous absorption profiles C that are read in from files. The scattering routines output C upwelling radiance at the top of the layer, and this is propagated C to the observation level for output. The multiple layers are dealt C with by computing the appropriate incident radiances/fluxes for each C layer, and then running the single layer routines for each layer from C the bottom up. For upward looking geometry, the downwelling radiance C is computed by just switching the incident radiation and the top/bottom C temperature of the layer, thus the "upwelling" radiance that the subroutine C outputs is actually the downwelling. Note: It is not necesary to switch C the "phi" functions since "+/-" doesn't mean up/down just forward/back. C For multiple layers, the single layer routines must be run from top down. C The cloud boundaries must be at levels in the input absorption profile. C The surface is assumed to reflect specularly with the specified C emissivity. The emissivity is given for the starting and ending C wavenumbers and linearly interpolated between. Multiple reflection C between the cloud (scattering layer) and surface is ignored. The C input includes the cloud properties, which are ice water path (IWP) C and particle size (Dme) for each layer. There may be multiple input C scattering table files for different cloud types. C The output is upwelling or downwelling brightness temperature or C radiance as a function of wavenumber at the desired observation level. C The radiance is calculated at the absorption file wavenumbers. C For output the radiance may be averaged with either a triangular C or rectangular bandpass or convolved with sinc or sinc squared. C The width of the rectangular bandpass or the full width half max C (FHHM) of the triangular bandpass is equal to the output wavenumber C spacing. The output wavenumbers, which are the center of the C bandpasses, are on the regular grid specified by the user. The C convolution averaging (for a Fourier Transform Spectrometer) is C done by FFTing and filtering with a rectangular (sinc), triangular C (sinc squared), or sinc^2 (out to first zero) windowing function. C Multiple spectral segments may be output, with potentially different C absorption profile input files, output units, and output bandpass C averaging for each one. Multiple "channels" may be output by making C the output spacing for each segment equal to zero, setting the C channel width, and using rectangular averaging. The channel width C extends equally from each side of the output wavenumber. If the channel C width is one absorption file wavenumber spacing or less then C monochromatic radiative transfer is done. A double sideband receiver C may be emulated, in which case two rectangular bandpasses equally C spaced from the receiver center frequency are simulated. C C The input ascii absorption file has the following format: C 1) header line C 2) starting, ending, delta, and number of wavenumbers in file C 3) number of levels (1 + number of atmosphere layers) C 4) heights of levels (from top to bottom) (in km) C 5) temperatures at levels (in Kelvin) C wavenumber (cm^-1) and layer optical depths for each Nlev-1 layers C . . . C A binary format absorption file may be input for faster loading C (the BINARYABSFILE logical is set from user input; see C READ_ABS_PROFILE subroutine). C C The input scattering tables are in the format produced by C sscatmie.f. This format also accomodates single scattering C information for randomly nonspherical oriented particles. C C Author: Frank Evans University of Colorado July 2001 C Co-authors: C Aaron Evans (multilayer code, spectral segments, double sidebands, C upward looking geometry, output convolution) C Merritt Deeter (single scattering layer routines) IMPLICIT NONE C The scattering tables are read in with READ_SSCATTAB. The scattering C table is 3D: wavenumber, particle size, and viewing angle. C Scattering table variables: C MUTAB is view angle values (cosine zenith), C DMETAB is particle size values (median mass diameter, micron), C WAVETAB is wavenumber values (cm^-1). C MUINC(2) are the mu values of the two incident angles C TABEXTINCT is extinction, TABSSALB is single scattering albedo, C TABASYM is the asymmetry parameter C TABPHI??? are phase function info for incident directions INTEGER MAXTAB, MAXGRID, MAXSCAT PARAMETER (MAXTAB=10*25*500, MAXGRID=10000, MAXSCAT=10) INTEGER NMUOBS(MAXSCAT), NDME(MAXSCAT), NWAVETAB(MAXSCAT) REAL MUTAB(MAXGRID,MAXSCAT) REAL DMETAB(MAXGRID,MAXSCAT), WAVETAB(MAXGRID,MAXSCAT) REAL MUINC(2) REAL TABEXTINCT(MAXTAB,MAXSCAT), TABSSALB(MAXTAB,MAXSCAT) REAL TABASYM(MAXTAB,MAXSCAT) REAL TABPHI1UP(MAXTAB,MAXSCAT), TABPHI1DN(MAXTAB,MAXSCAT) REAL TABPHI2UP(MAXTAB,MAXSCAT), TABPHI2DN(MAXTAB,MAXSCAT) C Radiative transfer variables: INTEGER MAXABSNU, MAXNZ, MAXSPEC PARAMETER (MAXABSNU=10000, MAXNZ=100, MAXSPEC=10) INTEGER NUMSPEC, NSCATTAB, NCLDLAY, NABSNU, NLEV INTEGER ICLDTOP, ICLDBOT, IOBS, ISCATTAB(MAXNZ) INTEGER I, J, J0, L, M, JNU1, JNU2, NBAND, IBAND INTEGER NWAVENO(MAXSPEC), NOUTNU LOGICAL TWOSIDEBAND, BINARYABSFILE REAL ZOBS, MUOBS(MAXSPEC), ZTOP, IWP(MAXNZ), DME(MAXNZ) REAL SFCTEMP, SFCEMIS, SFCEMIS1(MAXSPEC), SFCEMIS2(MAXSPEC) REAL RADOBS(2,MAXSPEC,MAXABSNU), OUTRAD(MAXABSNU) REAL HEIGHT(MAXNZ), TEMP(MAXNZ), ABSPROF(MAXNZ,MAXABSNU) REAL*8 ABSNU1, ABSNU2, ABSDELNU REAL*8 OUTNU1(MAXSPEC), OUTNU2(MAXSPEC), OUTDELNU(MAXSPEC) REAL*8 BANDCENTER(MAXSPEC), BANDOFFSET(MAXSPEC) REAL*8 OUTNU(MAXABSNU), WAVENO(2,MAXSPEC,MAXABSNU), XO CHARACTER*100 SCATFILE(MAXSCAT), ABSFILE(MAXSPEC), OUTPUTFILE CHARACTER*1 RTMODEL, ABSTYPE, OUTUNITS(MAXSPEC) CHARACTER*2 OUTAVERAGE(MAXSPEC) C Read in all the user input CALL INPUT_PARAMETERS (NUMSPEC, MAXSPEC, TWOSIDEBAND, OUTNU1, . OUTNU2, OUTDELNU, BANDCENTER, BANDOFFSET, MUOBS, . SFCEMIS1, SFCEMIS2, OUTUNITS, OUTAVERAGE, XO, . ABSFILE, NSCATTAB, MAXSCAT, SCATFILE, ZOBS, . SFCTEMP, NCLDLAY, MAXNZ, IWP, DME, ISCATTAB, ZTOP, . RTMODEL, ABSTYPE, BINARYABSFILE, OUTPUTFILE) C Read in scattering tables DO I = 1, NSCATTAB CALL READ_SSCATTAB (SCATFILE(I), MAXTAB, MAXGRID, . NMUOBS(I), MUTAB(1,I), NDME(I), DMETAB(1,I), . NWAVETAB(I), WAVETAB(1,I), . MUINC, TABEXTINCT(1,I), TABSSALB(1,I), TABASYM(1,I), . TABPHI1UP(1,I), TABPHI1DN(1,I), . TABPHI2UP(1,I), TABPHI2DN(1,I)) IF ((ABS(MUINC(1)-0.2113) .GT. 0.001) .OR. . (ABS(MUINC(2)-0.7887) .GT. 0.001)) . STOP 'RTSPEC: Coded for incident mu=0.2113,0.7887' ENDDO C Read in the first absorption profile file C Absorption file may be binary or ascii CALL READ_ABS_PROFILE (ABSFILE(1), BINARYABSFILE, . ABSNU1, ABSNU2, ABSDELNU, NABSNU, . MAXNZ, MAXABSNU, NLEV, HEIGHT, TEMP, ABSPROF) IF (OUTDELNU(1) .LT. ABSDELNU) THEN WRITE (*,*) 'RTSPEC: output wavenumber resolution less than', . ' in absorption file' STOP ENDIF C Find the levels for top of cloud and observation level ICLDTOP = 0 IOBS = 0 DO I = 1, NLEV IF (ABS(HEIGHT(I)-ZTOP) .LT. 0.010) ICLDTOP = I IF (ABS(HEIGHT(I)-ZOBS) .LT. 0.010) IOBS = I ENDDO ICLDBOT = ICLDTOP + NCLDLAY IF (ICLDTOP.EQ.0 .OR. ICLDBOT.EQ.0 .OR. IOBS.EQ.0) THEN WRITE (*,*) 'RTSPEC: Observer or cloud top height does', . ' not match absorption file levels' STOP ENDIF IF(SFCTEMP.EQ.-1) SFCTEMP = TEMP(NLEV) C if using double side band need to loop twice over RT/chan IF(.NOT.TWOSIDEBAND) NBAND = 1 IF(TWOSIDEBAND) NBAND = 2 C Loop over each segment of spectra DO M = 1,NUMSPEC C Read in gas absorption file for each segment IF (M .GT. 1 .AND. ABSFILE(M) .NE. ABSFILE(M-1)) THEN CALL READ_ABS_PROFILE (ABSFILE(M), BINARYABSFILE, . ABSNU1, ABSNU2, ABSDELNU, NABSNU, . MAXNZ, MAXABSNU, NLEV, HEIGHT, TEMP, ABSPROF) ENDIF C loop if have double side band DO IBAND = 1,NBAND C Need to find the Tb for each side band IF(TWOSIDEBAND) THEN OUTNU1(M) = BANDCENTER(M)+(2.*DBLE(IBAND)-3.)* . BANDOFFSET(M) OUTNU2(M) = OUTNU1(M) ENDIF C Get spectral range in absorption file for this spectral segment IF (OUTAVERAGE(M).EQ.'R'.OR.OUTAVERAGE(M).EQ.'r') THEN JNU1 = INT((OUTNU1(M)-0.5D0*OUTDELNU(M) - ABSNU1) . /ABSDELNU + 0.9999D0) + 1 JNU2 = INT((OUTNU2(M)+0.5D0*OUTDELNU(M) - ABSNU1) . /ABSDELNU + 0.0001D0) + 1 ELSE JNU1 = INT((OUTNU1(M)-OUTDELNU(M) - ABSNU1) . /ABSDELNU + 0.9999D0) + 1 JNU2 = INT((OUTNU2(M)+OUTDELNU(M) - ABSNU1) . /ABSDELNU + 0.0001D0) + 1 ENDIF IF (JNU1 .LE. 0) THEN PRINT *, ABSNU1,OUTNU1(M),OUTDELNU(M) WRITE (*,*) . 'RTSPEC: lower wavenumber out of absfile range.', . ' Segment:',M STOP ENDIF IF (JNU2 .GT. NABSNU) THEN WRITE (*,*) . 'RTSPEC: upper wavenumber out of absfile range.', . ' Segment:',M STOP ENDIF C Loop over a segment of wavenumbers input in the absorption file NWAVENO(M) = JNU2-JNU1+1 DO J = JNU1, JNU2 J0 = J - JNU1 + 1 WAVENO(IBAND,M,J0) = ABSNU1 + (J-1)*ABSDELNU C Interpolate the emissivity linearly along segement IF(OUTNU1(M).EQ.OUTNU2(M)) THEN SFCEMIS = SFCEMIS1(M) ELSE SFCEMIS = (SFCEMIS2(M)-SFCEMIS1(M))/(OUTNU2(M)-OUTNU1(M))* . (WAVENO(IBAND,M,J0)-OUTNU1(M)) + SFCEMIS1(M) ENDIF C Calculate the upwelling observed radiance for this wavenumber C Store radiance at absorption file wavenumber resolution CALL COMPUTE_RADIATIVE_TRANSFER (RTMODEL, . MUOBS(M), IOBS, SNGL(WAVENO(IBAND,M,J0)), . NLEV, TEMP, ABSPROF(1,J), SFCTEMP, SFCEMIS, . NCLDLAY, ICLDTOP, ICLDBOT, IWP, DME, ISCATTAB, . MAXTAB, MAXGRID, NSCATTAB, MUINC, . NMUOBS, MUTAB, NDME, DMETAB, NWAVETAB, WAVETAB, . TABEXTINCT, TABSSALB, TABASYM, . TABPHI1UP, TABPHI1DN, TABPHI2UP, TABPHI2DN, . RADOBS(IBAND,M,J0)) C End of loop over wavenumbers ENDDO C End of loop for sidebands ENDDO C End of loop over spectral segments ENDDO C Average the radiance over desired output bands IF (OUTAVERAGE(1) .EQ. 'R'.OR. OUTAVERAGE(1) .EQ. 'r' .OR. . OUTAVERAGE(1) .EQ. 'T'.OR. OUTAVERAGE(1) .EQ. 't') THEN CALL AVG_SPECTRUM (MAXSPEC, NUMSPEC, NWAVENO, NOUTNU, WAVENO, . OUTNU, RADOBS, OUTRAD, OUTUNITS, OUTAVERAGE, . NBAND, OUTNU1, OUTNU2, BANDCENTER, BANDOFFSET, OUTDELNU) ELSE CALL CONVOLVE_SPECTRUM (MAXSPEC, NWAVENO(1), WAVENO, RADOBS, . OUTUNITS, OUTAVERAGE, XO, ABSDELNU, . OUTNU1, OUTNU2, OUTDELNU, NOUTNU, OUTNU, OUTRAD) ENDIF C Output the spectra, header to file CALL OUTPUT_FILE (OUTPUTFILE, NUMSPEC, TWOSIDEBAND, OUTNU1, . OUTNU2, OUTDELNU, BANDCENTER, BANDOFFSET, OUTUNITS, . OUTAVERAGE, XO, ABSFILE, NSCATTAB, SCATFILE, MUOBS, ZOBS, . SFCTEMP, SFCEMIS1, SFCEMIS2, ZTOP, HEIGHT(ICLDBOT), . NCLDLAY,IWP,DME,ISCATTAB, RTMODEL, NOUTNU, OUTNU,OUTRAD) END SUBROUTINE COMPUTE_RADIATIVE_TRANSFER (RTMODEL, . MUOBS, IOBS, WAVENO, . NLEV, TEMP, TAUGAS, SFCTEMP, SFCEMIS, . NCLDLAY, ICLDTOP, ICLDBOT, IWP, DME, ISCATTAB, . MAXTAB, MAXGRID, NSCATTAB, MUINC, . NMUOBS, MUTAB, NDME, DMETAB, NWAVETAB, WAVETAB, . TABEXTINCT, TABSSALB, TABASYM, . TABPHI1UP, TABPHI1DN, TABPHI2UP, TABPHI2DN, . RADOBS) C Calculates the radiance at the desired level (IOBS) and angle (MUOBS) C for one wavenumber (WAVENO) for the atmosphere and cloud specified. C The temperature and molecular absorption profiles are input C (TEMP, TAUGAS) for the NLEV levels. The multilayer cloud is C between levels ICLDTOP and ICLDBOT. Each of the NCLDLAY cloud layers C has an ice water path (IWP), median particle diameter (DME), and C pointer to the scattering table (ISCATTAB). The scattering properties C are input as tables. Accepts mu > 0 (down looking) or mu < 0 (up looking). C The observed radiance (RADOBS) is output in W/(m^2 ster cm^-1) units. C RTMODEL : S,E,H for single scatter, Eddington, and Hybrid respectively C MUOBS : Cosine of observation angle (negative for looking up) C IOBS : The level number of observation height C WAVENO : Wave number (1/cm) for this monochromatic RT calculation. C NLEV : Number of levels in atmosphere C Note: Levels/layers numbering start at top of atmopshere C TEMP : Temperature array for levels in atmosphere C TAUGAS : Optical depth array for gases in layer of atmosphere C SFCTEMP : Temperature of surface C SFCEMIS : Emissivity of surface C NCLDLAY : Number of cloud layers, ie. # of atm layer occupied by cloud C ICLDTOP : The level number of cloud top C ICLDBOT : The level number of cloud bottom C IWP : Ice Water Path (g/m^2) array for cloud layers C DME : Median Particle Diameter (um) array for cloud layers C ISCATTAB : Scattering table file number array for cloud layers C Note: ISCATAB allows multiple scattering table files C MAXTAB : Maximum number of elements in scattering table C MAXGRID : Maximum number of points on scattring table axis C NSCATTAB : Number of scattering tables/files C MUINC : Cosines of two incident angles of single scattering solution C NMUOBS : Number of viewing cosines in scattering table C MUTAB : Viewing cosines making up grid of scattering table C NDME : Number of Median Particle Diameters in scattering table C DMETAB : Median Particle Diameters making up grid of scattering table C NWAVETAB : Number of wavenumbers in scattering table C WAVETAB : Wavenumbers making up grid of scattering table C TABEXTINCT : Extinction (1/km) matrix, dimensions are wavenumber/Dme C Note : Scattering table calculated for IWC = 1 g/m^2, thus C Extinction = TABEXTINCT*IWP/(layer thickness) C TABSSALB : Single Scattering albedo matrix, dimensions are Wavenumber/Dme C TABASYM : Asymmetry Parameter matrix, dimensions are Wavenumber/Dme C TABPHI : "Phi" function matrix, dimensions are Wavenumber/Dme/Mu C "Phi" represents the amount of incident radiation scattered. C TABPHI1UP : "Phi" up applied to forward scattering radiation, for MU1 C TABPHI2UP : "Phi" up applied to forward scattering radiation, for MU2 C TABPHI1DN : "Phi" down applied to back scattering radiation, for MU1 C TABPHI2DN : "Phi" down applied to back scattering radiation, for MU2 C RADOBS : Output of subroutine, Upwelling radiation at observation height IMPLICIT NONE INTEGER NLEV, IOBS INTEGER NCLDLAY, ICLDTOP, ICLDBOT, ISCATTAB(NCLDLAY) REAL MUOBS, WAVENO, SFCTEMP, SFCEMIS REAL TEMP(NLEV), TAUGAS(NLEV-1) REAL IWP(NCLDLAY), DME(NCLDLAY), IWPTOT REAL RADOBS CHARACTER*1 RTMODEL INTEGER NSCATTAB, MAXTAB, MAXGRID INTEGER NMUOBS(NSCATTAB), NDME(NSCATTAB), NWAVETAB(NSCATTAB) REAL MUTAB(MAXGRID,NSCATTAB) REAL DMETAB(MAXGRID,NSCATTAB), WAVETAB(MAXGRID,NSCATTAB) REAL MUINC(2) REAL TABEXTINCT(MAXTAB,NSCATTAB), TABSSALB(MAXTAB,NSCATTAB) REAL TABASYM(MAXTAB,NSCATTAB) REAL TABPHI1UP(MAXTAB,NSCATTAB), TABPHI1DN(MAXTAB,NSCATTAB) REAL TABPHI2UP(MAXTAB,NSCATTAB), TABPHI2DN(MAXTAB,NSCATTAB) C Local variables INTEGER MAXNZ PARAMETER (MAXNZ=150) INTEGER I, L, N REAL TAUTOT(MAXNZ), COALB(MAXNZ), EXTINCT REAL TAUC(MAXNZ), TAUCG(MAXNZ), SSALB(MAXNZ), ASYM(MAXNZ) REAL PHI1UP(MAXNZ), PHI1DN(MAXNZ) REAL PHI2UP(MAXNZ), PHI2DN(MAXNZ) REAL RADBNDRYUP(2), RADBNDRYDN(2) REAL RADBNDRYUP1(MAXNZ), RADBNDRYDN1(MAXNZ) REAL RADBNDRYUP2(MAXNZ), RADBNDRYDN2(MAXNZ) REAL RAD0UPOBS(MAXNZ), RAD0DNOBS(MAXNZ) REAL FLUXES(3,MAXNZ) REAL FLUXTOP, FLUXBOT, TTOP, TBOT, RADBOTTOM, RADTOP REAL FLUXUP, FLUXDN, FLUXUPSEDD, FLUXDNSEDD C Initialize the optical depths at each level IF (NLEV .GT. MAXNZ) . STOP 'COMPUTE_RADIATIVE_TRANSFER: MAXNZ exceeded' DO N = 1, NLEV TAUTOT(N) = TAUGAS(N) COALB(N) = 1.0 ENDDO IWPTOT = 0 DO N = 1,NCLDLAY IWPTOT = IWPTOT + IWP(N) ENDDO C If clear sky only do GASRT1 and GASRT2 IF (IWPTOT.EQ.0) THEN C Integrate radiance from surface to observer C First go from surface to cloud bottom CALL GASRT1 (ABS(MUOBS), WAVENO, SFCTEMP, SFCEMIS, . NLEV, TEMP, TAUTOT, COALB, ICLDTOP, ICLDBOT, . NCLDLAY, RAD0UPOBS, RAD0DNOBS) IF(MUOBS.GT.0) THEN C if looking down then go from cloud top to observer RADOBS = RAD0UPOBS(1) CALL GASRT2 (WAVENO, NLEV, TEMP, TAUGAS, . ICLDTOP+1, IOBS, MUOBS, RADOBS) ELSEIF(MUOBS.LT.0) THEN C Then go from cloud bottom to observer RADOBS = RAD0DNOBS(NCLDLAY) CALL GASRT2 (WAVENO, NLEV, TEMP, TAUGAS, . ICLDBOT-1, IOBS, ABS(MUOBS), RADOBS) ENDIF C If cloud present do Single-scatter,Eddington or hybrid routine ELSE C Get the optical properties for the cloud layers DO N = ICLDTOP, ICLDBOT - 1 L = N-ICLDTOP+1 C Interpolate to get values of extinction, s.s. albedo, and C phi function values at given obs. mu, waveno, and particle size. C Note: we don't need phi functions for Eddington only solution. I = ISCATTAB(L) IF (RTMODEL .EQ. 'E') THEN CALL INTERP_SCAT_TABLE2 (WAVENO, DME(L), . EXTINCT, SSALB(L), ASYM(L), . NDME(I), DMETAB(1,I), NWAVETAB(I), WAVETAB(1,I), . TABEXTINCT(1,I), TABSSALB(1,I), TABASYM(1,I)) ELSE CALL INTERP_SCAT_TABLE3 (ABS(MUOBS), WAVENO, DME(L), . EXTINCT, SSALB(L), ASYM(L), . PHI1UP(L), PHI1DN(L), PHI2UP(L), PHI2DN(L), . NMUOBS(I), MUTAB(1,I), . NDME(I), DMETAB(1,I), NWAVETAB(I), WAVETAB(1,I), . TABEXTINCT(1,I), TABSSALB(1,I), TABASYM(1,I), . TABPHI1UP(1,I), TABPHI1DN(1,I), . TABPHI2UP(1,I), TABPHI2DN(1,I)) ENDIF C Compute the optical depth of cloud layer, including gas TAUC(L) = IWP(L)*EXTINCT/1000. TAUCG(L) = TAUGAS(N) + TAUC(L) IF (TAUCG(L) .GT. 1.0E-5) THEN SSALB(L) = SSALB(L)*TAUC(L)/TAUCG(L) ELSE SSALB(L) = 0.0 ENDIF TAUTOT(N) = TAUCG(L) COALB(N) = 1.0 - SSALB(L) ENDDO C Calculate the gas/cloud emission only radiative transfer C for incident radiation on each layer CALL GASRT1 (MUINC(1), WAVENO, SFCTEMP, SFCEMIS, . NLEV, TEMP, TAUTOT, COALB, . ICLDTOP, ICLDBOT, NCLDLAY, . RADBNDRYUP1, RADBNDRYDN1) CALL GASRT1 (MUINC(2), WAVENO, SFCTEMP, SFCEMIS, . NLEV, TEMP, TAUTOT, COALB, . ICLDTOP, ICLDBOT, NCLDLAY, . RADBNDRYUP2, RADBNDRYDN2) CALL GASRT1 (ABS(MUOBS), WAVENO, SFCTEMP, SFCEMIS, . NLEV, TEMP, TAUTOT, COALB, . ICLDTOP, ICLDBOT, NCLDLAY, RAD0UPOBS, RAD0DNOBS) C Calculate the fluxes for the multilayer boundary conditions FLUXTOP = MUINC(1)*.5*RADBNDRYDN1(1) . + MUINC(2)*.5*RADBNDRYDN2(1) FLUXBOT = MUINC(1)*.5*RADBNDRYUP1(NCLDLAY) . + MUINC(2)*.5*RADBNDRYUP2(NCLDLAY) C Don't compute multilayer fluxes for a single cloud layer IF (NCLDLAY .LE. 1) THEN FLUXES(2,1) = FLUXTOP FLUXES(1,2) = FLUXBOT ELSE C For multilayer cloud calculate the multilayer Eddington fluxes CALL EDDRTF (NCLDLAY, TAUCG, SSALB, ASYM, . TEMP(ICLDTOP), WAVENO, . FLUXTOP, FLUXBOT, FLUXES) ENDIF C If looking down must loop over cloud layers from bottom to top IF(MUOBS.GT.0) THEN C Introduce scalar variable for radiance at bottom of cloud RADBOTTOM = RAD0UPOBS(NCLDLAY) C Loop over the layers of cloud, from bottom to top DO L = NCLDLAY,1, -1 TTOP = TEMP(L-1+ICLDTOP) TBOT = TEMP(L+ICLDTOP) FLUXUP = FLUXES(1,L+1) FLUXDN = FLUXES(2,L) C Do Single-scattering, Eddington, or Hybrid for each layer IF (RTMODEL .EQ. 'S') THEN RADBNDRYUP(1) = RADBNDRYUP1(L) RADBNDRYUP(2) = RADBNDRYUP2(L) RADBNDRYDN(1) = RADBNDRYDN1(L) RADBNDRYDN(2) = RADBNDRYDN2(L) CALL SSCATRT (MUOBS, TTOP, TBOT, WAVENO, MUINC, . RADBNDRYUP, RADBNDRYDN, RADBOTTOM, . TAUCG(L), SSALB(L), ASYM(L), . PHI1UP(L),PHI1DN(L),PHI2UP(L),PHI2DN(L), . RADOBS) ELSE IF (RTMODEL .EQ. 'E') THEN CALL EDDSCATRT (MUOBS, TTOP, TBOT, WAVENO, . FLUXUP, FLUXDN, RADBOTTOM, . TAUCG(L), SSALB(L), ASYM(L), RADOBS) ELSE FLUXDNSEDD = MUINC(1)*.5*RADBNDRYDN1(L) . + MUINC(2)*.5*RADBNDRYDN2(L) FLUXUPSEDD = MUINC(1)*.5*RADBNDRYUP1(L) . + MUINC(2)*.5*RADBNDRYUP2(L) RADBNDRYUP(1) = RADBNDRYUP1(L) RADBNDRYUP(2) = RADBNDRYUP2(L) RADBNDRYDN(1) = RADBNDRYDN1(L) RADBNDRYDN(2) = RADBNDRYDN2(L) CALL SEDDSCATRT (MUOBS, TTOP, TBOT, WAVENO, MUINC, . RADBNDRYUP, RADBNDRYDN, . FLUXUPSEDD, FLUXDNSEDD, . FLUXUP, FLUXDN, RADBOTTOM, . TAUCG(L), SSALB(L), ASYM(L), . PHI1UP(L),PHI1DN(L),PHI2UP(L),PHI2DN(L), . RADOBS) ENDIF C The calculated upwelling radiance at top of layer is used C as the incident radiance on bottom of next layer. RADBOTTOM = RADOBS ENDDO C Integrate radiance from cloud top to observer CALL GASRT2 (WAVENO, NLEV, TEMP, TAUGAS, . ICLDTOP, IOBS, MUOBS, RADOBS) C If looking up must loop over cloud layers from top to bottom ELSEIF(MUOBS.LT.0) THEN C Introduce scalar variable for radiance at top of cloud RADTOP = RAD0DNOBS(1) C Loop over the layers of cloud, from top to bottom DO L = 1,NCLDLAY TTOP = TEMP(L-1+ICLDTOP) TBOT = TEMP(L+ICLDTOP) FLUXUP = FLUXES(1,L+1) FLUXDN = FLUXES(2,L) C Do Single-scattering, Eddington, or Hybrid for each layer IF (RTMODEL .EQ. 'S') THEN RADBNDRYUP(1) = RADBNDRYUP1(L) RADBNDRYUP(2) = RADBNDRYUP2(L) RADBNDRYDN(1) = RADBNDRYDN1(L) RADBNDRYDN(2) = RADBNDRYDN2(L) CALL SSCATRT (ABS(MUOBS), TBOT, TTOP, WAVENO, MUINC, . RADBNDRYDN, RADBNDRYUP, RADTOP, . TAUCG(L), SSALB(L), ASYM(L), . PHI1UP(L),PHI1DN(L),PHI2UP(L),PHI2DN(L), . RADOBS) ELSE IF (RTMODEL .EQ. 'E') THEN CALL EDDSCATRT (ABS(MUOBS), TBOT, TTOP, WAVENO, . FLUXDN, FLUXUP, RADTOP, . TAUCG(L), SSALB(L), ASYM(L), RADOBS) ELSE FLUXDNSEDD = MUINC(1)*.5*RADBNDRYDN1(L) . + MUINC(2)*.5*RADBNDRYDN2(L) FLUXUPSEDD = MUINC(1)*.5*RADBNDRYUP1(L) . + MUINC(2)*.5*RADBNDRYUP2(L) RADBNDRYUP(1) = RADBNDRYUP1(L) RADBNDRYUP(2) = RADBNDRYUP2(L) RADBNDRYDN(1) = RADBNDRYDN1(L) RADBNDRYDN(2) = RADBNDRYDN2(L) CALL SEDDSCATRT (ABS(MUOBS), TBOT, TTOP, WAVENO, MUINC, . RADBNDRYDN, RADBNDRYUP, . FLUXDNSEDD, FLUXUPSEDD, . FLUXDN, FLUXUP, RADTOP, . TAUCG(L), SSALB(L), ASYM(L), . PHI1UP(L),PHI1DN(L),PHI2UP(L),PHI2DN(L), . RADOBS) ENDIF C The calculated upwelling radiance at bottom of layer is used C as the incident radiance on top of next layer. RADTOP = RADOBS ENDDO C Integrate radiance from cloud top to observer CALL GASRT2 (WAVENO, NLEV, TEMP, TAUGAS, . ICLDBOT, IOBS, ABS(MUOBS), RADOBS) ENDIF ENDIF RETURN END SUBROUTINE SSCATRT (MUOBS, TTOP, TBOT, WAVENO, MUINC, . IBNDRYUP, IBNDRYDN, I0UPOBS, . TAUTOT, SSALB, ASYM, . PHI1UP, PHI1DN, PHI2UP, PHI2DN, IUPOBS) C Calculates upwelling radiance from a single thermally emitting C homogeneous layer using the single scattering radiative transfer C method of Deeter and Evans (1998). C Input Parameters: C MUOBS observation mu (cosine of viewing zenith angle) C TTOP cloud top temperature (K) C TBOT cloud bottom temperature (K) C WAVENO wavenumber (cm^-1) C MUINC cosine zenith angles for the two incident directions C IBNDRYUP two incident radiances upwelling on bottom of layer C IBNDRYDN two incident radiance downwelling on top of layer C I0UPOBS incident radiance at bottom of layer at observation angle C TAUTOT optical depth of layer C SSALB single scattering albedo of layer C ASYM asymmetry parameter of layer C PHI* single scattering functions at the observation angle C C Output Parameters: C IUPOBS outgoing radiance at top at observation angle C C Input and output radiances are in units of W m^-2 sr^-1 cm and C fluxes in units of W m^-2 cm. IMPLICIT NONE REAL MUOBS, TTOP, TBOT, WAVENO, MUINC(2) REAL IBNDRYUP(2), IBNDRYDN(2), I0UPOBS, IUPOBS REAL TAUTOT, SSALB, ASYM REAL PHI1UP, PHI1DN, PHI2UP, PHI2DN INTEGER IMU REAL BETA, COPLNCK, PLANCK0, PLANCK1 REAL EXPBETA, EXPMUOBS, DMUOBS, EXPMAX REAL MU, DMUUP, DMUDN, EXPMU REAL S0UP, S1UP, IUP(2), IDN(2) REAL TERM1, TERM2, TERM3 EXPMAX = 1.E08 IF (TAUTOT .LT. 1.0E-6) THEN IUPOBS = I0UPOBS RETURN ENDIF C Planck function radiances in units of W m^-2 sr^-1 cm PLANCK0 = 1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/TTOP) - 1) PLANCK1 = 1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/TBOT) - 1) BETA = LOG(PLANCK1/PLANCK0)/TAUTOT COPLNCK = (1. - SSALB)*PLANCK0 EXPBETA = PLANCK1/PLANCK0 EXPMUOBS = EXP(-TAUTOT/MUOBS) IF (EXPMUOBS .LT. 1./EXPMAX) EXPMUOBS = 1./EXPMAX DMUOBS = 1. - BETA*MUOBS IF(ABS(DMUOBS).LT.1.E-4) DMUOBS = SIGN(1.E-4,DMUOBS) C Zero order scattering contribution S0UP = (1. - EXPBETA*EXPMUOBS)*COPLNCK/DMUOBS C IMU (1 or 2) is index for discrete mu's in simple model DO IMU = 1, 2 MU = MUINC(IMU) EXPMU = EXP(-TAUTOT/MU) IF (EXPMU .LT. 1./EXPMAX) EXPMU = 1./EXPMAX DMUUP = 1. - BETA*MU DMUDN = 1. + BETA*MU IF(ABS(DMUUP).LT.1.E-4) DMUUP = SIGN(1.E-4,DMUUP) IF(ABS(DMUDN).LT.1.E-4) DMUDN = SIGN(1.E-4,DMUDN) TERM1 = IBNDRYUP(IMU)*EXPMU - COPLNCK*EXPBETA*EXPMU/DMUUP TERM2 = (MU/(MUOBS - MU))*(1. - EXPMUOBS/EXPMU) TERM3 = COPLNCK*(1. - EXPBETA*EXPMUOBS)/(DMUUP*DMUOBS) IUP(IMU) = - (TERM1*TERM2 - TERM3) TERM1 = IBNDRYDN(IMU) - COPLNCK/DMUDN TERM2 = (MU/(MUOBS + MU))*(1. - EXPMUOBS*EXPMU) TERM3 = COPLNCK*(1. - EXPBETA*EXPMUOBS)/(DMUDN*DMUOBS) IDN(IMU) = (TERM1*TERM2 + TERM3) ENDDO C First order scattering contribution S1UP = (SSALB/2.)*(PHI1UP*IUP(1) + PHI2UP*IUP(2) . + PHI1DN*IDN(1) + PHI2DN*IDN(2)) IUPOBS = I0UPOBS*EXPMUOBS + S0UP + S1UP RETURN END SUBROUTINE EDDSCATRT (MUOBS, TTOP, TBOT, WAVENO, . FLUXUP, FLUXDN, I0UPOBS, . TAUTOT, SSALB, ASYM, IUPOBS) C Calculates upwelling radiance from a single thermally emitting C homogeneous layer using Eddington's second approximation C radiative transfer method. C C Input Parameters: C MUOBS observation mu (cosine of viewing zenith angle) C TTOP cloud top temperature (K) C TBOT cloud bottom temperature (K) C WAVENO wavenumber (cm^-1) C FLUXUP incident flux on bottom from multilayer Eddington solution C FLUXDN incident flux on top from multilayer Eddington solution C I0UPOBS incident radiance at bottom of layer at observation angle C TAUTOT optical depth of layer C SSALB single scattering albedo of layer C ASYM asymmetry parameter of layer C C Output Parameters: C IUPOBS outgoing radiance at top at observation angle C C Input and output radiances are in units of W m^-2 sr^-1 cm and C fluxes in units of W m^-2 cm. IMPLICIT NONE REAL MUOBS, TTOP, TBOT REAL WAVENO, FLUXUP,FLUXDN REAL I0UPOBS REAL IUPOBS REAL SSALB, ASYM, TAUTOT REAL BETA, PLANCK0, PLANCK1, COPLANCK REAL EXPBETA, EXPMUOBS, EXPLAMP, EXPLAMM, EXPMAX REAL T, R, LAMBDA REAL F1UP, F0DN, FP1UP, FP0DN, FH1UP, FH0DN REAL PARTFAC, PARTP, PARTM REAL CFAC, CP, CM, DB, D1, D2, DP, DM, SB, SP, SM EXPMAX = 1.E08 IF (TAUTOT .LT. 1.0E-6) THEN IUPOBS = I0UPOBS RETURN ENDIF C Planck function radiances in units of W m^-2 sr^-1 cm PLANCK0 = 1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/TTOP) - 1) PLANCK1 = 1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/TBOT) - 1) C Compute the Planck function variation coefficient BETA = LOG(PLANCK1/PLANCK0)/TAUTOT EXPBETA = PLANCK1/PLANCK0 C Find lambda so that it can be compared to beta C R, T are 2x2 matrix coefficents, lambda is eigenvalue T = (7-SSALB*(4+3*ASYM))*0.25 R = (1-SSALB*(4-3*ASYM))*0.25 C solution explodes if r = 0!! IF (ABS(R) .LT. 1.E-5) R = 1.E-5 LAMBDA = SQRT(T**2-R**2) IF(ABS(BETA**2-LAMBDA**2).LT.1.E-4) THEN BETA = BETA/ABS(BETA)*SQRT(LAMBDA**2+1.E-4) EXPBETA = EXP(BETA*TAUTOT) ENDIF COPLANCK = (1-SSALB)*PLANCK0 EXPMUOBS = EXP(-TAUTOT/MUOBS) IF (EXPMUOBS .LT. 1./EXPMAX) EXPMUOBS = 1./EXPMAX EXPLAMP = EXP(LAMBDA*TAUTOT) IF (EXPLAMP .GT. EXPMAX) EXPLAMP = EXPMAX EXPLAMM = 1.0/EXPLAMP F1UP = FLUXUP F0DN = FLUXDN C Compute particular solution at boundaries PARTFAC = COPLANCK/(LAMBDA**2-BETA**2) PARTP = PARTFAC*(T-R+BETA) PARTM = PARTFAC*(T-R-BETA) FP0DN = PARTM FP1UP = PARTP*EXPBETA C Subtract particular solution off to get homogeneous solution FH0DN = F0DN - FP0DN FH1UP = F1UP - FP1UP C Compute homogeneous coefficients C+ and C- CFAC = 1.0/(R**2*EXPLAMP - (LAMBDA-T)**2*EXPLAMM) CP = CFAC*(R*FH1UP - (LAMBDA-T)*EXPLAMM*FH0DN) CM = CFAC*(R*EXPLAMP*FH0DN - (LAMBDA-T)*FH1UP) C Make the constants for the 3 source function terms DB = COPLANCK + SSALB*PARTFAC*(2*(T-R) + 3*ASYM*BETA*MUOBS) D1 = R-T+LAMBDA D2 = 1.5*ASYM*(R+T-LAMBDA)*MUOBS DP = CP*SSALB*(D1+D2) DM = CM*SSALB*(D1-D2) C Finally do the answer IF(ABS(1-BETA*MUOBS).LE.1.E-5) THEN SB = DB*TAUTOT/MUOBS ELSE SB = DB*(1 - EXPBETA*EXPMUOBS)/(1-BETA*MUOBS) ENDIF IF(ABS(1-LAMBDA*MUOBS).LE.1.E-5) THEN SP = DP*TAUTOT/MUOBS ELSE SP = DP*(1 - EXPLAMP*EXPMUOBS)/(1-LAMBDA*MUOBS) ENDIF SM = DM*(1 - EXPLAMM*EXPMUOBS)/(1+LAMBDA*MUOBS) IUPOBS = I0UPOBS*EXPMUOBS + SB + SP + SM RETURN END SUBROUTINE SEDDSCATRT (MUOBS, TTOP, TBOT, WAVENO, . MUINC, IBNDRYUP, IBNDRYDN, . FLUXUPSEDD, FLUXDNSEDD, . FLUXUP, FLUXDN, I0UPOBS, . TAUTOT, SSALB, ASYM, . PHI1UP, PHI1DN, PHI2UP, PHI2DN, IUPOBS) C Calculates upwelling radiance from a single thermally emitting C homogeneous layer using the single scattering/Eddington second C approximation hybrid radiative transfer method of Deeter and Evans (1998). C C Input Parameters: C MUOBS observation mu (cosine of viewing zenith angle) C TTOP cloud top temperature (K) C TBOT cloud bottom temperature (K) C WAVENO wavenumber (cm^-1) C MUINC cosine zenith angles for the two incident directions C IBNDRYUP two incident radiances upwelling on bottom of layer C IBNDRYDN two incident radiance downwelling on top of layer C FLUXUPSEDD incident flux on bottom from single scattering/Eddington C FLUXDNSEDD incident flux on top from single scattering/Eddington C FLUXUP incident flux on bottom from multilayer Eddington solution C FLUXDN incident flux on top from multilayer Eddington solution C I0UPOBS incident radiance at bottom of layer at observation angle C TAUTOT optical depth of layer C SSALB single scattering albedo of layer C ASYM asymmetry parameter of layer C PHI* single scattering functions at the observation angle C C Output Parameters: C IUPOBS outgoing radiance at top at observation angle C C Input and output radiances are in units of W m^-2 sr^-1 cm and C fluxes in units of W m^-2 cm. IMPLICIT NONE REAL MUOBS, TTOP, TBOT REAL WAVENO, TAUTOT REAL IBNDRYUP(2), IBNDRYDN(2) REAL I0UPOBS REAL IUPOBS, FLUXUP, FLUXDN, FLUXUPSEDD, FLUXDNSEDD REAL MUINC(2) INTEGER IMU REAL SSALB, ASYM, PHI1UP, PHI1DN, PHI2UP, PHI2DN REAL BETA, COPLANCK, PLANCK0, PLANCK1 REAL EXPBETA, EXPMUOBS, DMUOBS, EXPMAX REAL MU, DMUUP, DMUDN, EXPMU REAL S1UP, IUP(2), IDN(2) REAL TERM1, TERM2, TERM3 REAL EXPLAMP, EXPLAMM REAL T, R, LAMBDA REAL F1UP, F0DN, FP1UP, FP0DN, FH1UP, FH0DN REAL PARTFAC, PARTP, PARTM REAL CFAC, CP, CM, DB, D1, D2, DP, DM, SB, SP, SM REAL SQR3, EXP3TP, EXP3TM, FAC3TP, FAC3TM, FAC3BETA REAL FACF0DN, FACF1UP, FACG, FACGP, FACGM REAL FACDENOM, DENOM1, DENOM2, FACNUM1, FACNUM2 REAL SEDTERM1, SEDTERM2, FACMUBET, FACNUM3, SEDTERM3, S1EDD IF (TAUTOT .LT. 1.0E-6) THEN IUPOBS = I0UPOBS RETURN ENDIF EXPMAX = 1.E08 SQR3 = SQRT(3.0) C Planck function radiances in units of W m^-2 sr^-1 cm PLANCK0 = 1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/TTOP) - 1) PLANCK1 = 1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/TBOT) - 1) BETA = LOG(PLANCK1/PLANCK0)/TAUTOT EXPBETA = PLANCK1/PLANCK0 C Find lambda so that it can be compared to beta C R, T are 2x2 matrix coefficents, lambda is eigenvalue T = (7-SSALB*(4+3*ASYM))*0.25 R = (1-SSALB*(4-3*ASYM))*0.25 C solution explodes if r = 0!! IF (ABS(R) .LT. 1.E-5) R = 1.E-5 LAMBDA = SQRT(T**2-R**2) IF(ABS(BETA**2-LAMBDA**2).LT.1.E-4) THEN BETA = BETA/ABS(BETA)*SQRT(LAMBDA**2+1.E-4) EXPBETA = EXP(BETA*TAUTOT) ENDIF COPLANCK = (1. - SSALB)*PLANCK0 EXPMUOBS = EXP(-TAUTOT/MUOBS) IF (EXPMUOBS .LT. 1./EXPMAX) EXPMUOBS = 1./EXPMAX DMUOBS = 1. - BETA*MUOBS IF(ABS(DMUOBS).LT.1.E-4) DMUOBS = SIGN(1.E-4,DMUOBS) C Accurate Single Scattering Part C IMU (1 or 2) is index for discrete mu's in simple model DO IMU = 1, 2 MU = MUINC(IMU) EXPMU = EXP(-TAUTOT/MU) IF (EXPMU .LT. 1./EXPMAX) EXPMU = 1./EXPMAX DMUUP = 1. - BETA*MU DMUDN = 1. + BETA*MU IF(ABS(DMUUP).LT.1.E-4) DMUUP = SIGN(1.E-4,DMUUP) IF(ABS(DMUDN).LT.1.E-4) DMUDN = SIGN(1.E-4,DMUDN) TERM1 = IBNDRYUP(IMU)*EXPMU - COPLANCK*EXPBETA*EXPMU/DMUUP TERM2 = (MU/(MUOBS - MU))*(1. - EXPMUOBS/EXPMU) TERM3 = COPLANCK*(1. - EXPBETA*EXPMUOBS)/(DMUUP*DMUOBS) IUP(IMU) = - (TERM1*TERM2 - TERM3) TERM1 = IBNDRYDN(IMU) - COPLANCK/DMUDN TERM2 = (MU/(MUOBS + MU))*(1. - EXPMUOBS*EXPMU) TERM3 = COPLANCK*(1. - EXPBETA*EXPMUOBS)/(DMUDN*DMUOBS) IDN(IMU) = (TERM1*TERM2 + TERM3) ENDDO C First order scattering contribution S1UP = (SSALB/2.)*(PHI1UP*IUP(1) + PHI2UP*IUP(2) . + PHI1DN*IDN(1) + PHI2DN*IDN(2)) C Eddington part EXPLAMP = EXP(LAMBDA*TAUTOT) IF (EXPLAMP .GT. EXPMAX) EXPLAMP = EXPMAX EXPLAMM = 1.0/EXPLAMP C Fluxes sent into subroutine from output of multilayer solution F1UP = FLUXUP F0DN = FLUXDN C Compute particular solution at boundaries PARTFAC = COPLANCK/(LAMBDA**2-BETA**2) PARTP = PARTFAC*(T-R+BETA) PARTM = PARTFAC*(T-R-BETA) FP0DN = PARTM FP1UP = PARTP*EXPBETA C Subtract particular solution off to get homogeneous solution FH0DN = F0DN - FP0DN FH1UP = F1UP - FP1UP C Compute homogeneous coefficients C+ and C- CFAC = 1.0/(R**2*EXPLAMP - (LAMBDA-T)**2*EXPLAMM) CP = CFAC*(R*FH1UP - (LAMBDA-T)*EXPLAMM*FH0DN) CM = CFAC*(R*EXPLAMP*FH0DN - (LAMBDA-T)*FH1UP) C Make the constants for the 3 source function terms DB = COPLANCK + SSALB*PARTFAC*(2*(T-R) + 3*ASYM*BETA*MUOBS) D1 = R-T+LAMBDA D2 = 1.5*ASYM*(R+T-LAMBDA)*MUOBS DP = CP*SSALB*(D1+D2) DM = CM*SSALB*(D1-D2) C Finally do the answer IF(ABS(1-BETA*MUOBS).LE.1.E-5) THEN SB = DB*TAUTOT/MUOBS ELSE SB = DB*(1 - EXPBETA*EXPMUOBS)/(1-BETA*MUOBS) ENDIF IF(ABS(1-LAMBDA*MUOBS).LE.1.E-5) THEN SP = DP*TAUTOT/MUOBS ELSE SP = DP*(1 - EXPLAMP*EXPMUOBS)/(1-LAMBDA*MUOBS) ENDIF SM = DM*(1 - EXPLAMM*EXPMUOBS)/(1+LAMBDA*MUOBS) C Compute the Eddington single scattering solution F0DN = FLUXDNSEDD F1UP = FLUXUPSEDD EXP3TP = EXP(SQR3*TAUTOT) IF (EXP3TP .GT. EXPMAX) EXP3TP = EXPMAX EXP3TM = 1./EXP3TP FAC3TP = 1 - EXPMUOBS*EXP3TP FAC3TM = 1 - EXPMUOBS*EXP3TM FAC3BETA = 1./(3. - BETA**2) FACF0DN = F0DN - (1.5 - BETA)*COPLANCK*FAC3BETA FACF1UP = F1UP - (1.5 + BETA)*COPLANCK*EXPBETA*FAC3BETA FACG = 1.5*(2. - SQR3)*ASYM*MUOBS FACGP = -1.5 + SQR3 + FACG FACGM = -1.5 + SQR3 - FACG FACDENOM = EXP3TM*(-7./4. + SQR3)**2 - EXP3TP/16. DENOM1 = FACDENOM*(1 + MUOBS*SQR3) DENOM2 = -FACDENOM*(1 - MUOBS*SQR3) FACNUM1 = -.25*EXP3TP*FACF0DN + (-7./4.+ SQR3)*FACF1UP FACNUM2 = -(-7./4. + SQR3)*EXP3TM*FACF0DN + .25*FACF1UP SEDTERM1 = FAC3TM*FACNUM1*FACGM/DENOM1 SEDTERM2 = FAC3TP*FACNUM2*FACGP/DENOM2 FACMUBET = 1. - EXPBETA*EXPMUOBS FACNUM3 = COPLANCK*3.*(1 + BETA*ASYM*MUOBS)*FAC3BETA SEDTERM3 = FACMUBET*FACNUM3/DMUOBS S1EDD = SSALB*(SEDTERM1 + SEDTERM2 + SEDTERM3) C Zero order scattering contribution C S0UP = (1. - EXPBETA*EXPMUOBS)*COPLANCK/DMUOBS IUPOBS = I0UPOBS*EXPMUOBS + SB + SP + SM - S1EDD + S1UP RETURN END SUBROUTINE EDDRTF (NLAYER, OPTDEPTHS, ALBEDOS, . ASYMMETRIES, TEMPS, WAVENO, . FLUXTOP, FLUXBOT, FLUXES) C EDDRTF computes the layer interface fluxes for a multilayer C plane-parallel atmosphere with thermal sources of radiation using the C Eddington approximation. The medium is specified by a number C of homogeneous layers. The Planck function is assumed linear with C optical depth (slightly inconsistent with the exponential assumption C in the single layer routines). The temperatures, optical depth, C single scattering albedo, and asymmetry parameter are specified C for each layer. The boundary conditions are the fluxes incident C on the top and bottom of the domain. The Eddington fluxes at each C level are returned. C The model works by calculating the reflection, transmission, and C source terms for each layer from the input properties. A C tri-diagonal matrix solver is then used to compute the fluxes C at each layer interface from the applied boundary conditions. C C Parameters: C Input: C NLAYER integer Number of homogenous layers C (layers are specified from the top down) C OPTDEPTHS real array Optical thickness of layers C ALBEDOS real array Single scattering albedos C ASYMMETRIES real array Asymmetry parameters C TEMPS real array Temperatures (K) at layer interfaces C (e.g. TEMPS(1) is at top of top layer, C TEMPS(2) is at bottom of top layer). C WAVENO real wavenumber (cm^-1, for Planck function) C C Output: C FLUXES real Eddington fluxes at layer interfaces. C FLUXES(1,L) is upwelling, C FLUXES(2,L) is downwelling, C L=1 is top, L=NUML+1 is bottom IMPLICIT NONE INTEGER NLAYER REAL TEMPS(NLAYER+1) REAL OPTDEPTHS(NLAYER) REAL ALBEDOS(NLAYER) REAL ASYMMETRIES(NLAYER) REAL WAVENO REAL FLUXES(3,(NLAYER+1)) C MAXLAY is maximum number of layers INTEGER MAXLAY, MAXN PARAMETER (MAXLAY=200, MAXN=2*MAXLAY+1) INTEGER N, L, I REAL*8 DELTAU, G, OMEGA REAL*8 LAMBDA, R, T, D, CM, CP, A, B, X1, X2 REAL*8 REFLECT, TRANS, SOURCEP, SOURCEM REAL*8 RADP1P, RADP1M, RADP2P, RADP2M REAL*8 PI, PLANCK1, PLANCK2, TAU REAL*8 EXLP, EXLM, V REAL*8 LOWER(MAXN), UPPER(MAXN), DIAG(MAXN), RHS(MAXN) REAL FLUXTOP, FLUXBOT PARAMETER (PI=3.1415926535) C Compute the reflection, transmission, and source C coefficients for each layer for the diffuse Eddington C two stream problem. N = 2*NLAYER+2 IF (N .GT. MAXN) STOP 'EDDRTF: Exceeded maximum number of layers' PLANCK1 = 0.5*1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/TEMPS(1)) - 1) TAU = 0.0 I = 2 DO L = 1, NLAYER DELTAU = OPTDEPTHS(L) IF (DELTAU .LT. 0.0) STOP 'EDDRTF: TAU<0' C Special case for zero optical depth IF (DELTAU .EQ. 0.0) THEN TRANS = 1.0 REFLECT = 0.0 SOURCEP = 0.0 SOURCEM = 0.0 ELSE OMEGA = ALBEDOS(L) G = ASYMMETRIES(L) R = ( 1.0 - OMEGA*(4.0-3.0*G) )/4.0 T = ( 7.0 - OMEGA*(4.0+3.0*G) )/4.0 LAMBDA = SQRT( 3.0*(1.0-OMEGA)*(1.0-OMEGA*G) ) C Special case for conservative scattering (lambda=0) IF (LAMBDA .EQ. 0.0) THEN D = 1.0/(1.0+T*DELTAU) TRANS = D REFLECT = -R*DELTAU*D ELSE X1 = -R X2 = LAMBDA + T EXLP = DEXP(MIN(LAMBDA*DELTAU,75.D0)) EXLM = 1.0/EXLP TRANS = 2.*LAMBDA/(X2*EXLP + (LAMBDA-T)*EXLM) REFLECT = X1*(EXLP - EXLM) *TRANS /(2.*LAMBDA) D = 1.0/(X2**2 *EXLP - X1**2 *EXLM) ENDIF C Calculate thermal source terms PLANCK2 = 0.5*1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/TEMPS(L+1)) - 1) V = 2.0*(PLANCK2-PLANCK1)/(3.0*(1.-OMEGA*G)*DELTAU) RADP1P = -V + PLANCK1 RADP2M = V + PLANCK2 RADP2P = -V + PLANCK2 RADP1M = V + PLANCK1 IF (LAMBDA .EQ. 0.0) THEN SOURCEP = 0 SOURCEM = 0 ELSE CP = (X1*EXLM*RADP1P - X2*RADP2M) *D CM = (-X2*EXLP*RADP1P + X1*RADP2M) *D SOURCEP = X1*CP*EXLP + X2*CM*EXLM + RADP2P SOURCEM = X2*CP + X1*CM + RADP1M ENDIF PLANCK1 = PLANCK2 FLUXES(3,L) = 0.0 ENDIF DIAG(I) = -REFLECT DIAG(I+1) = -REFLECT LOWER(I) = 1.0 LOWER(I+1) = -TRANS UPPER(I) = -TRANS UPPER(I+1) = 1.0 RHS(I) = SOURCEM RHS(I+1) = SOURCEP I = I + 2 ENDDO C Setup for and call the tri-diagonal matrix solver RHS(1) = FLUXTOP DIAG(1) = 0.0 UPPER(1) = 1.0 c DIAG(N) = -(1.0-GNDEMIS) DIAG(N) = 0.0 LOWER(N) = 1.0 RHS(N) = FLUXBOT CALL TRIDIAG (N, LOWER, DIAG, UPPER, RHS) C Put the fluxes in the output array I = 1 DO L = 1, NLAYER+1 FLUXES(1,L) = RHS(I) FLUXES(2,L) = RHS(I+1) I = I + 2 ENDDO RETURN END SUBROUTINE TRIDIAG (N, LOWER, DIAG, UPPER, RHS) C Computes the solution to a tridiagonal system. C N is order of the matrix. LOWER(2..N) is the subdiagonal, C DIAG(1..N) is the diagonal, and UPPER(1..N-1) is the C superdiagonal. On input RHS is the right hand side, while C on output it is the solution vector. Everything is destroyed. C Hacked from Linpack DGTSL. IMPLICIT NONE INTEGER N REAL*8 LOWER(*), DIAG(*), UPPER(*), RHS(*) INTEGER K, KB REAL*8 T IF (N .EQ. 1) THEN IF (DIAG(1) .EQ. 0.0) GOTO 990 RHS(1) = RHS(1)/DIAG(1) ENDIF LOWER(1) = DIAG(1) DIAG(1) = UPPER(1) UPPER(1) = 0.0 UPPER(N) = 0.0 DO K = 1, N-1 C Interchange this and next row to the get the largest pivot. IF (ABS(LOWER(K+1)) .GE. ABS(LOWER(K))) THEN T = LOWER(K+1) LOWER(K+1) = LOWER(K) LOWER(K) = T T = DIAG(K+1) DIAG(K+1) = DIAG(K) DIAG(K) = T T = UPPER(K+1) UPPER(K+1) = UPPER(K) UPPER(K) = T T = RHS(K+1) RHS(K+1) = RHS(K) RHS(K) = T ENDIF IF (LOWER(K) .EQ. 0.0) GOTO 990 T = -LOWER(K+1)/LOWER(K) LOWER(K+1) = DIAG(K+1) + T*DIAG(K) DIAG(K+1) = UPPER(K+1) + T*UPPER(K) UPPER(K+1) = 0.0 RHS(K+1) = RHS(K+1) + T*RHS(K) ENDDO IF (LOWER(N) .EQ. 0.0) GOTO 990 C Back substitute RHS(N) = RHS(N)/LOWER(N) RHS(N-1) = (RHS(N-1) - DIAG(N-1)*RHS(N))/LOWER(N-1) DO KB = 1, N-2 K = N - 2 - KB + 1 RHS(K) = (RHS(K) -DIAG(K)*RHS(K+1) -UPPER(K)*RHS(K+2))/LOWER(K) ENDDO RETURN 990 CONTINUE STOP 'Singular matrix in TRIDIAG' END SUBROUTINE GASRT1 (MU, WAVENO, SFCTEMP, SFCEMIS, . NLEV, TEMP, TAU, COALB, ICLDTOP, . ICLDBOT, NCLDLAY, RAD1UP, RAD0DN) C Compute the upwelling radiance at the bottom of each cloud layer C (RAD1UP) and the downwelling radiance at the top of each cloud C layer (RAD0DN). Integrates the nonscattering RTE through the domain C from top to bottom and then bottom to top. Only the absorption C part of the optical depth in the cloud is considered. The bottom C of the lowest cloud layer is at level ICLDBOT and the top of the C highest cloud layer is as level ICLDTOP. The radiance is computed C at the angle MU and wavenumber WAVENO. The temperature profile (TEMP), C layer optical depth profile (TAU), and single scattering coalbedo C (COALB=1-omega) profile are input. There is specular reflection from C the surface yet no multiple scattering effect due of the surface, ie. C The gas integrated down from TOA and reflected off surface once with C reflectivity equal to one minus emissivity, and this is the initial C condition when integrating gas incident on bottom of cloud. IMPLICIT NONE INTEGER NLEV, ICLDTOP, ICLDBOT, NCLDLAY REAL MU, WAVENO, SFCTEMP, SFCEMIS REAL TEMP(NLEV), TAU(NLEV), COALB(NLEV) REAL RAD1UP(NCLDLAY), RAD0DN(NCLDLAY) REAL RAD1UP0, RAD0DN0 INTEGER I, L, IEND REAL PLANCK0, PLANCK1, PLANCKSFC, TAUO, TRANS, DELPLANCK C Loop through layers, starting at top and going down C Continue past cloud layer to surface to compute specular C reflection. RAD0DN0 = 1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/2.7) - 1) c RAD0DN0 = 0.0 IF(ICLDBOT.EQ.2) RAD0DN(1) = RAD0DN0 PLANCK0 = 1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/TEMP(1)) - 1) L = 1 IF(SFCEMIS.LT.1) IEND = NLEV-1 IF(SFCEMIS.EQ.1) IEND = ICLDBOT-2 DO I = 1, IEND C Planck function radiances in units of W m^-2 sr^-1 cm PLANCK1 = 1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/TEMP(I+1)) - 1) TAUO = TAU(I)/MU IF (TAU(I) .LT. 0.001) THEN RAD0DN0 = RAD0DN0*(1-TAUO) . + TAUO*0.5*(PLANCK0+PLANCK1)*COALB(I) ELSE TRANS = EXP(-TAUO) DELPLANCK = (PLANCK1-PLANCK0)/TAUO RAD0DN0 = RAD0DN0*TRANS + COALB(I)*( PLANCK1-DELPLANCK . - TRANS*(PLANCK1-DELPLANCK*(1.0+TAUO)) ) ENDIF IF (I+1 .GE. ICLDTOP . AND. I .LE. (ICLDBOT-2)) THEN RAD0DN(L) = RAD0DN0 L = L + 1 ENDIF PLANCK0 = PLANCK1 ENDDO C Loop through layers, starting at bottom and going up L = NCLDLAY PLANCKSFC = 1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/SFCTEMP) - 1) RAD1UP0 = SFCEMIS*PLANCKSFC+RAD0DN0*(1-SFCEMIS) IF(ICLDTOP.EQ.(NLEV-1)) RAD1UP(1) = RAD1UP0 PLANCK1 = 1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/TEMP(NLEV)) - 1) DO I = NLEV-1, ICLDTOP+1, -1 PLANCK0 = 1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/TEMP(I)) - 1) TAUO = TAU(I)/MU IF (TAUO .LT. 0.001) THEN RAD1UP0 = RAD1UP0*(1-TAUO) . + TAUO*0.5*(PLANCK0+PLANCK1)*COALB(I) ELSE TRANS = EXP(-TAUO) DELPLANCK = (PLANCK1-PLANCK0)/TAUO RAD1UP0 = RAD1UP0*TRANS + COALB(I)*( PLANCK0+DELPLANCK . - TRANS*(PLANCK0+DELPLANCK*(1.0+TAUO)) ) ENDIF IF (I .LE. ICLDBOT) THEN RAD1UP(L) = RAD1UP0 L = L - 1 ENDIF PLANCK1 = PLANCK0 ENDDO RETURN END SUBROUTINE GASRT2 (WAVENO, NLEV, TEMP, TAU, ICLD, IOBS, MU, RAD) C Compute the upwelling radiance starting at the cloud top (ICLD) C and going to the observation level (IOBS) or the downwelling C radiance starting at cloud bottom (ICLD) and going down to C observation level (IOBS) depending on MU. C The radiance is computed at the angle MU and wavenumber WAVENO. C The temperature profile (TEMP) and layer optical depth profile (TAU) C are input. IMPLICIT NONE INTEGER NLEV, ICLD, IOBS REAL WAVENO, TEMP(NLEV), TAU(NLEV), MU, RAD INTEGER I REAL PLANCK0, PLANCK1, TAUO, TRANS, DELPLANCK C Loop through layers, if looking down start at bottom and go up IF(IOBS.LT.ICLD) THEN PLANCK1 = 1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/TEMP(ICLD)) - 1) DO I = ICLD-1, IOBS, -1 PLANCK0 = 1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/TEMP(I)) - 1) TAUO = TAU(I)/MU IF (TAUO .LT. 0.001) THEN RAD = RAD*(1-TAUO) + TAUO*0.5*(PLANCK0+PLANCK1) ELSE TRANS = EXP(-TAUO) DELPLANCK = (PLANCK1-PLANCK0)/TAUO RAD = RAD*TRANS + PLANCK0+DELPLANCK . - TRANS*(PLANCK0+DELPLANCK*(1.0+TAUO)) ENDIF PLANCK1 = PLANCK0 ENDDO C Loop through layers, if looking up start at top and go down ELSEIF(IOBS.GT.ICLD) THEN PLANCK0 = 1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/TEMP(ICLD)) - 1) DO I = ICLD, IOBS-1 C Planck function radiances in units of W m^-2 sr^-1 um^-1 PLANCK1 = 1.1911E-8 *WAVENO**3 . / (EXP(1.4388*WAVENO/TEMP(I+1)) - 1) TAUO = TAU(I)/MU IF (TAU(I) .LT. 0.001) THEN RAD = RAD*(1-TAUO) . + TAUO*0.5*(PLANCK0+PLANCK1) ELSE TRANS = EXP(-TAUO) DELPLANCK = (PLANCK1-PLANCK0)/TAUO RAD = RAD*TRANS + ( PLANCK1-DELPLANCK . - TRANS*(PLANCK1-DELPLANCK*(1.0+TAUO)) ) ENDIF PLANCK0 = PLANCK1 ENDDO ENDIF RETURN END SUBROUTINE READ_SSCATTAB (SCATFILE, MAXTAB, MAXGRID, . NMUOBS, MUTAB, NDME, DMETAB, NWAVE, WAVETAB, . MUINC, TABEXTINCT, TABSSALB, TABASYM, . TABPHI1UP, TABPHI1DN, TABPHI2UP, TABPHI2DN) C Reads in the single scattering table for a number of wavenumbers, C particle sizes, and viewing angles. The scattering properties are C computed for a IWC/LWC of 1 g/m^3. C Input parameters: C SCATFILE file name of scattering file C MAXTAB maximum array size for whole table C MAXGRID maximum array size for table grids C Output parameters: C NMUOBS number of viewing angle mu grid values C MUTAB viewing angle grid values C NDME number of Dme grid values C DMETAB Dme grid values C NWAVE number of wavenumber grid values C WAVETAB wavenumber grid values C MUINC(2) cosine zenith of two incident angles C TABEXTINCT tabulated extinction (km^-1) C TABSSALB tabulated single scattering albedo C TABASYM tabulated asymmetry parameter C TABPHI1UP, TABPHI1DN, TABPHI2UP, TABPHI2DN C tabulated phase function terms for incident radiance angles IMPLICIT NONE INTEGER MAXTAB, MAXGRID INTEGER NMUOBS, NDME, NWAVE REAL MUTAB(*), DMETAB(*), WAVETAB(*) REAL MUINC(2) REAL TABEXTINCT(*), TABSSALB(*), TABASYM(*) REAL TABPHI1UP(*), TABPHI1DN(*) REAL TABPHI2UP(*), TABPHI2DN(*) CHARACTER*(*) SCATFILE INTEGER IMU, ID, IW, K2, K3 OPEN (UNIT=2, STATUS='OLD', FILE=SCATFILE) READ (2,*) READ (2,*) NMUOBS READ (2,*) NDME READ (2,*) NWAVE IF (MAX(NMUOBS,NDME,NWAVE) .GT. MAXGRID) . STOP 'READ_SSCATTAB: MAXGRID exceeded' IF (NMUOBS*NDME*NWAVE .GT. MAXTAB) . STOP 'READ_SSCATTAB: MAXTAB exceeded' READ (2,*) MUINC(1), MUINC(2) READ (2,*) READ (2,*) READ (2,*) READ (2,*) READ (2,*) READ (2,*) (MUTAB(IMU), IMU = 1, NMUOBS) DO IW = 1, NWAVE DO ID = 1, NDME K2 = IW-1 + NWAVE*(ID-1) K3 = NMUOBS*K2 READ(2,*) READ(2,*) READ(2,*) DMETAB(ID), WAVETAB(IW), TABEXTINCT(K2+1), . TABSSALB(K2+1), TABASYM(K2+1) READ(2,*) (TABPHI1UP(IMU+K3), IMU = 1, NMUOBS) READ(2,*) (TABPHI2UP(IMU+K3), IMU = 1, NMUOBS) READ(2,*) (TABPHI1DN(IMU+K3), IMU = 1, NMUOBS) READ(2,*) (TABPHI2DN(IMU+K3), IMU = 1, NMUOBS) ENDDO ENDDO CLOSE (UNIT=2) RETURN END SUBROUTINE INTERP_SCAT_TABLE2 (WAVENO, DME, . EXTINCT, SSALB, ASYM, . NDME, DMETAB, NWAVE, WAVETAB, . TABEXTINCT, TABSSALB, TABASYM) C Interpolates the scattering properties from the table for C a particular wavenumber and particle size. Does a bilinear C interpolation, but optimized for fixed particle size and slowly C varying wavenumber. If the DME is the same as last time then we C can just linearly interpolate in wavenumber between stored C scattering values. If the DME has changed then we linearly C interpolate between the DMETAB grid lines for each of the two C wavenumber grid lines. IMPLICIT NONE REAL WAVENO, DME REAL EXTINCT, SSALB, ASYM INTEGER NDME, NWAVE REAL DMETAB(NDME), WAVETAB(NWAVE) REAL TABEXTINCT(NWAVE,NDME), TABSSALB(NWAVE,NDME) REAL TABASYM(NWAVE,NDME) INTEGER IW0, IW1, ID, IL, IU, IM LOGICAL NEWIW REAL FWAV, FDME, FLDME, F REAL OLDDME, EXT0, EXT1, ALB0, ALB1, ASYM0, ASYM1 SAVE IW0, IW1, ID, OLDDME, FDME, FLDME SAVE EXT0,EXT1, ALB0,ALB1, ASYM0,ASYM1 DATA IW0/1/, IW1/2/ C Check that parameter are in range of table IF (WAVENO .LT. WAVETAB(1) .OR. WAVENO .GT. WAVETAB(NWAVE)) . STOP 'INTERP_SCAT_TABLE: wavenumber out of range.' IF (DME .LT. DMETAB(1) .OR. DME .GT. DMETAB(NDME)) . STOP 'INTERP_SCAT_TABLE: particle Dme out of range.' C See if wavenumber is within last wavenumber grid, otherwise C find the grid location and interpolation factor for WAVENO NEWIW = .FALSE. IF (WAVENO .LT. WAVETAB(IW0) .OR. WAVENO .GT. WAVETAB(IW1)) THEN IL=1 IU=NWAVE DO WHILE (IU-IL .GT. 1) IM = (IU+IL)/2 IF (WAVENO .GE. WAVETAB(IM)) THEN IL = IM ELSE IU = IM ENDIF ENDDO IW0 = MAX(IL,1) IW1 = IW0+1 NEWIW = .TRUE. ENDIF IF (DME .NE. OLDDME) THEN C Find the grid location and interpolation factor for DME IL=1 IU=NDME DO WHILE (IU-IL .GT. 1) IM = (IU+IL)/2 IF (DME .GE. DMETAB(IM)) THEN IL = IM ELSE IU = IM ENDIF ENDDO ID = MAX(IL,1) FDME = (DME-DMETAB(ID))/(DMETAB(ID+1)-DMETAB(ID)) FLDME = LOG(DME/DMETAB(ID))/LOG(DMETAB(ID+1)/DMETAB(ID)) ENDIF IF (DME .NE. OLDDME .OR. NEWIW) THEN C If not the same Dme or a new wavenumber grid, then C linearly interpolate omega and g and log interpolate extinction EXT0 = EXP( (1-FLDME)*LOG(TABEXTINCT(IW0,ID)) . + FLDME*LOG(TABEXTINCT(IW0,ID+1)) ) EXT1 = EXP( (1-FLDME)*LOG(TABEXTINCT(IW1,ID)) . + FLDME*LOG(TABEXTINCT(IW1,ID+1)) ) ALB0 = (1-FDME)*TABSSALB(IW0,ID) + FDME*TABSSALB(IW0,ID+1) ALB1 = (1-FDME)*TABSSALB(IW1,ID) + FDME*TABSSALB(IW1,ID+1) ASYM0 = (1-FDME)*TABASYM(IW0,ID) + FDME*TABASYM(IW0,ID+1) ASYM1 = (1-FDME)*TABASYM(IW1,ID) + FDME*TABASYM(IW1,ID+1) ENDIF C Linearly interpolate the scattering properties in wavenumber FWAV = (WAVENO-WAVETAB(IW0))/(WAVETAB(IW1)-WAVETAB(IW0)) F = 1-FWAV EXTINCT = F*EXT0 + FWAV*EXT1 SSALB = F*ALB0 + FWAV*ALB1 ASYM = F*ASYM0 + FWAV*ASYM1 OLDDME = DME RETURN END SUBROUTINE INTERP_SCAT_TABLE3 (MU, WAVENO, DME, . EXTINCT, SSALB, ASYM, . PHI1UP, PHI1DN, PHI2UP, PHI2DN, . NMU, MUTAB, NDME, DMETAB, NWAVE, WAVETAB, . TABEXTINCT, TABSSALB, TABASYM, . TABPHI1UP, TABPHI1DN, TABPHI2UP, TABPHI2DN) C Interpolates the scattering properties from the table for C a particular observation angle, wavenumber, and particle size. C Does a trilinear interpolation, but optimized for fixed viewing C angle and particle size and slowly varying wavenumber. If the C DME and MU are the same as last time then we can just linearly C interpolate in wavenumber between stored scattering values. C If the DME or MU has changed then we bilinearly interpolate between C the DMETAB and MUTAB grid lines for each of the two wavenumber C grid lines. IMPLICIT NONE REAL MU, WAVENO, DME REAL EXTINCT, SSALB, ASYM, PHI1UP, PHI1DN, PHI2UP, PHI2DN INTEGER NMU, NDME, NWAVE REAL MUTAB(NMU), DMETAB(NDME), WAVETAB(NWAVE) REAL TABEXTINCT(NWAVE,NDME), TABSSALB(NWAVE,NDME) REAL TABASYM(NWAVE,NDME) REAL TABPHI1UP(NMU,NWAVE,NDME), TABPHI1DN(NMU,NWAVE,NDME) REAL TABPHI2UP(NMU,NWAVE,NDME), TABPHI2DN(NMU,NWAVE,NDME) INTEGER IW0, IW1, ID, IMU, IL, IU, IM LOGICAL NEWIW REAL FWAV, FDME, FLDME, FMU, F1, F2, F3, F4 REAL OLDMU, OLDDME, EXT0, EXT1, ALB0, ALB1, ASYM0, ASYM1 REAL PHI1UP0, PHI1UP1, PHI1DN0, PHI1DN1 REAL PHI2UP0, PHI2UP1, PHI2DN0, PHI2DN1 SAVE IW0, IW1, IMU,ID, OLDMU, OLDDME, FDME, FLDME, FMU SAVE EXT0, EXT1, ALB0, ALB1, ASYM0, ASYM1 SAVE PHI1UP0, PHI1UP1, PHI1DN0, PHI1DN1 SAVE PHI2UP0, PHI2UP1, PHI2DN0, PHI2DN1 DATA IW0/1/, IW1/2/ C Check that parameter are in range of table IF (WAVENO .LT. WAVETAB(1) .OR. . WAVENO .GT. WAVETAB(NWAVE)) THEN PRINT *, WAVENO,WAVETAB(1),WAVETAB(NWAVE) STOP 'INTERP_SCAT_TABLE3: wavenumber out of range.' ENDIF IF (DME .LT. DMETAB(1) .OR. DME .GT. DMETAB(NDME)) THEN WRITE(*,*) 'Dme: ', DME, DMETAB(1), DMETAB(NDME) STOP 'INTERP_SCAT_TABLE: particle Dme out of range.' ENDIF c IF (MU .LT. MUTAB(1) .OR. MU .GT. MUTAB(NMU)) c . STOP 'INTERP_SCAT_TABLE3: viewing angle mu out of range.' C See if wavenumber is within last wavenumber grid, otherwise C find the grid location and interpolation factor for WAVENO NEWIW = .FALSE. IF (WAVENO .LT. WAVETAB(IW0) .OR. WAVENO .GT. WAVETAB(IW1)) THEN IL=1 IU=NWAVE DO WHILE (IU-IL .GT. 1) IM = (IU+IL)/2 IF (WAVENO .GE. WAVETAB(IM)) THEN IL = IM ELSE IU = IM ENDIF ENDDO IW0 = MAX(IL,1) IW1 = IW0+1 NEWIW = .TRUE. ENDIF IF (MU .NE. OLDMU) THEN C Find the grid location and interpolation factor for MU IL=1 IU=NMU DO WHILE (IU-IL .GT. 1) IM = (IU+IL)/2 IF (MU .GE. MUTAB(IM)) THEN IL = IM ELSE IU = IM ENDIF ENDDO IMU = MAX(IL,1) FMU = (MU-MUTAB(IMU))/(MUTAB(IMU+1)-MUTAB(IMU)) ENDIF IF (DME .NE. OLDDME) THEN C Find the grid location and interpolation factor for DME IL=1 IU=NDME DO WHILE (IU-IL .GT. 1) IM = (IU+IL)/2 IF (DME .GE. DMETAB(IM)) THEN IL = IM ELSE IU = IM ENDIF ENDDO ID = MAX(IL,1) FDME = (DME-DMETAB(ID))/(DMETAB(ID+1)-DMETAB(ID)) FLDME = LOG(DME/DMETAB(ID))/LOG(DMETAB(ID+1)/DMETAB(ID)) ENDIF C If not the same Dme and mu, then bilinearly interpolate things C Logarithmically interpolate extinction IF (DME .NE. OLDDME .OR. MU .NE. OLDMU .OR. NEWIW) THEN F1 = (1-FMU)*(1-FDME) F2 = (1-FMU)*FDME F3 = FMU*(1-FDME) F4 = FMU*FDME EXT0 = EXP( (1-FLDME)*LOG(TABEXTINCT(IW0,ID)) . + FLDME*LOG(TABEXTINCT(IW0,ID+1)) ) EXT1 = EXP( (1-FLDME)*LOG(TABEXTINCT(IW1,ID)) . + FLDME*LOG(TABEXTINCT(IW1,ID+1)) ) ALB0 = (1-FDME)*TABSSALB(IW0,ID) + FDME*TABSSALB(IW0,ID+1) ALB1 = (1-FDME)*TABSSALB(IW1,ID) + FDME*TABSSALB(IW1,ID+1) ASYM0 = (1-FDME)*TABASYM(IW0,ID) + FDME*TABASYM(IW0,ID+1) ASYM1 = (1-FDME)*TABASYM(IW1,ID) + FDME*TABASYM(IW1,ID+1) PHI1UP0 = F1*TABPHI1UP(IMU,IW0,ID) + F2*TABPHI1UP(IMU,IW0,ID+1) . + F3*TABPHI1UP(IMU+1,IW0,ID) + F4*TABPHI1UP(IMU+1,IW0,ID+1) PHI1UP1 = F1*TABPHI1UP(IMU,IW1,ID) + F2*TABPHI1UP(IMU,IW1,ID+1) . + F3*TABPHI1UP(IMU+1,IW1,ID) + F4*TABPHI1UP(IMU+1,IW1,ID+1) PHI1DN0 = F1*TABPHI1DN(IMU,IW0,ID) + F2*TABPHI1DN(IMU,IW0,ID+1) . + F3*TABPHI1DN(IMU+1,IW0,ID) + F4*TABPHI1DN(IMU+1,IW0,ID+1) PHI1DN1 = F1*TABPHI1DN(IMU,IW1,ID) + F2*TABPHI1DN(IMU,IW1,ID+1) . + F3*TABPHI1DN(IMU+1,IW1,ID) + F4*TABPHI1DN(IMU+1,IW1,ID+1) PHI2UP0 = F1*TABPHI2UP(IMU,IW0,ID) + F2*TABPHI2UP(IMU,IW0,ID+1) . + F3*TABPHI2UP(IMU+1,IW0,ID) + F4*TABPHI2UP(IMU+1,IW0,ID+1) PHI2UP1 = F1*TABPHI2UP(IMU,IW1,ID) + F2*TABPHI2UP(IMU,IW1,ID+1) . + F3*TABPHI2UP(IMU+1,IW1,ID) + F4*TABPHI2UP(IMU+1,IW1,ID+1) PHI2DN0 = F1*TABPHI2DN(IMU,IW0,ID) + F2*TABPHI2DN(IMU,IW0,ID+1) . + F3*TABPHI2DN(IMU+1,IW0,ID) + F4*TABPHI2DN(IMU+1,IW0,ID+1) PHI2DN1 = F1*TABPHI2DN(IMU,IW1,ID) + F2*TABPHI2DN(IMU,IW1,ID+1) . + F3*TABPHI2DN(IMU+1,IW1,ID) + F4*TABPHI2DN(IMU+1,IW1,ID+1) ENDIF C Linearly interpolate the scattering properties in wavenumber FWAV = (WAVENO-WAVETAB(IW0))/(WAVETAB(IW1)-WAVETAB(IW0)) F1 = 1-FWAV EXTINCT = F1*EXT0 + FWAV*EXT1 SSALB = F1*ALB0 + FWAV*ALB1 ASYM = F1*ASYM0 + FWAV*ASYM1 PHI1UP = F1*PHI1UP0 + FWAV*PHI1UP1 PHI1DN = F1*PHI1DN0 + FWAV*PHI1DN1 PHI2UP = F1*PHI2UP0 + FWAV*PHI2UP1 PHI2DN = F1*PHI2DN0 + FWAV*PHI2DN1 OLDDME = DME OLDMU = MU RETURN END SUBROUTINE READ_ABS_PROFILE (ABSFILE, BINARYABSFILE, . ABSNU1, ABSNU2, ABSDELNU, . NABSNU, MAXNZ, MAXABSNU, NLEV, HEIGHT, TEMP, ABSPROF) C Reads in the gaseous absorption file which contains layer C optical depths for NLEV-1 layers and NABSNU wavenumbers. C The height (km) and temperature (K) of the profile are also returned. C The file may be in ascii text or binary format. IMPLICIT NONE INTEGER NABSNU, MAXNZ, MAXABSNU, NLEV LOGICAL BINARYABSFILE REAL*8 ABSNU1, ABSNU2, ABSDELNU REAL HEIGHT(*), TEMP(*), ABSPROF(MAXNZ,*) CHARACTER ABSFILE*(*) INTEGER I, J REAL NU IF (BINARYABSFILE) THEN OPEN (UNIT=1,FILE=ABSFILE, STATUS='OLD',FORM='UNFORMATTED') READ (1) ABSNU1, ABSNU2, ABSDELNU, NABSNU C WRITE (*,*) ABSNU1, ABSNU2, ABSDELNU, NABSNU IF (NABSNU .GT. MAXABSNU) STOP 'RTSPEC: MAXABSNU exceeded' READ (1) NLEV c WRITE (*,*) NLEV IF (NLEV .GT. MAXNZ) STOP 'RTSPEC: MAXNZ exceeded' READ (1) (HEIGHT(I), I=1,NLEV) c WRITE (*,'(52(1X,F5.1),:)') (HEIGHT(I), I=1,NLEV) READ (1) (TEMP(I), I=1,NLEV) c WRITE (*,'(52(1X,F5.1),:)') (TEMP(I), I=1,NLEV) DO J = 1, NABSNU READ (1) (ABSPROF(I,J), I=1,NLEV-1) ENDDO CLOSE (1) ELSE OPEN (UNIT=1, FILE=ABSFILE, STATUS='OLD') READ (1,*) READ (1,*) ABSNU1, ABSNU2, ABSDELNU, NABSNU IF (NABSNU .GT. MAXABSNU) THEN PRINT *, 'NABSNU',NABSNU,MAXABSNU STOP 'RTSPEC: MAXABSNU exceeded' ENDIF READ (1,*) NLEV IF (NLEV .GT. MAXNZ) STOP 'RTSPEC: MAXNZ exceeded' READ (1,*) (HEIGHT(I), I=1,NLEV) READ (1,*) (TEMP(I), I=1,NLEV) READ (1,*) DO J = 1, NABSNU READ (1,*) NU, (ABSPROF(I,J), I=1,NLEV-1) ENDDO CLOSE (1) ENDIF RETURN END SUBROUTINE INPUT_PARAMETERS (NUMSPEC, MAXSPEC, TWOSIDEBAND, . OUTNU1, OUTNU2, OUTDELNU, BANDCENTER, BANDOFFSET, MUOBS, . SFCEMIS1, SFCEMIS2, OUTUNITS, OUTAVERAGE, XO, . ABSFILE, NSCATTAB, MAXSCAT, SCATFILE, ZOBS, SFCTEMP, NCLDLAY, . MAXNZ, IWP, DME, ISCATTAB, ZTOP, RTMODEL, ABSTYPE, . BINARYABSFILE, OUTPUTFILE) C This subroutine just reads all the user input needed to run C the RT. IMPLICIT NONE INTEGER NUMSPEC, NSCATTAB, NCLDLAY, ISCATTAB(*) INTEGER MAXSPEC, MAXSCAT, MAXNZ INTEGER I, L, M LOGICAL TWOSIDEBAND, BINARYABSFILE REAL*8 OUTDELNU(*), OUTNU1(*), OUTNU2(*), XO REAL*8 BANDCENTER(*), BANDOFFSET(*) REAL MUOBS(*), SFCEMIS1(*), SFCEMIS2(*), SFCTEMP, ZOBS REAL ZTOP, IWP(*), DME(*) CHARACTER*100 OUTPUTFILE, ABSFILE(*), SCATFILE(*) CHARACTER*2 OUTAVERAGE(*) CHARACTER*1 RTMODEL, ABSTYPE, OUTUNITS(*) WRITE (*,*) WRITE (*,*) 'RTSPEC radiative transfer code' WRITE(*,*) 'Number of spectral segments to output' READ (*,*) NUMSPEC WRITE (*,*) NUMSPEC WRITE (*,*) 'Using a double side band radiometer? (T or F)' READ (*,*) TWOSIDEBAND WRITE (*,*) TWOSIDEBAND IF (NUMSPEC .GT. MAXSPEC) STOP 'RTSPEC: MAXSPEC exceeded' DO M = 1, NUMSPEC WRITE (*,'(A,I3)') ' Segment: ', M IF(TWOSIDEBAND) THEN WRITE (*,*) ' Band center, offset, and band pass' READ (*,*) BANDCENTER(M), BANDOFFSET(M), OUTDELNU(M) WRITE (*,*) BANDCENTER(M), BANDOFFSET(M), OUTDELNU(M) ELSE WRITE(*,*) ' Starting, ending, and delta wavenumber (cm^-1)' READ (*,*) OUTNU1(M), OUTNU2(M), OUTDELNU(M) WRITE (*,*) OUTNU1(M), OUTNU2(M), OUTDELNU(M) ENDIF WRITE(*,*) ' Viewing angle cosine for this segment' READ (*,*) MUOBS(M) WRITE (*,*) MUOBS(M) WRITE (*,*) ' Starting and ending emissivity for this segment' READ (*,*) SFCEMIS1(M),SFCEMIS2(M) WRITE (*,*) SFCEMIS1(M),SFCEMIS2(M) IF (OUTNU1(M) .GT. OUTNU2(M)) . STOP 'RTSPEC: ending wavenumber must be greater than starting' WRITE (*,*) . ' Output units (T-Tb(K) or R-radiance (mW/(m^2 ster cm^-1))' READ (*,'(A)') OUTUNITS(M) WRITE (*,*) OUTUNITS(M) WRITE (*,*) . ' Output bandpass averaging (R-rectangular, T-triangular,' WRITE (*,*) . ' CR-convolution rectangular, CT-convolution triangular,' WRITE (*,*) . ' CS-convolution sinc^2)' READ (*,'(A)') OUTAVERAGE(M) WRITE (*,*) OUTAVERAGE(M) IF (OUTAVERAGE(M).EQ.'CR' .OR. OUTAVERAGE(M).EQ.'cr' .OR. . OUTAVERAGE(M).EQ.'CT' .OR. OUTAVERAGE(M).EQ.'ct' .OR. . OUTAVERAGE(M).EQ.'CS' .OR. OUTAVERAGE(M).EQ.'cs') THEN IF (NUMSPEC .GT. 1) . STOP 'RTSPEC: can only have one spectrum if convolving' WRITE (*,*) . ' The FFT window half width (cm) (usually ~.5/OUTDELNU)' READ (*,*) XO WRITE (*,*) XO ENDIF WRITE (*,*) ' Absorption profile file name' READ (*,'(A)') ABSFILE(M) WRITE (*,*) ABSFILE(M) ENDDO WRITE(*,*) 'Number of scattering tables' READ (*,*) NSCATTAB WRITE (*,*) NSCATTAB IF (NSCATTAB .GT. MAXSCAT) STOP 'RTSPEC: MAXSCAT exceeded' WRITE(*,*) 'Input scattering table file name for each table' DO I = 1, NSCATTAB READ (*,'(A)') SCATFILE(I) WRITE (*,*) SCATFILE(I) ENDDO WRITE(*,*) 'Observation level (km)' READ (*,*) ZOBS WRITE (*,*) ZOBS WRITE(*,*) 'Surface Temperature (K) (-1 for same as atmosphere)' READ (*,*) SFCTEMP WRITE (*,*) SFCTEMP WRITE(*,*) 'Number of cloud layers' READ (*,*) NCLDLAY WRITE (*,*) NCLDLAY IF (NCLDLAY .GT. MAXNZ) STOP 'RTSPEC: MAXNZ exceeded by NCLDLAY' WRITE (*,*) 'IWP/LWP (g/m^2), Dme (micron), scattering table #', . ' for each layer' DO L = 1, NCLDLAY READ (*,*) IWP(L), DME(L), ISCATTAB(L) WRITE (*,*) IWP(L), DME(L), ISCATTAB(L) ENDDO WRITE (*,*) 'Input cloud top height (km)' READ (*,*) ZTOP WRITE (*,*) ZTOP WRITE (*,*) . 'Single-scattering, Eddington or hybrid RT model (S, E or H)' READ (*,*) RTMODEL WRITE (*,*) RTMODEL WRITE (*,*) 'Absorption file type (A-ascii text, B-binary)' READ (*,'(A1)') ABSTYPE WRITE (*,*) ABSTYPE BINARYABSFILE = ABSTYPE .EQ. 'B' .OR. ABSTYPE .EQ. 'b' WRITE(*,*) 'Output file name' READ (*,'(A)') OUTPUTFILE WRITE (*,*) OUTPUTFILE RETURN END SUBROUTINE OUTPUT_FILE (OUTPUTFILE, NUMSPEC, TWOSIDEBAND, OUTNU1, . OUTNU2, OUTDELNU, BANDCENTER, BANDOFFSET, OUTUNITS, . OUTAVERAGE, XO, ABSFILE, NSCATTAB, SCATFILE, MUOBS, ZOBS, . SFCTEMP, SFCEMIS1, SFCEMIS2, ZTOP, ZBOT, NCLDLAY, IWP, DME, . ISCATTAB, RTMODEL, NONU, ONU, OUTRAD) C This subroutine just opens the output file and write the header C and spectra. IMPLICIT NONE INTEGER NONU, NUMSPEC, NSCATTAB, NCLDLAY, ISCATTAB(*) INTEGER I, L, M LOGICAL TWOSIDEBAND, OUTTB, OUTTRIANG REAL*8 OUTDELNU(*), OUTNU1(*), OUTNU2(*) REAL*8 BANDCENTER(*), BANDOFFSET(*), ONU(*), XO REAL MUOBS(*), ZOBS, SFCTEMP, SFCEMIS1(*), SFCEMIS2(*) REAL ZTOP, ZBOT REAL IWP(*), DME(*) REAL OUTRAD(*) CHARACTER*100 OUTPUTFILE, ABSFILE(*), SCATFILE(*) CHARACTER*2 OUTAVERAGE(*) CHARACTER*1 RTMODEL, OUTUNITS(*) CHARACTER*48 UNITS, AVERAGE OPEN (UNIT=4, STATUS='UNKNOWN', FILE=OUTPUTFILE) WRITE (4,'(A)') '! RTSPEC Output' WRITE (4,'(A)') '! ' DO M = 1, NUMSPEC WRITE (4,'(A,I2)') '! Segment: ',M IF (.NOT.TWOSIDEBAND) THEN WRITE (4,'(A,3(1X,F9.3)))') '! Wavenumber range and step = ', . OUTNU1(M), OUTNU2(M), OUTDELNU(M) ELSE WRITE (4,'(A,3(1X,F9.3)))') . '! Band Center, offset, and width = ', . BANDCENTER(M),BANDOFFSET(M), OUTDELNU(M) ENDIF IF (OUTUNITS(M).EQ.'T'.OR.OUTUNITS(M).EQ.'t') THEN UNITS = '(K)' ELSE UNITS = '(mW/(m^2 ster cm^-1))' ENDIF IF (OUTAVERAGE(M).EQ.'T'.OR.OUTAVERAGE(M).EQ.'t') THEN AVERAGE = 'Triangular' ELSEIF (OUTAVERAGE(M).EQ.'R'.OR.OUTAVERAGE(M).EQ.'r') THEN AVERAGE = 'Rectangular' ELSEIF (OUTAVERAGE(M).EQ.'CR'.OR.OUTAVERAGE(M).EQ.'cr') THEN AVERAGE = 'Convolution: rectangular FFT window, XO = ' ELSEIF (OUTAVERAGE(M).EQ.'CT'.OR.OUTAVERAGE(M).EQ.'ct') THEN AVERAGE = 'Convolution: triangular FFT window, XO = ' ELSEIF (OUTAVERAGE(M).EQ.'CS'.OR.OUTAVERAGE(M).EQ.'cs') THEN AVERAGE = 'Convolution: sinc^2 FFT window, XO = ' ENDIF WRITE (4,'(A,A22))') '! Output units = ',UNITS IF (OUTAVERAGE(M).EQ.'CR'.OR.OUTAVERAGE(M).EQ.'cr' .OR. . OUTAVERAGE(M).EQ.'CT'.OR.OUTAVERAGE(M).EQ.'ct' .OR. . OUTAVERAGE(M).EQ.'CS'.OR.OUTAVERAGE(M).EQ.'cs') THEN WRITE (4,'(A,A42,F8.3))') . '! Output averaging = ',AVERAGE,XO ELSE WRITE (4,'(A,A12))') . '! Output averaging = ',AVERAGE ENDIF WRITE (4,'(A,A55)') '! Absorption file = ',ABSFILE(M) WRITE (4,'(A,F7.4)') . '! Observation angle mu = ',MUOBS(M) WRITE (4,'(A,F5.3,2X,F5.3)') . '! Surface emissivity range = ',SFCEMIS1(M),SFCEMIS2(M) ENDDO DO I = 1, NSCATTAB WRITE (4,'(A,I2,A,A55)') . '! Scattering table ',I,' = ',SCATFILE(I) ENDDO WRITE (4,'(A,F6.2,A,F6.2)') . '! Zobs = ',ZOBS,' Surface temperature = ',SFCTEMP WRITE (4,'(2(A,F7.3),A,F6.2)') . '! Cloud Ztop = ',ZTOP, ' Zbot = ',ZBOT DO L = 1, NCLDLAY WRITE (4,'(A,I3,A,F7.2,A,F6.1,A,I2)') '! Layer:',L, . ' IWP=',IWP(L), ' Dme=',DME(L), ' Iscat=',ISCATTAB(L) ENDDO WRITE (4,'(A,A1)') '! RT model = ',RTMODEL WRITE (4,'(A)') '! Waveno Tb/Radiance' IF (OUTUNITS(1).EQ.'T'.OR.OUTUNITS(1).EQ.'t') THEN DO I = 1,NONU WRITE (4,'(F9.3,1X,F8.3)') ONU(I), OUTRAD(I) ENDDO ELSE DO I = 1,NONU WRITE (4,'(F9.3,1X,E12.5)') ONU(I), OUTRAD(I) ENDDO ENDIF CLOSE (4) RETURN END SUBROUTINE AVG_SPECTRUM (MAXSPEC, NUMSPEC, NWAVENO, NOUTNU, . WAVENO, OUTNU, RADOBS, OUTRAD, OUTUNITS, OUTAVERAGE, . NBAND, OUTNU1, OUTNU2, BANDCENTER, BANDOFFSET, OUTDELNU) C This subroutine averages the high resolution RT over the C various segments either using rectangular or triangular C averaging. IMPLICIT NONE INTEGER NUMSPEC, NWAVENO(*), NOUTNU, NBAND, MAXSPEC REAL RADOBS(2,MAXSPEC,*), OUTRAD(*) REAL*8 WAVENO(2,MAXSPEC,*), OUTNU(*) REAL*8 OUTNU1(*), OUTNU2(*), OUTDELNU(*) REAL*8 BANDCENTER(*), BANDOFFSET(*) CHARACTER*1 OUTUNITS(*) CHARACTER*2 OUTAVERAGE(*) LOGICAL OUTTB, OUTTRIANG INTEGER I, J, M, IBAND REAL*8 SUMRAD1, SUMRAD2, SUMWT1, SUMWT2 REAL*8 INVOUTDELNU, ONU1, ONU2, WT1, WT2 REAL OUTRAD0 I = 1 C Loop over each segment of spectra DO M = 1,NUMSPEC IF (OUTAVERAGE(M).EQ.'T' .OR. OUTAVERAGE(M).EQ.'t') THEN OUTTRIANG = .TRUE. ELSE OUTTRIANG = .FALSE. ENDIF IF (OUTUNITS(M).EQ.'T' .OR. OUTUNITS(M).EQ.'t') THEN OUTTB = .TRUE. ELSE OUTTB = .FALSE. ENDIF C Loop if have double side band DO IBAND = 1,NBAND C Need to find the Tb for each side band IF(NBAND.EQ.2) THEN OUTNU1(M) = BANDCENTER(M)+(2.*DBLE(IBAND)-3.)* . BANDOFFSET(M) OUTNU2(M) = OUTNU1(M) ENDIF C Set up for the output spectral averaging SUMRAD1 = 0.0D0 SUMWT1 = 0.0D0 SUMRAD2 = 0.0D0 SUMWT2 = 0.0D0 IF (OUTTRIANG) THEN INVOUTDELNU = 1.0D0/OUTDELNU(M) ONU1 = OUTNU1(M)-OUTDELNU(M) ELSE ONU1 = OUTNU1(M)-0.5D0*OUTDELNU(M) ENDIF ONU2 = ONU1 + OUTDELNU(M) C Loop over a segment of wavenumbers input in the absorption file DO J = 1,NWAVENO(M) C Do the spectral window averaging of the radiances C If past window for current output wavenumber (ONU1) then do C the weighted average, output the radiance, and set up for next. IF ((WAVENO(IBAND,M,J)-ONU2) .GE. (-.00001D0) . .OR. J .EQ. NWAVENO(M)) THEN IF(NWAVENO(M).EQ.1) THEN OUTRAD(I) = RADOBS(IBAND,M,J) ELSE IF (SUMWT1 .LE. 0.0) . STOP 'RTSPEC: zero weight in averaging' OUTRAD(I) = SUMRAD1/SUMWT1 ENDIF IF (OUTTRIANG) THEN OUTNU(I) = ONU1 ELSE OUTNU(I) = 0.5*(ONU1+ONU2) ENDIF IF (OUTTB) THEN OUTRAD(I) = 1.4388*OUTNU(I)/ . LOG(1.0+(1.1911E-8*OUTNU(I)**3)/OUTRAD(I)) ELSE OUTRAD(I) = 1000*OUTRAD(I) ENDIF ONU1 = ONU2 ONU2 = ONU1 + OUTDELNU(M) SUMRAD1 = SUMRAD2 SUMWT1 = SUMWT2 SUMRAD2 = 0.0D0 SUMWT2 = 0.0D0 IF ((.NOT.OUTTRIANG.OR.OUTNU(I).GE.OUTNU1(M)).AND. . (NBAND.EQ.1)) I = I + 1 ENDIF IF (OUTTRIANG) THEN C Triangular bandpass averaging: C Do the weighted average of the radiance: keep two sums, one C for each side of triangular weighting function WT2 = (WAVENO(IBAND,M,J)-ONU1)*INVOUTDELNU WT1 = 1.0D0 - WT2 SUMRAD1 = SUMRAD1 + WT1*RADOBS(IBAND,M,J) SUMWT1 = SUMWT1 + WT1 SUMRAD2 = SUMRAD2 + WT2*RADOBS(IBAND,M,J) SUMWT2 = SUMWT2 + WT2 ELSE C Rectangular bandpass averaging: SUMRAD1 = SUMRAD1 + RADOBS(IBAND,M,J) SUMWT1 = SUMWT1 + 1.0D0 ENDIF C End of loop over wavenumbers ENDDO IF(NBAND.EQ.2) THEN IF(IBAND.EQ.1) OUTRAD0 = OUTRAD(I) IF (IBAND.EQ.2) THEN OUTNU(I) = SNGL(BANDCENTER(M)+BANDOFFSET(M)) OUTRAD(I) = (OUTRAD0+OUTRAD(I))/2 I = I + 1 ENDIF ENDIF C End of loop for sidebands ENDDO C End of loop over spectral segments ENDDO NOUTNU = I - 1 RETURN END SUBROUTINE CONVOLVE_SPECTRUM (MAXSPEC,NWAVENO, WAVENOIN,RADOBSIN, . OUTUNITS, OUTAVERAGE,XO, DELWAVENO, . OUTNU1, OUTNU2, OUTDELNU, . NOUTNU, OUTNU, OUTRAD) C Subroutine convolves spectra with either sinc, sinc^2, or the C convolution of sinc with a triangle function. C The convolution is done by using two FFT's: C FFT^-1(f(x)*g(x)) = conv(F(k),G(k)) C f(x) = FFT(F(k)) g(x) = FFT(G(k)) C thus g(x) is a window function like a boxcar or triangle. The FFT of C the input spectra, f(x), is multiplied by the window function. This C product is then FFTed back to produce the convolved spectrum. The C window function can be a boxcar (OUTAVERAGE=CR), a triangle C (OUTAVERAGE=CT), or a sinc^2 (out to first zero, OUTAVERAGE=CS). C The width of the window function is given by XO (+/-XO is where the C window function becomes zero). The FFT arrays are padded with zeros C to twice the next largest power of two greater than NWAVENO, to C prevent aliasing. The convolution is performed on radiance and then C converted to brightness temperature at the end. The convolution is C also performed on spectra with the absorption file resolution. Once C the convolution is complete the high resolution spectra is interpolated C to desired low resolution output wavenumbers. C The input spectrum is the first radiance spectrum stored in RADOBSIN C with NWAVENO points at WAVENOIN wavenumbers. The output spectrum C is returned in OUTRAD with NOUTNU points at OUTNU wavenumbers C (determined from OUTNU1, OUTNU2, OUTDELNU). IMPLICIT NONE INTEGER MAXSPEC, NWAVENO INTEGER NOUTNU REAL*8 XO, DELWAVENO, OUTNU1, OUTNU2, OUTDELNU REAL*8 WAVENOIN(2,MAXSPEC,NWAVENO) REAL RADOBSIN(2,MAXSPEC,NWAVENO) REAL*8 OUTNU(*) REAL OUTRAD(*) CHARACTER*1 OUTUNITS CHARACTER*2 OUTAVERAGE INTEGER ND PARAMETER (ND=16384) INTEGER I, J, N1, N2, NFFT REAL DATA(2*ND), FILTER REAL*8 WAVENO(ND), X, XMAX, AS C Calculate size of FFT array. Put the input spectrum in the FFT C array, padding at either end with zeros (also put zeros in Im part). NFFT = 2* 2**(INT(LOG(REAL(NWAVENO-1))/LOG(2.))+1) IF (NFFT .GT. ND) THEN STOP 'CONVOLVE_SPECTRUM: ND exceeded' ENDIF N1 = INT((NFFT-2*NWAVENO)/4) DO I = 1, NFFT DATA(2*I-1) = 0.0 DATA(2*I) = 0.0 ENDDO N2 = N1 + NWAVENO J = 1 DO I = N1+1,N2 DATA(2*I-1) = RADOBSIN(1,1,J) DATA(2*I) = 0.0 J = J + 1 ENDDO C Transform spectrum from wavenumber space to x-space signal CALL FFT(DATA, NFFT, +1) C Create filter function and multiply x-space signal by filter XMAX = 1.0/DELWAVENO AS = ACOS(-1.0D0)/XO DO I = 1, NFFT X = XMAX*DBLE(I-1)/DBLE(NFFT) X = MIN(X,XMAX-X) IF (OUTAVERAGE .EQ. 'CR') THEN IF (X .LE. XO) THEN FILTER = 1.0 ELSE FILTER = 0.0 ENDIF ELSE IF (OUTAVERAGE .EQ. 'CT') THEN IF (X .LE. XO) THEN FILTER = SNGL((XO-X)/XO) ELSE FILTER = 0.0 ENDIF ELSE IF (OUTAVERAGE .EQ. 'CS') THEN IF (X .EQ. 0.0) THEN FILTER = 1.0 ELSE IF (X .LE. XO) THEN FILTER = ( SIN(AS*X) / (AS*X) )**2 ELSE FILTER = 0.0 ENDIF ELSE PRINT *, 'CONVOLVE_SPECTRUM: Unknown type: ',OUTAVERAGE STOP ENDIF DATA(2*I-1) = DATA(2*I-1)*FILTER DATA(2*I) = DATA(2*I)*FILTER ENDDO C Transform from x-space back to radiance in k-space CALL FFT(DATA, NFFT, -1) C Divide by N since FFT routine doesn't for inverse J = 1 DO I = N1+1,N2 DATA(J) = DATA(2*I-1)/NFFT WAVENO(J) = WAVENOIN(1,1,J) J = J + 1 ENDDO C Interpolate to desired output wavenumbers C and convert to brightness temperature NOUTNU = INT((OUTNU2-OUTNU1)/OUTDELNU) + 1 DO I = 1, NOUTNU OUTNU(I) = OUTNU1+OUTDELNU*(I-1) CALL INTERPOLATE_PROFILE (NWAVENO, WAVENO, DATA, . OUTNU(I), OUTRAD(I)) IF (OUTUNITS .EQ. 'T' .OR. OUTUNITS .EQ. 't') THEN OUTRAD(I) = 1.4388*OUTNU(I)/ . LOG(1.0+(1.1911E-8*OUTNU(I)**3)/OUTRAD(I)) ELSE OUTRAD(I) = 1000*OUTRAD(I) ENDIF ENDDO RETURN END SUBROUTINE FFT (DATA,N,ISIGN) C Fast Fourier transform routine. First the input array C must be "bit-reversed". This is due to the fact that if one kept C breaking up discrete FFT in halves (Danielson-Lanczos) made up C of even and odd points, we would eventually be left with series of C one-point FFT which are just original series in bit-reversed order. C The algorithm goes backwards from this so we must begin with the C array bit-reversed so that we combine correctly. C Basically uses the method of combining adjacent blocks to form C 2 point transforms, then combining 2 point transforms to get 4 point C transforms, etc, until desired N point transform is reached. At each C step an N/n n-point transforms are produced which repeat with period n. C K is the index for the x-axis in Fourier space, F_k = F_k+n. C J is the index for the number N/n n-point transforms created C in each step. For each step there is a constant Wn = exp(2*pi*i/n) C that is used to go from n/2-point transform to n-point transform: C Fk,n = Fk,n/2,j + Wn^k*Fk,n/2,j+1 C and since Fk,n/2 = Fk+n/2,n/2 we can produced two n-point C transforms at k and k+n/2 since we know Wn^k = Wn^(k+n/2) and C Wn^k = -Wn^(k+n/2). Thus, two n-point transforms is calculated C from only one Wn calculation for two points in "k" given a pair C of n/2-point transforms with only one point in "k". Wn^k is C calculated using recurrence relations, so the J loop is inside C the K loop limiting the number of calculations. C The sequence is saved in a 2*N length array. 2*N since first spot C is the real part and second spot imaginary. For each step N/n C n-point transforms with k = 1 to n are stored in the array. There C is no need to store k = n+1 to N since n point transforms are periodic C in n. The first spot in array is k = 1 the next spot k = 2 up to C k = n for the 1st n-point transform. This is repeated for the C second n-point transform up to the N/n th n-point transform. C ie. for N = 4. C intial : F(k=1,n=1,j=1), F(k=1,n=1,j=2), F(k=1,n=1,j=3), F(k=1,n=1,j=4) C 1st step: F(k=1,n=2,j=1), F(k=2,n=2,j=1), F(k=1,n=2,j=2), F(k=2,n=2,j=2) C 2nd step: F(k=1,n=4,j=1), F(k=2,n=4,j=1), F(k=3,n=4,j=1), F(k=4,n=4,j=1) IMPLICIT NONE INTEGER N, ISIGN REAL DATA(*), TMPR, TMPI REAL*8 WR, WI, WTMP, THETA, ALPHA, BETA INTEGER I, IREV, M, J, K, M0, M1 INTEGER POWER, POWER2, POWER4 REAL*8 PI PARAMETER (PI=3.14159265359) C This is the bit reversal process IF (N .LE. 1) RETURN IREV=0 DO I = 0, N-1 IF (I .GT. IREV) THEN M0 = 2*I+1 M1 = 2*IREV+1 TMPR = DATA(M0) TMPI = DATA(M0+1) DATA(M0) = DATA(M1) DATA(M0+1) = DATA(M1+1) DATA(M1) = TMPR DATA(M1+1) = TMPI ENDIF M = N 50 CONTINUE M = M/2 IF (IREV.LT.M) GO TO 70 IREV = IREV - M IF (M .GT. 1) GOTO 50 70 IREV = IREV + M ENDDO POWER = 1 POWER2 = POWER*2 POWER4 = POWER*4 C Loop over the log2(N) combination steps 200 CONTINUE THETA = ISIGN*PI/POWER ALPHA = 2*(SIN(THETA/2))**2 BETA = SIN(THETA) WR = 1. WI = 0. M = 1 C Loop over the Fourier space index DO K = 1,POWER M0 = M M1 = M0 + POWER2 C Loop over the N/n n-point transforms summing adjacent blocks DO J = 1, N/POWER2 TMPR = SNGL(WR)*DATA(M1)-SNGL(WI)*DATA(M1+1) TMPI = SNGL(WR)*DATA(M1+1)+SNGL(WI)*DATA(M1) DATA(M1) = DATA(M0) - TMPR DATA(M1+1) = DATA(M0+1) - TMPI DATA(M0) = DATA(M0) + TMPR DATA(M0+1) = DATA(M0+1) + TMPI M0 = M0 + POWER4 M1 = M1 + POWER4 ENDDO M = M + 2 C Use recurssion to find new "Wk" WTMP = WR WR = WR - ALPHA*WR - BETA*WI WI = WI - ALPHA*WI + BETA*WTMP ENDDO POWER = POWER2 POWER2 = POWER4 POWER4 = 2*POWER4 IF (POWER .LT. N) GOTO 200 RETURN END FUNCTION RNDUP(X) IMPLICIT NONE REAL X, MOD INTEGER RNDUP RNDUP = INT(X) MOD = X-REAL(RNDUP) IF(MOD.GT..01) RNDUP=RNDUP+1 RETURN END SUBROUTINE INTERPOLATE_PROFILE (N, XO, YO, XN, YN) C Interpolates Xo,Yo data to Xn finding Yn IMPLICIT NONE INTEGER N REAL*8 XO(*), XN REAL YO(*), YN INTEGER I, IL, IU, IM REAL X c Find XO value that is closest to XN IL=0 IU=N+1 10 IF (IU-IL.GT.1) THEN IM=(IU+IL)/2 IF (XN.LT.XO(IM)) THEN IU=IM ELSE IL=IM ENDIF GO TO 10 ENDIF I=IL C If XN out of XO range, don't extrapolate, use min/max XO IF (I .LE. 0) THEN I = 1 X = 0.0 ELSE IF (I .GE. N) THEN I = N-1 X = 1.0 ELSE X = (XN-XO(I))/(XO(I+1)-XO(I)) ENDIF YN = YO(I) + (YO(I+1)-YO(I))*X RETURN END