!------------------------------------------------------------------------------ ! ! NAME: ! ncep_26level ! ! PURPOSE: ! Program to interpolate temperature and humidity profiles onto 26 level ! profile. ! ! ! 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 ! ! FILES ACCESSED: ! Input: Binary, direct access profile data file. ! upper stratospheric temperature climatology file ! Output: Binary, direct access profile data file on 26 levels. ! ! SIDE EFFECTS: ! Output file(s) are overwritten. ! ! RESTRICTIONS: ! None. ! ! CREATION HISTORY: ! Darren Jackson ETL/CIRES December-2001 ! Darren.L.Jackson@noaa.gov ! !-------------------------------------------------------------------------- PROGRAM ncep_26level !#----------------------------------------------------------------------------# !# -- MODULE USAGE -- # !#----------------------------------------------------------------------------# USE type_kinds USE error_handler USE file_utility USE dljf90lib !#----------------------------------------------------------------------------# !# -- TYPE DECLARATIONS -- # !#----------------------------------------------------------------------------# ! --------------------------- ! Disable all implicit typing ! --------------------------- IMPLICIT NONE ! ------------------ ! Program parameters ! ------------------ ! -- Name CHARACTER( * ), PARAMETER :: PROGRAM_NAME = 'ncep_26level' ! -- Profile definitions INTEGER, PARAMETER :: NLON=144,NLAT=73,RLEVELS=26,TLEVELS=17,SLEVELS=8 LOGICAL, PARAMETER :: logint=.true. REAL, PARAMETER :: MRMIN = 0.0000001 CHARACTER (*), PARAMETER :: DPATH = '/data/dlj/tovs1b/pathfinder2/ncep_data/' ! ----------------- ! 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 INTEGER :: record_number INTEGER :: tp_lun,sh_lun,tp_out_lun,sh_out_lun,in_lun INTEGER :: iyr4,ntime,yyddd,iyr,imon,jday ! -- 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) :: temperature_grid REAL, DIMENSION( NLON, NLAT, RLEVELS) :: water_vapor_grid REAL, DIMENSION( 40 ) :: T_Hal REAL, DIMENSION( 6 ) :: tmp 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( fp_kind ), DIMENSION( NLON, NLAT, TLEVELS ) :: level_temperature REAL( fp_kind ), DIMENSION( NLON, NLAT, TLEVELS ) :: level_water_vapor REAL :: xlat CHARACTER( 80 ) :: TP_FILE CHARACTER( 80 ) :: SH_FILE CHARACTER( 80 ) :: O3_OUT_FILE CHARACTER( 80 ) :: TP_OUT_FILE CHARACTER( 80 ) :: SH_OUT_FILE ! ---------- ! 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 ! -------------------------------------- ! Input information !--------------------------------------- in_lun = get_lun() OPEN(in_lun,FILE='./ncep_26level.input') READ(in_lun,*) iyr4 READ(in_lun,*) jday CLOSE(in_lun) !############################################################################## !############################################################################## !############################################################################## !# # !# -- 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 tp_out_record_length = NLAT * NLON * RLEVELS * 4 sh_out_record_length = NLAT * NLON * RLEVELS * 4 ! -------------------------------------------------- ! Open input temperature and specific humidity files ! -------------------------------------------------- WRITE(TP_FILE,'(a,a,i4,a)') dpath,'air.',iyr4,'.bin' WRITE(SH_FILE,'(a,a,i4,a)') dpath,'shum.',iyr4,'.bin' 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 !------------------------------------ WRITE(TP_OUT_FILE,'(a,a,i4,a,i3.3,a)') dpath,'air.',iyr4,'.D',jday,'.v22.bin' WRITE(SH_OUT_FILE,'(a,a,i4,a,i3.3,a)') dpath,'shum.',iyr4,'.D',jday,'.v22.bin' 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 !############################################################################ !---------------------------------------------------------------------------- ! Begin Time Loop for NCEP Data !---------------------------------------------------------------------------- !############################################################################ orec=1 DO itime = jday*4-2, jday*4, 2 IF(iyr4 < 2000) THEN yyddd=INT(itime/4)+1+(iyr4-1900)*1000 ELSE yyddd=INT(itime/4)+1+(iyr4-2000)*1000 ENDIF ! ------------------------------------------------------------------ ! 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 IF(sh_data(i,j,k) <= 0.) sh_data(i,j,k) = mrmin 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 xlat= 92.5 - (j*2.5) !----------------------------------------------------------------------- ! Retrieve Hal's upper stratospheric temperature data !----------------------------------------------------------------------- CALL CLIMATM(yyddd,xlat,T_Hal) !----------------------------------------------------------------------- ! Retrieve CPC upper stratospheric temperature data !----------------------------------------------------------------------- ! CALL CPC_Temp(xlat,imon,T_CPC) ! write(*,*) T_CPC ! stop m_lon_loop: DO i = 1, NLON !-------------------------------------------------------------------------- ! Interpolate CPC 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 Hal's upper stratospheric temperature into temp profile !----------------------------------------------------------------- temperature_grid(i,j,1:4) = T_Hal(1:4) temperature_grid(i,j,5:9) = T_Hal(6:10) !-------------------------------------------------------------- ! Create NCEP Reanalysis temperature and h2o grid arrays !-------------------------------------------------------------- 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,:) !-------------------------------------------------------------------- ! Smooth the temperature profile between upper and lower stratosphere !-------------------------------------------------------------------- CALL smooth(temperature_grid(i,j,7:12), 6, 3, .false., tmp) temperature_grid(i,j,7:12) = tmp(:) 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 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 ) END PROGRAM ncep_26level !------------------------------------------------------------------------------- ! -- MODIFICATION HISTORY -- !------------------------------------------------------------------------------- ! ! $Id:$ ! ! $Date:$ ! ! $Revision:$ ! ! $State:$ ! ! $Log:$ ! ! !