! ! NAME: ! misc_darren_code ! ! PURPOSE: ! Module contains subroutines added by Darren Jackson to the NCEP RTM code. ! ! CALLING SEQUENCE: ! USE misc_darren_code ! ! MODULES ! ! type_kinds: Module to define kind types for variable declaration. ! ! error_handler: Module to define error codes and handle error conditions ! ! file_utility: Module containing generic file utility routines ! ! profile_conversion: Module containing functions and routines to allow ! for conversion of atmospheric profile concentration units ! ! CONTAINS: ! ! level2layer: Subroutine that converts level data to layer averaged data. ! ncep_channel_index: Function that returns the channel index number for a given ! satellite, instrument and channel number. ! interpolate: Linear interpolation routine ! ! CREATION HISTORY: ! Written by: Darren Jackson, NOAA/CIRES/ETL November 2001 ! MODULE misc_darren_code !--------------- ! Module usage !---------------- USE type_kinds USE error_handler USE file_utility USE profile_conversion IMPLICIT NONE CONTAINS !------------------------------------------------------------------------------ ! ! NAME: ! level2layer ! ! PURPOSE: ! Program to convert profile data located on vertical levels to layer ! averages. Layer averages are required input for the NCEP RTM code. ! ! LANGUAGE: ! Fortran-90 ! ! INPUTS: ! ! I -> number of profile points, J - number of gases ! ! level_p: pressure level data (I) (hPa) ! level_t: temperature level data (I) (K) ! level_aa: absorber array (IxJ) (water vapor g/kg, j=1) ! ! OUTPUTS: ! ! layer_p: pressure layer data (I-1) (hPa) ! layer_t: temperature layer data (I-1) (K) ! layer_aa: absorber array ((I-1)xJ) (water vapor g/kg, j=1) ! ! 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. ! ! profile_conversion: Module containing functions and routines to allow ! for conversion of atmospheric profile concentration ! units. ! ! CONTAINS: ! ! SIDE EFFECTS: ! ! RESTRICTIONS: ! ! (1) Requires that layer_ arrays be dimensioned from call statement with one ! less than the level_ arrays for the profile number dimension. ! (2) Total number of input levels can not exceed 200 without modification of ! parameter MAX_N_LAYERS. ! ! CREATION HISTORY: ! Darren Jackson, CIRES-ETL, September 2001 ! !--------------------------------------------------------------------- ! SUBROUTINE level2layer(level_p,level_t,level_absorb,layer_p,layer_t,layer_absorb) ! !--------------------------------------------------------------------------- ! Begin parameter section ! INTEGER, PARAMETER :: N_PER_LAYER = 10 REAL (fp_kind), PARAMETER :: Z0 = 0.0 CHARACTER( * ), PARAMETER :: PROGRAM_NAME = 'level2layer' ! ! End parameters section ! !---------------------------------------------------------- ! INTEGER :: N_LEVELS, N_LAYERS, N_SUBLEVELS, N_ABSORBERS INTEGER :: error_status, allocate_status CHARACTER( 128 ) :: message ! -- Inputs for level2layer subroutine REAL( fp_kind ), INTENT(IN), DIMENSION( : ) :: level_p REAL( fp_kind ), INTENT(IN), DIMENSION( : ) :: level_t REAL( fp_kind ), INTENT(IN), DIMENSION( : , : ) :: level_absorb ! -- Outputs from integrate_sublevels REAL( fp_kind ), INTENT(OUT), DIMENSION( : ) :: layer_p REAL( fp_kind ), INTENT(OUT), DIMENSION( : ) :: layer_t REAL( fp_kind ), INTENT(OUT), DIMENSION( : , : ) :: layer_absorb ! -- Outputs from create_sublevels REAL( fp_kind ), DIMENSION( : ), ALLOCATABLE :: sublevel_p, sublevel_t REAL( fp_kind ), DIMENSION( :, : ), ALLOCATABLE :: sublevel_aa REAL( fp_kind ), DIMENSION( : ), ALLOCATABLE :: sublevel_z !------------------------------------ ! Determine size of input arrays !----------------------------------- N_LEVELS = SIZE ( level_p, DIM = 1) N_ABSORBERS = SIZE ( level_absorb, DIM = 2) N_LAYERS = SIZE ( layer_p, DIM = 1) N_SUBLEVELS = ( N_LAYERS * N_PER_LAYER ) + 1 !----------------------------------------- ! Test to see if N_LAYERS exceeds Maximum !----------------------------------------- IF ( N_LAYERS /= N_LEVELS - 1 ) THEN WRITE(*,*) 'Number of layers must be one less than number of levels' WRITE(*,*) 'N_LAYERS = ',N_LAYERS,' N_LEVELS = ',N_LEVELS WRITE(*,*) 'Terminate program level2layer' STOP ENDIF !#---------------------------------------------------------------------------# ! -- ALLOCATE SPACE FOR LEVEL AND LAYER ARRAYS -- !#---------------------------------------------------------------------------# ALLOCATE( sublevel_p( N_SUBLEVELS ), & sublevel_t( N_SUBLEVELS ), & sublevel_z ( N_SUBLEVELS ), & sublevel_aa( N_SUBLEVELS, N_ABSORBERS ), & STAT = allocate_status ) IF (allocate_status /= 0) THEN WRITE( message, '( "Error allocating level2layer arrays ", & &"STAT = ",i5 )' ) allocate_status CALL display_message( PROGRAM_NAME, & TRIM( message ), & FAILURE ) STOP ENDIF !#----------------------------------------------------------------------------# !# -- INTERPOLATE THE DATA ONTO THE AIRS SUBLEVELS -- # !#----------------------------------------------------------------------------# ! ------------------------------------- ! PROFILE CONVERSION INTERPOLATION CODE ! ------------------------------------- ! -- Call to sublayer creation function ! WRITE( *, '( 10x, "Calling CREATE_SUBLEVELS...." )' ) error_status = create_sublevels( level_p, & level_t, & level_absorb, & N_PER_LAYER, & sublevel_p, & sublevel_t, & sublevel_aa ) IF ( error_status /= SUCCESS ) THEN CALL display_message( PROGRAM_NAME, & 'Error in CREATE_SUBLEVELS', & FAILURE ) STOP END IF !#----------------------------------------------------------------------------# !# -- INTEGRATE THE PROFILE DATA -- # !#----------------------------------------------------------------------------# ! WRITE( *, '( /5x, "Profile integration stage...." )' ) ! ----------------------------------- ! PROFILE CONVERSION INTEGRATION CODE ! ----------------------------------- ! -- First calculate heights. I decided to have a separate ! -- function for this as it's a useful thingy to have. ! WRITE( *, '( 10x, "Calling GEOPOTENTIAL_HEIGHT...." )' ) error_status = geopotential_height( sublevel_p, & sublevel_t, & sublevel_aa(:,1), & ! water vapor 1, & ! water vapor units are g/kg sublevel_z, & surface_height = z0, & gravity_correction = 1 ) IF ( error_status /= SUCCESS ) THEN CALL display_message( PROGRAM_NAME, & 'Error in GEOPOTENTIAL_HEIGHT', & FAILURE ) STOP END IF ! -- Now integrate the sublevels ! WRITE( *, '( 10x, "Calling INTEGRATE_SUBLEVELS...." )' ) error_status = integrate_sublevels( sublevel_z, & sublevel_p, & sublevel_t, & sublevel_aa, & n_per_layer, & 1, & ! index of h2o in sublevel_aa layer_p, & layer_t, & layer_absorb ) IF ( error_status /= SUCCESS ) THEN CALL display_message( PROGRAM_NAME, & 'Error in INTEGRATE_SUBLEVELS', & FAILURE ) STOP END IF !#---------------------------------------------------------------------------# ! -- DEALLOCATE SPACE FOR LEVEL AND LAYER ARRAYS -- !#---------------------------------------------------------------------------# DEALLOCATE( sublevel_p, & sublevel_t, & sublevel_z, & sublevel_aa, & STAT = allocate_status ) IF (allocate_status /= 0) THEN WRITE( message, '( "Error deallocating level2layer arrays ", & &"STAT = ",i5 )' ) & allocate_status CALL display_message( PROGRAM_NAME, & TRIM( message ), & FAILURE ) STOP ENDIF RETURN END SUBROUTINE level2layer ! !------------------------------------------------------------------------------ ! ! NAME: ! NCEP_CHANNEL_INDEX ! ! PURPOSE: ! Function that finds the channel index number used by the NCEP RT model. ! ! LANGUAGE: ! Fortran-90 ! ! INPUTS: ! ! satellite: NOAA polar orbit satellite number (5=TIROS-N, 6=NOAA-6, etc.) ! TYPE: INTEGER ! DIMENSION: SCALAR ! ATTRIBUTES: INTENT(IN) ! ! instrument: instrument type (1 = HIRS, 2 = MSU, 3 = AMSU) ! TYPE: INTEGER ! DIMENSION: SCALAR ! ATTRIBUTES: INTENT(IN) ! ! inst_channel: channel number for instrument (1-19 HIRS, 1-4 MSU, 1-20 AMSU) ! TYPE: INTEGER ! DIMENSION: SCALAR ! ATTRIBUTES: INTENT(IN) ! ! ! FUNCTION RESULT: ! ! Returns the NCEP RT channel index value as an SCALAR INTEGER ! ! NOTES: ! MSU channels have same index numbers for all NOAA satellites (5-14) ! AMSU channels only apply for satellite 15 and 16. ! ! CREATION HISTORY: ! Darren Jackson, CIRES-ETL, November 2001 ! !--------------------------------------------------------------------- ! FUNCTION NCEP_CHANNEL_INDEX(satellite,instrument,inst_channel) INTEGER :: satellite,instrument,inst_channel INTEGER :: begin_channel INTEGER :: ncep_channel_index ! -- Determine begin and end channel indices SELECT CASE (instrument) CASE ( 1 ) IF(satellite .ge. 5 .and. satellite .lt. 14) THEN begin_channel = (satellite - 5)*19 + 1 ! end_channel = begin_channel + 18 ELSE IF (satellite .eq. 14) THEN begin_channel = 153 ! end_channel = 171 ELSE IF (satellite .eq. 15) THEN begin_channel = 176 ! end_channel = 194 ELSE IF (satellite .eq. 16) THEN begin_channel = 215 ! end_channel = 233 ELSE WRITE(*,*) 'INVALID SATELLITE ID NUMBER FOR HIRS: SATELLITE = ',satellite STOP ENDIF CASE (2) begin_channel = 172 ! end_channel = 175 CASE (3) IF(satellite .eq. 15) THEN begin_channel = 195 ! end_channel = 214 ELSE IF (satellite .eq. 16) THEN begin_channel = 234 ! end_channel = 253 ELSE WRITE(*,*) 'INVALID SATELLITE ID NUMBER FOR AMSU: SATELLITE = ',satellite STOP ENDIF CASE DEFAULT WRITE(*,*) 'INVALID INSTRUMENT NUMBER: INSTRUMENT = ',instrument END SELECT ncep_channel_index = begin_channel + inst_channel - 1 RETURN END FUNCTION NCEP_CHANNEL_INDEX ! !------------------------------------------------------------------------------ ! ! NAME: ! INTERPOLATE ! ! PURPOSE: ! Subroutine that performs linear interpolation ! ! LANGUAGE: ! Fortran-90 ! ! INPUTS: ! n..........number of values in vector v ! v..........vector of values to be interpolated ! x..........abscissae locations of values in v. Must monotomically increase ! or decrease. ! m..........number of values for interpolated vector ! u..........new abscissae values. Length of m. Must monotomically increase ! logint.....true interpolates abscissae logarithically. ! or decrease. Must be same direction as x. ! ! OUTPUT: ! vnew.......interpolated vector values corresponding to u. ! ! CREATION HISTORY: ! Darren Jackson, CIRES-ETL, November 2001 ! !--------------------------------------------------------------------- subroutine interpolate(n,v,x,m,u,logint,vnew) integer n,m,i,j real v(n),x(n),u(m),vnew(m),slope,dir logical logint,increase increase=.false. dir=x(2)-x(1) if(abs(dir) .eq. dir) increase=.true. do i=2,n-1 dir=x(i+1)-x(i) if(increase) then if(abs(dir) .ne. dir) then write(*,*) 'X not monotomic' stop endif else if(abs(dir) .eq. dir) then write(*,*) 'X not monotomic' stop endif endif enddo increase=.false. dir=u(2)-u(1) if(abs(dir) .eq. dir) increase=.true. do i=2,m-1 dir=u(i+1)-u(i) if(increase) then if(abs(dir) .ne. dir) then write(*,*) 'U not monotomic' stop endif else if(abs(dir) .eq. dir) then write(*,*) 'U not monotomic' stop endif endif enddo if(u(1) .lt. u(2) .and. x(1) .gt. x(2)) then write(*,*) 'X and U must increase or decrease together' write(*,*) 'U increases and X decreases' stop else if(u(1) .gt. u(2) .and. x(1) .lt. x(2)) then write(*,*) 'X and U must increase or decrease together' write(*,*) 'X increases and U decreases' stop endif if(u(1) .lt. u(2)) increase=.true. do j=1,m if((increase .and. (u(j) .lt. x(1) .or. u(j) .gt. x(n))) .or. & (.not. increase .and. & (u(j) .gt. x(1) .or. u(j) .lt. x(n)))) then write(*,*) 'Can not interpolate outside abscissa range' stop endif enddo do i=1,n-1 do j=1,m if((increase .and. (u(j).ge.x(i) .and. u(j).le.x(i+1))) .or. & (.not. increase .and. & (u(j).le.x(i) .and. u(j).ge.x(i+1)))) then if(logint) then slope=log(u(j)/x(i))/log(x(i+1)/x(i)) else slope=(u(j)-x(i))/(x(i+1)-x(i)) endif vnew(j)=v(i)+slope*(v(i+1)-v(i)) endif enddo enddo return end subroutine interpolate END MODULE misc_darren_code