!------------------------------------------------------------------------------ ! ! NAME: ! nceprtm_hirs_v2 ! ! PURPOSE: ! Program to compute HIRS brightness temperatures from NCEP reanalysis ! temperature and water vapor annual files containing 6-hourly data. ! Intended to be used on TOVS Pathfinder Radiance processing V2.2 ! ! Uses ozone climatology from Fortuin and Kelder, JGR, 1999, 103, 31709-? ! ! LANGUAGE: ! Fortran-90 ! ! MODULES: ! type_kinds: Module to hold specification kinds for variable ! declaration. ! ! file_utility: Module containing generic file utility routines ! ! error_handler: Module to define simple error codes and handle ! error conditions ! ! parameters: Module to hold RT model parameter constants ! ! initialize: Module for RT model initialisation. ! ! forward_model: Module containing the RT forward model function. ! ! spectral_coefficients: Module to hold the RT model spectral coefficients ! and their access routines. ! NOTE: This routine is USEd only to access a single ! data entity for test purposes. Does not need ! to be USEd to call the RT routines. ! misc_darren_code: Routines added by Darren Jackson ! ! INPUTS: ! satellite -> NOAA series satellite number ! UNITS: none ! DIMENSION: scalar ! TYPE: integer ! channel -> HIRS channel number ! UNITS: none ! DIMENSION: 1 x MAX_CHANNEL ! TYPE: integer ! NOTE: Loops for each desired channel number until ! you enter 99. You can enter any number of ! channels less than MAX_CHANNEL ! year -> Four digit year ! UNITS: none ! DIMENSION: scalar ! TYPE: integer ! jday -> Day of year ! UNITS: none ! DIMENSION: scalar ! TYPE: integer ! OUTPUTS: ! Binary file contain temperature and transmission profiles, radiance and ! brightness temperature data for selected HIRS channels and 26 predefined ! levels. Output file name has format NCEP.satellite.year.day.hour.RT. Example ! NCEP.N11.Y91.D200.06Z.RT represents RT calculations from NCEP reanalyis ! data for NOAA-11 on day 200 of 1991 at hour 6Z. ! ! CONTAINS: ! None. ! ! INCLUDE FILES: ! None. ! ! EXTERNALS: ! None. ! ! COMMON BLOCKS: ! None. ! ! FILES ACCESSED: ! Input: Binary, direct access profile data file. ! Output: Binary, direct access forward model output file. ! ! RESTRICTIONS: ! None. ! ! CREATION HISTORY: ! Written by: Paul van Delst, CIMSS/SSEC 20-Aug-2001 ! paul.vandelst@ssec.wisc.edu ! ! Modified for NCEP reanalysis profiles and HIRS coefficients: ! ! Darren Jackson ETL/CIRES Sep-2001 ! Darren.L.Jackson@noaa.gov ! ! Version 2 allows specify each channel at input rather than range of channels and ! changed the number angles used from original veriosn. DLJ,june 2002 ! ! !-------------------------------------------------------------------------- PROGRAM nceprtm_hirs_v2 !#----------------------------------------------------------------------------# !# -- MODULE USAGE -- # !#----------------------------------------------------------------------------# USE type_kinds USE error_handler USE file_utility USE parameters USE initialize USE misc_darren_code USE forward_model, ONLY: compute_rtm ! -- For access to IS_MICROWAVE_CHANNEL only. USE spectral_coefficients !#----------------------------------------------------------------------------# !# -- TYPE DECLARATIONS -- # !#----------------------------------------------------------------------------# ! --------------------------- ! Disable all implicit typing ! --------------------------- IMPLICIT NONE ! ------------------ ! Program parameters ! ------------------ ! -- Name CHARACTER( * ), PARAMETER :: PROGRAM_NAME = 'nceprtm_hirstb' ! -- Angle definitions ! INTEGER, PARAMETER :: N_ANGLES = 4 ! INTEGER, PARAMETER :: N_ANGLES = 1 ! -- Profile definitions INTEGER, PARAMETER :: NLON=144,NLAT=73,RLEVELS=26,TLEVELS=17,SLEVELS=8 INTEGER, PARAMETER :: NBANDS=17 INTEGER, PARAMETER :: N_LAYERS=25 LOGICAL, PARAMETER :: logint=.true. CHARACTER( * ), PARAMETER :: OFILE2 = 'NCEP_RT.TXT' CHARACTER( * ), PARAMETER :: O3_FILE = & '/data/dlj/tovs1b/pathfinder2/ncep_data/o3-mean_26level.dat' CHARACTER( * ), PARAMETER :: TRAN_FILE = 'POES_transmittance_coefficients' CHARACTER( * ), PARAMETER :: SPEC_FILE = 'POES_spectral_coefficients' CHARACTER( * ), PARAMETER :: DPATH = '/data/dlj/tovs1b/pathfinder2/ncep_data/' ! -- Other dimension parameters INTEGER, PARAMETER :: N_ABSORBERS = MAX_N_ABSORBERS INTEGER, PARAMETER :: N_PREDICTORS = MAX_N_PREDICTORS INTEGER, PARAMETER :: MAX_CHANNELS = 50 INTEGER, PARAMETER :: MAX_ANGLES = 10 ! -- Emissivity parameters REAL( fp_kind ), PARAMETER :: DEFAULT_MW_EMISSIVITY = 0.6_fp_kind REAL( fp_kind ), PARAMETER :: DEFAULT_IR_EMISSIVITY = 0.96_fp_kind ! -- Default solar angle secant (> 11.47 means no solar) REAL( fp_kind ), PARAMETER :: DEFAULT_SECANT_SOLAR_ANGLE = 12.0_fp_kind ! ----------------- ! Program variables ! ----------------- ! -- Error message CHARACTER( 128 ) :: message ! INTEGER :: COUNT ! -- Status variables INTEGER :: error_status INTEGER :: allocate_status INTEGER :: io_status ! -- Loop counters, dimensions and file lun INTEGER :: i, j, k, l, m, il, itime, irec INTEGER :: n_available_channels, l1, l2, n_channels, n_angles INTEGER :: temp_record_length, sh_record_length, tp_record_length, sf_record_length INTEGER :: record_number,out_record_length,v2_record_number INTEGER :: lun,tp_lun,sh_lun,sf_lun,ps_lun,o3_lun INTEGER :: lun_BIN,lun_BIN_10,lun_BIN_09,lun_BIN_08,lun_BIN_07 INTEGER :: instrument,satellite,inst_channel INTEGER, DIMENSION(1) :: idx INTEGER :: lun_Wylie,lun_O3,in_lun INTEGER :: jday,iday,imon,iyr4,iyr2,latidx,ilev,iband,ii,ihr,i_chan,jdx CHARACTER (1) :: dum CHARACTER (3) :: angle_name CHARACTER( 80 ) :: OFILE CHARACTER( 80 ) :: TP_FILE CHARACTER( 80 ) :: SH_FILE CHARACTER( 80 ) :: PS_FILE CHARACTER( 80 ) :: SF_FILE ! -- Secant angle arrays INTEGER :: i_angle ! REAL( fp_kind ) :: d_angle ! REAL( fp_kind ), DIMENSION( N_ANGLES ) :: view_angle = & ! (/ 0.0_fp_kind, 25.8419_fp_kind, 36.8699_fp_kind, 45.5730_fp_kind /) ! (/ 0.0_fp_kind /) ! REAL( fp_kind ), DIMENSION( N_ANGLES ) :: secant_view_angle ! REAL( fp_kind ), DIMENSION( N_ANGLES ) :: secant_solar_angle INTEGER, DIMENSION( MAX_CHANNELS ) :: channel REAL, DIMENSION( MAX_ANGLES ) :: view_angle ! -- Profile read array INTEGER ( kind=2 ), DIMENSION( NLON, NLAT ) :: isf_data INTEGER ( kind=2 ), DIMENSION( NLON, NLAT ) :: ips_data REAL, DIMENSION( NLON, NLAT, RLEVELS) :: temperature_grid REAL, DIMENSION( NLON, NLAT, RLEVELS) :: water_vapor_grid REAL, DIMENSION( NLAT, RLEVELS, NBANDS):: o3_climo REAL ( fp_kind ), DIMENSION( TLEVELS ) :: level_pressure = & (/ 10.,20.,30.,50.,70.,100.,150.,200.,250.,300.,400.,500.,600.,700.,850.,925.,1000. /) REAL ( fp_kind ), DIMENSION( RLEVELS ) :: level_pressure2 = & (/ 0.1, 0.2, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 7.0, 10., 20., 30., 50., 70., 100., 150., & 200.,250.,300.,400.,500.,600.,700.,850.,925.,1000. /) ! -- Data types for Hals ozone code ! REAL , DIMENSION(40) :: hal_pressure = & ! (/.1,.2,.5,1.0,1.5,2.0,3.0,4.0,5.0,7.0,10.0, & ! 15.,20.,25.,30.,50.,60.,70.,85.,100.,115., & ! 135.,150.,200.,250.,300.,350.,400.,430.,475., & ! 500.,570.,620.,670.,700.,780.,850.,920.,950.,1000. /) ! REAL , DIMENSION(40) :: o3mix ! REAL , DIMENSION(RLEVELS) :: o3mix2 LOGICAL , DIMENSION( RLEVELS ) :: mask REAL( fp_kind ), DIMENSION( RLEVELS-1 ) :: layer_pressure,layer_temperature REAL( fp_kind ), DIMENSION( RLEVELS-1 ) :: layer_water_vapor REAL( fp_kind ), DIMENSION( RLEVELS ) :: level_temperature2 REAL( fp_kind ), DIMENSION( RLEVELS, 2 ) :: level_aa REAL( fp_kind ), DIMENSION( RLEVELS-1, 2) :: layer_aa ! -- Number of channels processed for each profile ! INTEGER, DIMENSION( N_ANGLES ) :: n_channels_per_profile ! -- Allocatable arrays ! -- Profile data arrays INTEGER, DIMENSION( : ) , ALLOCATABLE :: n_channels_per_profile REAL( fp_kind ), DIMENSION( : ) , ALLOCATABLE :: secant_view_angle REAL( fp_kind ), DIMENSION( : ) , ALLOCATABLE :: secant_solar_angle REAL( fp_kind ), DIMENSION( : ) , ALLOCATABLE :: mw_emissivity REAL( fp_kind ), DIMENSION( : ) , ALLOCATABLE :: ir_emissivity REAL( fp_kind ), DIMENSION( :, : ), ALLOCATABLE :: level_p REAL( fp_kind ), DIMENSION( :, : ), ALLOCATABLE :: layer_p REAL( fp_kind ), DIMENSION( :, : ), ALLOCATABLE :: layer_t REAL( fp_kind ), DIMENSION( :, : ), ALLOCATABLE :: layer_w REAL( fp_kind ), DIMENSION( :, : ), ALLOCATABLE :: layer_o REAL( fp_kind ), DIMENSION( : ), ALLOCATABLE :: sfc_t REAL( fp_kind ), DIMENSION( : ), ALLOCATABLE :: surface_emissivity REAL( fp_kind ), DIMENSION( : ), ALLOCATABLE :: surface_reflectivity INTEGER, DIMENSION( : ), ALLOCATABLE :: channel_index REAL( fp_kind ), DIMENSION( :, : ), ALLOCATABLE :: tau REAL( fp_kind ), DIMENSION( :, : ), ALLOCATABLE :: flux_tau REAL( fp_kind ), DIMENSION( :, : ), ALLOCATABLE :: solar_tau REAL( fp_kind ), DIMENSION( : ), ALLOCATABLE :: upwelling_radiance REAL( fp_kind ), DIMENSION( : ), ALLOCATABLE :: brightness_temperature REAL( fp_kind ), DIMENSION( :, : ), ALLOCATABLE :: tau2 REAL( fp_kind ), DIMENSION( : ), ALLOCATABLE :: upr2 REAL( fp_kind ), DIMENSION( : ), ALLOCATABLE :: brt2 ! -- Stuff for emissivities REAL( fp_kind ) :: angle_modifier REAL(single) :: xlon,xlat,dtdp ! INTEGER :: count,count_rate,count_max ! ---------- ! Intrinsics ! ---------- INTRINSIC COS, & MIN, MAX, & REAL, & TRIM, & MOD !############################################################################## !############################################################################## !############################################################################## !# # !# -- INITIALIZE THE RTM -- # !# # !############################################################################## !############################################################################## !############################################################################## WRITE( *, '( /5x, "Initialising RTM...." )' ) error_status = initialize_rtm( tau_file = TRAN_FILE, spectral_file = SPEC_FILE) IF ( error_status /= SUCCESS ) THEN CALL display_message( PROGRAM_NAME, & 'Error initialzing RTM', & error_status ) STOP END IF !############################################################################## !############################################################################## !############################################################################## !# # !# -- SETUP DATA AND ARRAYS FOR CALCULATIONS -- # !# # !############################################################################## !############################################################################## !############################################################################## !#----------------------------------------------------------------------------# !# -- ALLOCATE ARRAYS FOR RTM MODEL -- # !#----------------------------------------------------------------------------# ! ---------------------------------- ! Get the number of channels read in ! from the coefficient data files ! ---------------------------------- !!! !!! CALL get_max_n_channels( n_channels ) !!! begin_channel = 1 !!! end_channel = n_channels ! Replace the above code with the following to process channel subsets ! -- Input Satellite ID and Instrument in_lun=get_lun() OPEN(in_lun,FILE='nceprtm_hirs_v2.input') ! WRITE(*,*) 'Input NOAA satellite number (5=TN,6=NOAA-6,etc.):' READ(in_lun,*) satellite ! WRITE(*,*) 'Input instrument type (1=HIRS,2=MSU,3=AMSU):' READ(in_lun,*) instrument n_channels=1 DO ! WRITE(*,*) 'Input channel number: (99 ends input)' READ(in_lun,*) channel(n_channels) IF(channel(n_channels) == 99) EXIT IF(channel(n_channels) > 20) THEN WRITE(*,*) 'ENTERED CHANNEL NUMBER EXCEEDS ATOVS LIMIT OF 20' WRITE(*,*) 'CHANNEL ENTERED = ',channel(n_channels) WRITE(*,*) 'TERMINATE PROGRAM' STOP ENDIF n_channels = n_channels + 1 IF(n_channels > MAX_CHANNELS) THEN WRITE(*,*) 'Number of channels entered exceeds MAX_CHANNELS' WRITE(*,*) 'MAX_CHANNELS = ',max_channels WRITE(*,*) 'Terminate program' STOP ENDIF ENDDO n_channels = n_channels - 1 n_angles = 1 DO ! WRITE(*,*) 'Input view angle: (99. ends input)' READ(in_lun,*) view_angle(n_angles) IF(view_angle(n_angles) == 99.) EXIT IF(view_angle(n_angles) > 60.) THEN WRITE(*,*) 'VIEW ANGLE = ',view_angle(n_angles),' TOO LARGE' WRITE(*,*) 'MAXIMUM VIEW ANGLE = 60 DEGRESS' WRITE(*,*) 'terminate program' STOP ENDIF IF(view_angle(n_angles) < 0.) THEN WRITE(*,*) 'VIEW ANGLE = ',view_angle(n_angles),' LESS THAN ZERO' WRITE(*,*) 'MINIMUM VIEW ANGLE = 0 DEGRESS' WRITE(*,*) 'terminate program' STOP ENDIF n_angles = n_angles + 1 IF(n_angles > MAX_ANGLES) THEN WRITE(*,*) 'Number of angles entered exceeds MAX_ANGLES' WRITE(*,*) 'MAX_ANGLES = ',max_angles WRITE(*,*) 'terminate program' STOP ENDIF ENDDO n_angles = n_angles - 1 ! WRITE(*,*) 'Input four digit year:' READ(in_lun,*) iyr4 ! WRITE(*,*) 'Input day of year:' READ(in_lun,*) jday CLOSE(in_lun) ! ---------------------------------- ! Allocate angle dependent vectors !----------------------------------- ALLOCATE( secant_view_angle( N_ANGLES ), & secant_solar_angle( N_ANGLES ), & mw_emissivity( N_ANGLES ), & ir_emissivity( N_ANGLES ), & STAT = allocate_status) IF ( allocate_status /= 0 ) THEN WRITE( message, '( "Error allocating angle ", & &"dependent arrays. STAT = ", i5 )' ) & allocate_status CALL display_message( PROGRAM_NAME, & TRIM( message ), & FAILURE ) STOP END IF ! --------------------------------------------------- ! Allocate the Forward model channel dependent arrays ! --------------------------------------------------- ALLOCATE( surface_emissivity( n_channels * N_ANGLES ), & ! Input, L*M surface_reflectivity( n_channels * N_ANGLES ), & ! Input, L*M channel_index( n_channels * N_ANGLES ), & ! Input, L*M n_channels_per_profile ( N_ANGLES ), & ! Input level_p( N_LAYERS, n_angles ), & ! Input layer_p( N_LAYERS, n_angles ), & ! Input layer_t( N_LAYERS, n_angles ), & ! Input layer_w( N_LAYERS, n_angles ), & ! Input layer_o( N_LAYERS, n_angles ), & ! Input sfc_t( n_angles ), & ! Input tau( N_LAYERS, n_channels * N_ANGLES ), & ! Output, K x L*M flux_tau( N_LAYERS, n_channels * N_ANGLES ), & ! Output, K x L*M solar_tau( N_LAYERS, n_channels * N_ANGLES ), & ! Output, K x L*M upwelling_radiance( n_channels * N_ANGLES ), & ! Output, L*M brightness_temperature( n_channels * N_ANGLES ), & ! Output, L*M tau2( N_LAYERS, n_channels), & ! Output, K x L*M upr2( n_channels ), & ! Output, L*M brt2( n_channels ), & ! Output, L*M STAT = allocate_status ) IF ( allocate_status /= 0 ) THEN WRITE( message, '( "Error allocating forward model channel ", & &"dependent arrays. STAT = ", i5 )' ) & allocate_status CALL display_message( PROGRAM_NAME, & TRIM( message ), & FAILURE ) STOP END IF !#----------------------------------------------------------------------------# !# -- FILL PROFILE INDEPENDENT INPUT ARRAYS -- # !#----------------------------------------------------------------------------# ! -- Compute delta angle value ! IF(N_ANGLES > 1) then ! d_angle = ( MAX_ANGLE - MIN_ANGLE ) / REAL( N_ANGLES-1, fp_kind ) ! ELSE ! d_angle = 0.0 ! ENDIF ! -- Initialise angle x channel counter il = 0 ! ---------------- ! Loop over angles ! ---------------- DO i_angle = 1, N_ANGLES ! -- Fill angle arrays ! view_angle( i_angle ) = MIN_ANGLE + ( REAL( i_angle-1, fp_kind ) & ! * d_angle ) secant_view_angle( i_angle ) = ONE / COS( view_angle( i_angle ) & * DEGREES_TO_RADIANS ) secant_solar_angle( i_angle ) = DEFAULT_SECANT_SOLAR_ANGLE ! -- Fill channel count array, i.e. do them all n_channels_per_profile( i_angle ) = n_channels ! ------------------ ! Loop over channels ! ------------------ DO l = 1, N_CHANNELS ! -- Increment angle x channel counter ! -- and set the channel index il = il + 1 inst_channel=channel(l) channel_index( il ) = NCEP_CHANNEL_INDEX(satellite,instrument,inst_channel) write(*,*) 'channel index = ',channel_index(il) ! -- Assign pretend surface emissivity and reflectivity ! -- For uW, r = specular; for IR, r = isotropic ! -- The angle modifier is just something to provide ! -- a little bit of angular variation in the surface ! -- emissivities and reflectivities. angle_modifier = (COS( view_angle( i_angle ) * & DEGREES_TO_RADIANS ))**(0.1_fp_kind) IF ( is_microwave_channel( channel_index( il ) ) == 1 ) THEN surface_emissivity( il ) = DEFAULT_MW_EMISSIVITY * angle_modifier surface_reflectivity( il ) = ONE - surface_emissivity( il ) mw_emissivity( i_angle ) = surface_emissivity( il ) ELSE surface_emissivity( il ) = DEFAULT_IR_EMISSIVITY * angle_modifier surface_reflectivity( il ) = ( ONE - surface_emissivity( il ) ) / PI ir_emissivity( i_angle ) = surface_emissivity( il ) END IF END DO END DO !############################################################################## !############################################################################## !############################################################################## !# # !# -- OPEN THE INPUT PROFILE DATA FILE -- # !# # !############################################################################## !############################################################################## !############################################################################## ! ------------------------------ ! Set the record length in bytes ! ------------------------------ sf_record_length = NLAT * NLON * n_bytes_for_Short_kind tp_record_length = NLAT * NLON * RLEVELS * 4 sh_record_length = NLAT * NLON * RLEVELS * 4 ! -------------------------------------------- ! Open temperature and specific humidity files ! -------------------------------------------- WRITE(TP_FILE,'(a,i4,a,i3.3,a)') 'air.',iyr4,'.D',jday,'.v22.bin' WRITE(SH_FILE,'(a,i4,a,i3.3,a)') 'shum.',iyr4,'.D',jday,'.v22.bin' WRITE(PS_FILE,'(a,i4,a)') 'pres.sfc.',iyr4,'.bin' WRITE(SF_FILE,'(a,i4,a)') 'air.sig995.',iyr4,'.bin' ps_lun = get_lun() OPEN( ps_lun, FILE = DPATH//PS_FILE, & STATUS = 'OLD', & FORM = 'UNFORMATTED', & ACCESS = 'DIRECT', & RECL = sf_record_length, & IOSTAT = io_status ) IF ( io_status /= 0 ) THEN error_status = FAILURE CALL display_message( PROGRAM_NAME, & 'Error occurred opening file '//PS_FILE, & error_status ) STOP END IF sf_lun = get_lun() OPEN( sf_lun, FILE = DPATH//SF_FILE, & STATUS = 'OLD', & FORM = 'UNFORMATTED', & ACCESS = 'DIRECT', & RECL = sf_record_length, & IOSTAT = io_status ) IF ( io_status /= 0 ) THEN error_status = FAILURE CALL display_message( PROGRAM_NAME, & 'Error occurred opening file '//SF_FILE, & error_status ) STOP END IF tp_lun = get_lun() OPEN( tp_lun, FILE = DPATH//TP_FILE, & STATUS = 'OLD', & FORM = 'UNFORMATTED', & ACCESS = 'DIRECT', & RECL = tp_record_length, & IOSTAT = io_status ) IF ( io_status /= 0 ) THEN error_status = FAILURE CALL display_message( PROGRAM_NAME, & 'Error occurred opening file '//TP_FILE, & error_status ) STOP END IF sh_lun = get_lun() OPEN( sh_lun, FILE = DPATH//SH_FILE, & STATUS = 'OLD', & FORM = 'UNFORMATTED', & ACCESS = 'DIRECT', & RECL = sh_record_length, & IOSTAT = io_status ) IF ( io_status /= 0 ) THEN error_status = FAILURE CALL display_message( PROGRAM_NAME, & 'Error occurred opening file '//SH_FILE, & error_status ) STOP END IF o3_lun = get_lun() OPEN( o3_lun, FILE = O3_FILE, & STATUS = 'OLD', & IOSTAT = io_status ) IF ( io_status /= 0 ) THEN error_status = FAILURE CALL display_message( PROGRAM_NAME, & 'Error occurred opening file '//O3_FILE, & error_status ) STOP END IF !########################################################################### !--------------------------------------------------------------------------- ! Read ozone climatology !--------------------------------------------------------------------------- !########################################################################### DO IMON = 1, 12 READ(o3_lun,*) dum DO ILEV = 1, RLEVELS READ(o3_lun,'(17f9.4)') (o3_climo(IMON,ILEV,IBAND),IBAND=1,NBANDS) ENDDO ENDDO CLOSE(o3_lun) !############################################################################ !---------------------------------------------------------------------------- ! Begin Time Loop for NCEP Data for one day (6 and 18 Z only) !---------------------------------------------------------------------------- !############################################################################ v2_record_number = 1 ! DO itime = jday*4-2, jday*4, 2 time: DO itime = jday*4-4, jday*4 ! ------------------------------------ ! Define hour of NCEP obs (0,6,12,18Z) ! ------------------------------------ ihr = MOD( itime + 3 , 4 ) * 6 ! ----------------------------- ! Find month and day number !------------------------------ CALL JULIAN2 (jday,iyr4,imon,iday) ! ------------------------------------------------------------------ ! Read the temperature data (K) and specific humidity (kg/kg) arrays ! ------------------------------------------------------------------ record_number = itime READ( ps_lun, REC = record_number, IOSTAT = io_status ) ips_data IF ( io_status /= 0 ) THEN error_status = FAILURE CALL display_message( PROGRAM_NAME, & 'Error occurred reading interface pressure from '//& PS_FILE, & error_status ) STOP END IF READ( sf_lun, REC = record_number, IOSTAT = io_status ) isf_data IF ( io_status /= 0 ) THEN error_status = FAILURE CALL display_message( PROGRAM_NAME, & 'Error occurred reading interface pressure from '//& SF_FILE, & error_status ) STOP END IF READ( tp_lun, REC = v2_record_number, IOSTAT = io_status ) temperature_grid IF ( io_status /= 0 ) THEN error_status = FAILURE CALL display_message( PROGRAM_NAME, & 'Error occurred reading interface pressure from '//& TP_FILE, & error_status ) STOP END IF READ( sh_lun, REC = v2_record_number, IOSTAT = io_status ) water_vapor_grid IF ( io_status /= 0 ) THEN error_status = FAILURE CALL display_message( PROGRAM_NAME, & 'Error occurred reading interface pressure from '//& SH_FILE, & error_status ) STOP END IF v2_record_number = v2_record_number + 1 !############################################################################## !# # !# -- OPEN THE OUTPUT DATA FILES -- # !# # !############################################################################## IF(iyr4 < 2000) THEN iyr2=iyr4-1900 ELSE iyr2=iyr4-2000 ENDIF DO i_angle = 1, N_ANGLES lun_BIN = get_lun() SELECT CASE (i_angle) CASE(1) lun_BIN_10 = lun_BIN angle_name = 'A10' CASE(2) lun_BIN_09 = lun_BIN angle_name = 'A09' CASE(3) lun_BIN_08 = lun_BIN angle_name = 'A08' CASE(4) lun_BIN_07 = lun_BIN angle_name = 'A07' END SELECT WRITE(OFILE,'(a,i2.2,a,i2.2,a,i3.3,a,i2.2,a,a,a)') 'NCEP.N',satellite,'.Y', & iyr2,'.D',jday,'.',ihr,'Z.',angle_name,'.RT' out_record_length = (3 + N_LAYERS + N_LAYERS*N_CHANNELS + & N_CHANNELS*2) * 4 OPEN( lun_BIN, & FILE = DPATH//OFILE, & ACCESS = 'DIRECT', & RECL = out_record_length, & IOSTAT = io_status) IF ( io_status /= 0 ) THEN error_status = FAILURE CALL display_message( PROGRAM_NAME, & 'Error occurred opening output binary file', & error_status ) STOP END IF ENDDO ! lun_Wylie = get_lun() ! OPEN( lun_Wylie, FILE = OFILE2,IOSTAT = io_status) ! IF ( io_status /= 0 ) THEN ! error_status = FAILURE ! CALL display_message( PROGRAM_NAME, & ! 'Error occurred opening file output ascii file', & ! error_status ) ! STOP ! END IF !############################################################################## !# # !# -- BEGIN LOOP OVER EACH PROFILE -- # !# # !############################################################################## irec=1 m_lat_loop: DO j = 1, NLAT ! CALL SYSTEM_CLOCK(count,count_rate,count_max) ! write(*,*) 'Begin lat, time = ',count xlat= 92.5 - (j*2.5) !-------------------------------------------- ! Compute latitude index for Fortuin o3 data !-------------------------------------------- latidx = NINT(xlat/10.) + 9 IF(latidx < 1) latidx = 1 IF(latidx > 17) latidx = 17 !------------------------------------------------- ! Hal's ozone data - interpolate to 26 level data !------------------------------------------------- ! CALL clmozo(abs(xlat),imon,o3mix) ! CALL interpolate(40,o3mix,hal_pressure,rlevels,real(level_pressure2), & ! .true.,o3mix2) m_lon_loop: DO i = 1, NLON ! xlon= (i*2.5) - 2.5 !------------------------------------------------- ! Move profile grid data into RT input arrays !------------------------------------------------- ! temperature_grid(i,j,1) = 250. ! temperature_grid(i,j,2) = 252.5 ! temperature_grid(i,j,3) = 257.5 ! temperature_grid(i,j,4) = 260. ! dtdp=(temperature_grid(i,j,10)-temperature_grid(i,j,4))/9. ! DO ii = 5,9 ! temperature_grid(i,j,ii)=temperature_grid(i,j,4)+dtdp* & ! (level_pressure2(ii)-level_pressure2(4)) ! ENDDO level_temperature2(:) = temperature_grid(i,j,:) level_aa(:,1) = water_vapor_grid(i,j,:) level_aa(:,2) = o3_climo(imon,:,latidx) ! level_aa(:,2) = o3mix2 !----------------------------------------------------------------- ! Compute Layer-averaged quantities for Pressure, Temp, H2O and O3 !----------------------------------------------------------------- CALL level2layer( & level_pressure2, & ! level pressure data (hPa) level_temperature2, & ! level temperature data (K) level_aa, & ! (1) water vapor (g/kg), (2) ozone (ppmv) layer_pressure, & ! layer-averaged pressure layer_temperature, & ! layer-averaged temperature layer_aa) ! (1) layer-averaged water vapor, (2) ozone ! --------------------------------------------------------- ! Load RTM input profile arrays ! (Using same profile for all angles - not efficient since ! all quantities are recalculated for each angle, but, eh, ! it's a test.) ! --------------------------------------------------------- DO i_angle = 1, N_ANGLES level_p( :, i_angle ) = level_pressure2(2:RLEVELS) layer_p( :, i_angle ) = layer_pressure layer_t( :, i_angle ) = layer_temperature layer_w( :, i_angle ) = layer_aa(:,1) layer_o( :, i_angle ) = layer_aa(:,2) sfc_t ( i_angle ) = REAL(isf_data(i,j)/100.+100.) END DO !------------------------------------------------------------- ! Locate index of lower boundary using surface pressure data !------------------------------------------------------------- WHERE(REAL(ips_data(i,j)/10.) >= level_pressure2) mask = .true. ELSEWHERE mask = .false. END WHERE idx = maxloc(level_pressure2,mask) - 1 !#---------------------------------------------------------------# !# -- FORWARD MODEL -- # !#---------------------------------------------------------------# error_status = compute_rtm( & ! -- Forward inputs level_p(:idx(1),:), & ! Input, K (layers) x M (profiles) layer_p(:idx(1),:), & ! Input, K x M layer_t(:idx(1),:), & ! Input, K x M layer_w(:idx(1),:), & ! Input, K x M layer_o(:idx(1),:), & ! Input, K x M sfc_t, & ! Input, M surface_emissivity, & ! Input, L (channels) * M surface_reflectivity, & ! Input, L*M secant_view_angle, & ! Input, M secant_solar_angle, & ! Input, M n_channels_per_profile, & ! Input, M channel_index, & ! Input, L*M ! -- Forward outputs tau(:idx(1),:), & ! Output, K x L*M flux_tau(:idx(1),:), & ! Output, K x L*M solar_tau(:idx(1),:), & ! Output, K x L*M upwelling_radiance, & ! Output, L*M brightness_temperature ) ! Output, L*M IF ( error_status /= SUCCESS ) THEN CALL display_message( PROGRAM_NAME, & 'Error occured in COMPUTE_RTM', & error_status ) STOP END IF !----------------------------------------------------------- ! Output ASCII file of profile information for D. Wylie !---------------------------------------------------------- ! write( lun_Wylie, ('2f7.2,i5')) xlon,xlat,idx(1) ! write( lun_Wylie, ('2f7.2')) real(ips_data(i,j)/10.),sfc_t ! write( lun_Wylie, ('4f8.4,4f7.2')) upwelling_radiance(:), & ! brightness_temperature(:) ! do k=1,idx(1) ! write( lun_Wylie, ('f8.2,f7.2,2f11.7,4f7.4')) level_pressure2(k+1), & ! level_temperature2(k+1),level_aa(k+1,1),level_aa(k+1,2),tau(k,:) ! enddo ! write( lun_Wylie, *) ' ' DO i_angle = 1, N_ANGLES DO i_chan = 1, N_CHANNELS jdx = i_chan + (i_angle - 1) * 4 tau2(:,i_chan) = tau(:,jdx) upr2(i_chan) = upwelling_radiance(jdx) brt2(i_chan) = brightness_temperature(jdx) ENDDO SELECT CASE (i_angle) CASE(1) lun_BIN = lun_BIN_10 CASE(2) lun_BIN = lun_BIN_09 CASE(3) lun_BIN = lun_BIN_08 CASE(4) lun_BIN = lun_BIN_07 END SELECT write( lun_BIN, rec=irec ) idx(1), & REAL(sfc_t(1)), & REAL(ips_data(i,j)/10.), & REAL(level_temperature2(2:)), & REAL(tau2(:,:)), & REAL(upr2(:)), & REAL(brt2(:)) ENDDO irec=irec+1 END DO m_lon_loop write(*,*) 'Finished processing latitude = ',j ! CALL SYSTEM_CLOCK(count,count_rate,count_max) ! write(*,*) 'End lat loop, time = ',count END DO m_lat_loop ! -------------------------------- ! Output the forward model results ! -------------------------------- ! CLOSE( lun_Wylie) END DO time write(*,*) 'Time to close RT files' DO i_angle = 1, N_ANGLES SELECT CASE (i_angle) CASE(1) CLOSE(lun_BIN_10) CASE(2) CLOSE(lun_BIN_09) CASE(3) CLOSE(lun_BIN_08) CASE(4) CLOSE(lun_BIN_07) END SELECT ENDDO CLOSE( sf_lun ) CLOSE( tp_lun ) CLOSE( sh_lun ) !#----------------------------------------------------------------------------# !# -- DESTROY THE RTM -- # !#----------------------------------------------------------------------------# ! --------------------------------------- ! Deallocate the channel dependent arrays ! --------------------------------------- ! -- Forward model and ozone arrays DEALLOCATE( surface_emissivity, & ! Input, L*M surface_reflectivity, & ! Input, L*M channel_index, & ! Input, L*M n_channels_per_profile, & ! Input level_p, & ! Input layer_p, & ! Input layer_t, & ! Input layer_w, & ! Input layer_o, & ! Input sfc_t, & ! Input tau, & ! Output, K x L*M flux_tau, & ! Output, K x L*M solar_tau, & ! Output, K x L*M upwelling_radiance, & ! Output, L*M brightness_temperature, & ! Output, L*M tau2, & ! Output, K x L*M upr2, & ! Output, L*M brt2, & ! Output, L*M STAT = allocate_status ) IF ( allocate_status /= 0 ) THEN WRITE( message, '( "Error deallocating forward model channel ", & &"dependent arrays. STAT = ", i5 )' ) & allocate_status CALL display_message( PROGRAM_NAME, & TRIM( message ), & WARNING ) END IF ! --------------------------------- ! Deallocate the coefficient arrays ! --------------------------------- WRITE( *, '( /5x, "Destroying RTM...." )' ) error_status = destroy_rtm() IF ( error_status /= SUCCESS ) THEN CALL display_message( PROGRAM_NAME, & 'Error destroying RTM', & error_status ) STOP END IF END PROGRAM nceprtm_hirs_v2 !------------------------------------------------------------------------------- ! -- MODIFICATION HISTORY -- !------------------------------------------------------------------------------- ! ! $Id:$ ! ! $Date:$ ! ! $Revision:$ ! ! $State:$ ! ! $Log:$ ! ! !