!------------------------------------------------------------------------------ ! ! NAME: ! ncep_ozone ! ! PURPOSE: ! Program to compute ozone profiles from ncep reanalysis temperature profiles. ! ! CATEGORY: ! ! ! 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 ! ! CONTAINS: ! None. ! ! INCLUDE FILES: ! None. ! ! EXTERNALS: ! None. ! ! COMMON BLOCKS: ! None. ! ! FILES ACCESSED: ! Input: Binary, direct access profile data file. ! Output: ASCII forward model output file. ! ! SIDE EFFECTS: ! Output file(s) are overwritten. ! ! RESTRICTIONS: ! None. ! ! CREATION HISTORY: ! Darren Jackson ETL/CIRES NOV-2001 ! Darren.L.Jackson@noaa.gov ! !-------------------------------------------------------------------------- PROGRAM ncep_ozone !#----------------------------------------------------------------------------# !# -- MODULE USAGE -- # !#----------------------------------------------------------------------------# USE type_kinds USE error_handler USE file_utility USE noaa88b USE coordinate_space !#----------------------------------------------------------------------------# !# -- TYPE DECLARATIONS -- # !#----------------------------------------------------------------------------# ! --------------------------- ! Disable all implicit typing ! --------------------------- IMPLICIT NONE ! ------------------ ! Program parameters ! ------------------ ! -- Name CHARACTER( * ), PARAMETER :: PROGRAM_NAME = 'ncep_ozone' ! -- Profile definitions INTEGER, PARAMETER :: NLON=144,NLAT=73,RLEVELS=26,TLEVELS=17,SLEVELS=8 INTEGER, PARAMETER :: NTIME=365 LOGICAL, PARAMETER :: logint=.true. INTEGER, PARAMETER :: N_LAYERS=25 REAL, PARAMETER :: MRMIN = 0.0000001 CHARACTER( * ), PARAMETER :: TP_FILE = & '/data/dlj/tovs1b/pathfinder2/ncep_data/air.1991.bin' CHARACTER( * ), PARAMETER :: SH_FILE = & '/data/dlj/tovs1b/pathfinder2/ncep_data/shum.1991.bin' CHARACTER( * ), PARAMETER :: O3_OUT_FILE = & '/data/dlj/tovs1b/pathfinder2/ncep_data/o3.1991.D200.v2.bin' CHARACTER( * ), PARAMETER :: TP_OUT_FILE = & '/data/dlj/tovs1b/pathfinder2/ncep_data/air.1991.D200.v2.bin' CHARACTER( * ), PARAMETER :: SH_OUT_FILE = & '/data/dlj/tovs1b/pathfinder2/ncep_data/shum.1991.D200.v2.bin' ! ----------------- ! Program variables ! ----------------- ! -- Error message CHARACTER( 128 ) :: message ! -- 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, orec INTEGER :: sh_in_record_length, tp_in_record_length INTEGER :: tp_out_record_length,sh_out_record_length,o3_out_record_length INTEGER :: record_number INTEGER :: tp_lun,sh_lun,tp_out_lun,sh_out_lun,o3_out_lun ! -- Profile read array INTEGER ( kind=2 ), DIMENSION( NLON, NLAT, TLEVELS ) :: itp_data INTEGER ( kind=2 ), DIMENSION( NLON, NLAT, SLEVELS ) :: ish_data REAL, DIMENSION( NLON, NLAT, RLEVELS) :: o3_grid REAL, DIMENSION( NLON, NLAT, RLEVELS) :: temperature_grid REAL, DIMENSION( NLON, NLAT, RLEVELS) :: water_vapor_grid REAL ( fp_kind ), DIMENSION( NLON, NLAT, TLEVELS ) :: sh_data = & reshape((/ (((mrmin, i=1,NLON), j=1,NLAT), k=1,TLEVELS) /), & (/ NLON,NLAT,TLEVELS /)) 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. /) REAL, DIMENSION( 2 ) :: v REAL, DIMENSION( 2 ) :: x = (/ 300.,100. /) REAL, DIMENSION( 3 ) :: u = (/ 250.,200.,150. /) REAL, DIMENSION( 3 ) :: vnew REAL( fp_kind ), DIMENSION( 15 ) :: h2o_stratosphere REAL, DIMENSION( RLEVELS ) :: level_ozone REAL, DIMENSION( 9 ) :: level_o3_temperature REAL( fp_kind ), DIMENSION( NLON, NLAT, TLEVELS ) :: level_temperature REAL( fp_kind ), DIMENSION( NLON, NLAT, TLEVELS ) :: level_water_vapor REAL( fp_kind ), DIMENSION( TLEVELS+2 ) :: level_tmp REAL( fp_kind ), DIMENSION( TLEVELS+2 ) :: level_prs ! -- Ozone profile arrays CHARACTER( * ), PARAMETER :: DATABASE_FILE = & '/home/dlj/rtmodel/nceprtm/NCEPfrtm/noaa88b.nc' INTEGER :: n_levels,n_profiles,profile_number REAL( fp_kind ), ALLOCATABLE, DIMENSION( : ) :: ozone_pressure REAL( fp_kind ), ALLOCATABLE, DIMENSION( :, : ) :: ozone_temperature REAL( fp_kind ), ALLOCATABLE, DIMENSION( :, : ) :: ozone REAL( fp_kind ), ALLOCATABLE, DIMENSION( : ) :: level_tmp_airs REAL, ALLOCATABLE, DIMENSION( : ) :: level_temperature_airs ! ---------- ! Intrinsics ! ---------- INTRINSIC COS, & MIN, MAX, & REAL, & TRIM, & MOD ! -- Assume stratospheric h2o = 3 ppmv (1.87e-6 kg/kg) from 100 mb to 0.1 mb h2o_stratosphere(:) = 0.00000187 !############################################################################## !############################################################################## !############################################################################## !# # !# -- SETUP DATA AND ARRAYS FOR CALCULATIONS -- # !# # !############################################################################## !############################################################################## !############################################################################## !############################################################################## !############################################################################## !############################################################################## !# # !# -- OPEN THE INPUT PROFILE DATA FILE -- # !# # !############################################################################## !############################################################################## !############################################################################## ! ------------------------------ ! Set the record length in bytes ! ------------------------------ tp_in_record_length = NLAT * NLON * TLEVELS * n_bytes_for_Short_kind sh_in_record_length = NLAT * NLON * SLEVELS * n_bytes_for_Short_kind o3_out_record_length = NLAT * NLON * RLEVELS * 4 tp_out_record_length = NLAT * NLON * RLEVELS * 4 sh_out_record_length = NLAT * NLON * RLEVELS * 4 ! -------------------------------------------- ! Open input temperature and specific humidity files ! -------------------------------------------- tp_lun = get_lun() OPEN( tp_lun, FILE = TP_FILE, & STATUS = 'OLD', & FORM = 'UNFORMATTED', & ACCESS = 'DIRECT', & RECL = tp_in_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 = SH_FILE, & STATUS = 'OLD', & FORM = 'UNFORMATTED', & ACCESS = 'DIRECT', & RECL = sh_in_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 !------------------------------------ ! Open output data files !------------------------------------ o3_out_lun = get_lun() OPEN( o3_out_lun, FILE = O3_OUT_FILE, & FORM = 'UNFORMATTED', & ACCESS = 'DIRECT', & RECL = o3_out_record_length, & IOSTAT = io_status ) IF ( io_status /= 0 ) THEN error_status = FAILURE CALL display_message( PROGRAM_NAME, & 'Error occurred opening file '//O3_OUT_FILE, & error_status ) STOP END IF tp_out_lun = get_lun() OPEN( tp_out_lun, FILE = TP_OUT_FILE, & FORM = 'UNFORMATTED', & ACCESS = 'DIRECT', & RECL = tp_out_record_length, & IOSTAT = io_status ) IF ( io_status /= 0 ) THEN error_status = FAILURE CALL display_message( PROGRAM_NAME, & 'Error occurred opening file '//TP_OUT_FILE, & error_status ) STOP END IF sh_out_lun = get_lun() OPEN( sh_out_lun, FILE = SH_OUT_FILE, & FORM = 'UNFORMATTED', & ACCESS = 'DIRECT', & RECL = sh_out_record_length, & IOSTAT = io_status ) IF ( io_status /= 0 ) THEN error_status = FAILURE CALL display_message( PROGRAM_NAME, & 'Error occurred opening file '//SH_OUT_FILE, & error_status ) STOP END IF !############################################################################ !############################################################################ !############################################################################ ! ! READ OZONE DATABASE AND ESTABLISH OZONE PROFILES ! !############################################################################ !############################################################################ !############################################################################ !#----------------------------------------------------------------------------# !# -- INQUIRE THE DATABASE -- # !#----------------------------------------------------------------------------# WRITE( *, '( /5x, "Inquiring the database ", a, "..." )' ) DATABASE_FILE error_status = inquire_noaa88b( database_file, & n_profiles = n_profiles, & n_levels = n_levels ) IF ( error_status /= SUCCESS ) THEN CALL display_message( PROGRAM_NAME, & 'Error inquiring database file '//DATABASE_FILE, & FAILURE ) STOP END IF ! -- Output results WRITE( *, '( 10x, "Number of profiles : ", i6, & &/10x, "Number of atmospheric layers : ", i6 )' ) & n_profiles, n_levels !#----------------------------------------------------------------------------# !# -- ALLOCATE THE REQUIRED OZONE DATA ARRAYS -- # !#----------------------------------------------------------------------------# WRITE( *, '( /5x, "Allocating the data read arrays..." )' ) ALLOCATE( ozone_pressure( n_levels ), & ozone_temperature( n_levels, n_profiles ), & ozone( n_levels, n_profiles ), & level_temperature_airs( n_levels ), & level_tmp_airs( n_levels ), & STAT = allocate_status) IF ( allocate_status /= 0 ) THEN WRITE( message, '( "Error allocating data arrays. STAT = ", i5 )' ) & allocate_status CALL display_message( PROGRAM_NAME, & TRIM( message ), & FAILURE ) STOP END IF !#----------------------------------------------------------------------------# !# -- READ THE OZONE PROFILE DATABASE -- # !#----------------------------------------------------------------------------# WRITE( *, '( /5x, "Reading the database ", a, "..." )' ) DATABASE_FILE error_status = read_noaa88b( database_file, & pressure = ozone_pressure, & temperature = ozone_temperature, & ozone = ozone ) IF ( error_status /= SUCCESS ) THEN CALL display_message( PROGRAM_NAME, & 'Error reading database file '//DATABASE_FILE, & FAILURE ) STOP END IF !#----------------------------------------------------------------------------# !# -- CREATE OZONE COORDINATE SPACE -- # !#----------------------------------------------------------------------------# WRITE( *, '( /5x, "Initializing the coordinate space..." )' ) error_status = initialize_coordinate_space( TRANSPOSE( ozone_temperature ), & debug = 1 ) IF ( error_status /= SUCCESS ) THEN CALL display_message( PROGRAM_NAME, & 'Error initialising coordinate space.', & FAILURE ) error_status = destroy_coordinate_space() STOP END IF !############################################################################ !---------------------------------------------------------------------------- ! Begin Time Loop for NCEP Data !---------------------------------------------------------------------------- !############################################################################ orec=1 ! DO itime = 1, NTIME ! DO itime = 200, 200 DO itime = 801, 804 ! ------------------------------------------------------------------ ! Read the temperature data (K) and specific humidity (kg/kg) arrays ! ------------------------------------------------------------------ record_number = itime READ( tp_lun, REC = record_number, IOSTAT = io_status ) itp_data 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 = record_number, IOSTAT = io_status ) ish_data 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 !--------------------------------------------------------- ! Extrapolate water vapor profile in upper troposphere !--------------------------------------------------------- sh_data(:,:,1:8)=ish_data(:,:,1:8)*0.000001+0.025 DO i=1,NLON DO j=1,NLAT v(1)=sh_data(i,j,8) v(2)=real(h2o_stratosphere(1)) call interpolate(2,v,x,3,u,logint,vnew) sh_data(i,j,9:11)=vnew sh_data(i,j,12:17)=h2o_stratosphere(10:15) ENDDO ENDDO !------------------------------------------------------------------------------ ! Convert specific humidity to mixing ratio (g/kg); rotate direction of profile !------------------------------------------------------------------------------ DO i=1,NLON DO j=1,NLAT DO K=1,TLEVELS level_water_vapor(i,j,TLEVELS+1-k)=sh_data(i,j,k)/(1.-sh_data(i,j,k))*1000. ENDDO ENDDO ENDDO !------------------------------------------------------------------------------- ! Scale and offset the input temperature field (K); rotate direction of profile !------------------------------------------------------------------------------- DO i=1,NLON DO j=1,NLAT DO K=1,TLEVELS level_temperature(i,j,TLEVELS+1-k)=itp_data(i,j,k)*0.01+100. ENDDO ENDDO ENDDO !############################################################################## !# # !# -- BEGIN LOOP OVER EACH PROFILE -- # !# # !############################################################################## m_lat_loop: DO j = 1, NLAT m_lon_loop: DO i = 1, NLON !-------------------------------------------------------------------- ! Interpolate the NCEP temperature data onto the AIRS pressure levels !-------------------------------------------------------------------- level_prs(2:tlevels+1)=level_pressure level_prs(1)=0.015 level_prs(tlevels+2)=1100.0 level_tmp(2:tlevels+1)=level_temperature(i,j,:) level_tmp(1)=level_temperature(i,j,1) level_tmp(tlevels+2)=level_temperature(i,j,tlevels) call interpolate(tlevels+2,real(level_tmp),real(level_prs), & n_levels,real(ozone_pressure),logint,level_temperature_airs) ! !---------------------------------------------------------- ! Find profile for Ozone using closest temperature profile !---------------------------------------------------------- ! level_tmp_airs=level_temperature_airs error_status = find_closest_vector( level_tmp_airs, & profile_number ) IF ( error_status /= SUCCESS ) THEN CALL display_message( PROGRAM_NAME, & 'Error matching temperature vector.', & error_status ) error_status = destroy_coordinate_space() STOP END IF !------------------------------------------------------------------------ ! Interpolate "best" ozone profile to RT vertical profile !------------------------------------------------------------------------ CALL interpolate(100,real(ozone(:,profile_number)),real(ozone_pressure), & rlevels,real(level_pressure2),logint,level_ozone) !-------------------------------------------------------------------------- ! Interpolate AIRS temperature profile above 10 mb onto RT vertical profile !--------------------------------------------------------------------------- CALL interpolate(19,real(ozone_temperature(:19,profile_number)), & real(ozone_pressure(:19)), & 9,real(level_pressure2(:9)),logint,level_o3_temperature) !-------------------------------------------------------------- ! Put temperature, h2o, o3 profiles into grid arrays !-------------------------------------------------------------- temperature_grid(i,j,1:9) =level_o3_temperature temperature_grid(i,j,10:) =level_temperature(i,j,:) water_vapor_grid(i,j,:9) = h2o_stratosphere(:9)*1000. water_vapor_grid(i,j,10:) = level_water_vapor(i,j,:) o3_grid(i,j,:) = level_ozone END DO m_lon_loop END DO m_lat_loop ! -------------------------------------------- ! Output the forward model results for one day ! -------------------------------------------- WRITE( tp_out_lun, rec=orec ) temperature_grid WRITE( sh_out_lun, rec=orec ) water_vapor_grid WRITE( o3_out_lun, rec=orec ) o3_grid orec=orec+1 WRITE(*,*) 'Completed processing time index',itime END DO CLOSE( tp_lun ) CLOSE( sh_lun ) CLOSE( tp_out_lun ) CLOSE( sh_out_lun ) CLOSE( o3_out_lun ) DEALLOCATE( ozone_pressure, & ozone_temperature, & ozone, & level_temperature_airs, & level_tmp_airs, & 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 END PROGRAM ncep_ozone !------------------------------------------------------------------------------- ! -- MODIFICATION HISTORY -- !------------------------------------------------------------------------------- ! ! $Id:$ ! ! $Date:$ ! ! $Revision:$ ! ! $State:$ ! ! $Log:$ ! ! !