! $Id$ ! module avhrr_read_mdl USE ISO_C_BINDING IMPLICIT NONE private public :: get_avhrr_data enum , bind(C) enumerator :: ETchannel_inavlid = 0 enumerator :: ETchannel_FIRST = 1 enumerator :: ETchannel_FIRST_VIS = 1 enumerator :: ETchannel_1 = 1 enumerator :: ETchannel_2 = 2 enumerator :: ETchannel_3a = 3 enumerator :: ETchannel_LAST_VIS = 3 enumerator :: ETchannel_FIRST_NIR = 4 enumerator :: ETchannel_3b = 4 enumerator :: ETchannel_4 = 5 enumerator :: ETchannel_5 = 6 enumerator :: ETchannel_LAST_NIR = 6 enumerator :: ETchannel_last = 6 END enum INTEGER, parameter , public:: NUM_CHAN = 6 ! - input configure structure TYPE , public :: avhrr_data_config LOGICAL :: chan_on_rfl ( NUM_CHAN ) LOGICAL :: chan_on_rad ( NUM_CHAN ) LOGICAL :: chan_on_bt ( NUM_CHAN ) CHARACTER ( len = 256 ) :: file INTEGER, dimension( 2 ) :: offset INTEGER, dimension( 2 ) :: count CHARACTER ( len =255 ) :: dir_1b CHARACTER ( len =255 ) :: abhrr_file_base CHARACTER (len = 200 ) :: Ancil_Data_Dir LOGICAL :: use_intern_calib_refl LOGICAL :: use_intern_calib_thermal END TYPE avhrr_data_config ! -output structures and substructures TYPE :: band_str REAL,dimension (:,:) , allocatable :: ref REAL,dimension (:,:) , allocatable :: rad REAL,dimension (:,:) , allocatable :: bt END TYPE band_str TYPE :: geo_str REAL , allocatable :: solzen (:,:) REAL , dimension (:,:) , allocatable :: satzen REAL , dimension (:,:) , allocatable :: solaz REAL , dimension (:,:) , allocatable :: sataz REAL , dimension (:,:) , allocatable :: relaz REAL , dimension (:,:) , allocatable :: lat REAL , dimension (:,:) , allocatable :: lon REAL , dimension (:) , allocatable :: scan_time INTEGER, dimension (:) , allocatable :: ascEND END TYPE geo_str TYPE, public :: avhrr_data_out LOGICAL :: is_allocated LOGICAL :: is_set = .false. TYPE (band_str), dimension(NUM_CHAN) :: band TYPE ( geo_str) :: geo contains procedure :: alloc_it => alloc_out procedure :: dealloc_it => dealloc_out END TYPE avhrr_data_out ! - enumerators TYPE :: ET_avhrr_t_TYPE INTEGER:: MISSING = 0 INTEGER:: FIRST = 1 INTEGER:: KLM = 1 INTEGER:: AAPP = 2 INTEGER:: PRE_KLM =3 END TYPE ET_avhrr_t_TYPE TYPE ( ET_avhrr_t_TYPE ) :: ET_avhrr_t ! some params INTEGER, parameter :: INT1 = selected_int_kind(R=2) INTEGER, parameter :: INT2 = selected_int_kind(R=4) INTEGER, parameter :: INT4 = selected_int_kind(R=9) INTEGER, parameter, public:: REAL4 = selected_REAL_kind(6,37) ! local informations needed TYPE avhrr_meta_TYPE LOGICAL :: is_GAC_not_LAC INTEGER:: l1b_rec_length INTEGER:: nrec_header INTEGER:: nword_clavr_start INTEGER:: nword_clavr INTEGER:: num_pix INTEGER:: avhrr_TYPE INTEGER:: sensor_id INTEGER( int2 ) :: ver_1b INTEGER:: num_data_bytes INTEGER:: num_digtel_bytes INTEGER:: num_house_bytes INTEGER:: num_clavr_bytes INTEGER:: first_data_byte END TYPE avhrr_meta_TYPE TYPE avhrr_line_TYPE INTEGER(kind = int2):: scan_year INTEGER(kind = int2):: scan_day INTEGER(kind=int4) :: scan_time LOGICAL :: is_ascEND LOGICAL :: ch3a_on INTEGER:: scan_number REAL :: vis_slope_1(3) REAL :: vis_intercept_1(3) REAL :: vis_slope_2(3) REAL :: vis_intercept_2(3) REAL :: vis_intersection(3) REAL :: ir_coef_1_1b(4:6) REAL :: ir_coef_2_1b(4:6) REAL :: ir_coef_3_1b(4:6) REAL :: solzen_anchor (51) REAL :: satzen_anchor (51) REAL :: relaz_anchor (51) INTEGER(kind=int2), dimension(5,10) :: Space_Count_Temp REAL (kind = REAL4) :: lat_anchor(51) REAL (kind = REAL4) :: lon_anchor(51) INTEGER( kind = int2) , allocatable :: chan_counts (:,:) INTEGER( kind = int1) , allocatable :: cld_mask_aux (:) INTEGER( kind = int1 ) :: clavr_status END TYPE avhrr_line_TYPE LOGICAL, parameter :: IS_SIGNED = .true. LOGICAL, parameter :: IS_UNSIGNED = .false. TYPE instr_const_TYPE CHARACTER (len =256) :: file CHARACTER ( len = 7 ) :: sat_name REAL(kind=REAL4) :: Gain_Low_0 (3) REAL(kind=REAL4) :: Gain_High_0 (3) REAL(kind=REAL4) :: Degrad_Low_1 (3) REAL(kind=REAL4) :: Degrad_High_1 (3) REAL(kind=REAL4) :: Degrad_Low_2 (3) REAL(kind=REAL4) :: Degrad_High_2 (3) REAL(kind=REAL4) :: Dark_Count(3) REAL(kind=REAL4) :: switch_Count(3) REAL (kind=REAL4) :: c1 = 1.191062e-5, & c2 = 1.4387863 REAL (kind=REAL4) :: A1_20,A2_20,Nu_20 REAL (kind=REAL4) :: A1_31,A2_31,Nu_31 REAL (kind=REAL4) :: A1_32,A2_32,Nu_32 REAL(kind=REAL4) :: sun_Earth_distance REAL(kind=REAL4) :: launch_date REAL(kind=REAL4) :: solar_ch20 REAL(kind=REAL4) :: solar_ch20_nu REAL(kind=REAL4) :: Ew_ch20 ! equivalent width REAL (kind=REAL4) :: & space_Rad_3b,space_Rad_4,space_Rad_5 REAL (kind=REAL4) :: b0 (4:6) , b1 (4:6) , b2(4:6) REAL(kind=REAL4), dimension(0:4,4):: Prt_Coef REAL(kind=REAL4), dimension(4):: Prt_Weight REAL(kind=REAL4) :: B1_day_mask,B2_day_mask, & B3_day_mask,B4_day_mask contains procedure :: READ_AVHRR_INSTR_CONSTANTS END TYPE instr_const_TYPE TYPE ( instr_const_TYPE ) :: instr_cons REAL, parameter :: Missing_Value_REAL4 = -999. contains !======================================================================================= ! ! The main program to read and prepare avhrr data ! ! !======================================================================================== SUBROUTINE get_avhrr_data ( config , out ) use file_tools TYPE ( avhrr_data_config ), INTENT (in out) :: config TYPE ( avhrr_data_out ), INTENT ( out ) :: out CHARACTER ( len =255) :: file TYPE ( avhrr_meta_TYPE) :: meta INTEGER(kind=int4) :: nx_start , nx_END , ny_start , ny_END INTEGER( kind = int4) :: dim_seg(2) INTEGER:: line_idx INTEGER( kind=int1) , dimension(:), allocatable ::segment_buffer INTEGER( kind=int1) , dimension(:), allocatable ::line_buffer INTEGER:: word_start INTEGER:: word_END INTEGER:: bytes_per_word INTEGER:: number_of_words INTEGER:: number_of_words_read TYPE ( avhrr_line_TYPE) :: data_line REAL :: year_since_launch LOGICAL :: Ref_Cal_1b = .false. LOGICAL :: therm_cal_1b= .true. INTEGER(kind=int2), dimension(:,:,:), allocatable :: Chan_Counts_segm CHARACTER ( len =256) :: Instr_Const_file REAL (kind= REAL4 ), allocatable, dimension(:,:) :: rad_line INTEGER:: i ! ======== executable ================================== ! check IF compressed and uncompress IF needed call uncompress_file ( config % file , file & , config % dir_1b & , "./temp/") ! - 1.get avhrr meta data call get_avhrr_meta ( trim(file) , meta ) ! - dimensions nx_start is always 1 nx_start = 1 ny_start = config % offset(2) nx_END = meta % num_pix ny_END = ny_start + config % count(2) - 1 dim_seg = [ meta % num_pix , ny_END - ny_start+1 ] !- allocate locals allocate ( segment_buffer ( meta % l1b_rec_length * dim_seg ( 2) )) allocate ( Chan_Counts_segm ( num_chan, dim_seg(1), dim_seg(2) ) ) allocate ( data_line % chan_counts ( num_chan , meta % num_pix) ) allocate ( data_line % cld_mask_aux ( meta % num_pix) ) ! - allocate output call out % alloc_it ( dim_seg(1), dim_seg(2) ) ! - populate instrument constant structure Instr_Const_file = trim ( config % Ancil_Data_Dir)//"avhrr_data/"//"avhrr_16_instr.dat" call instr_cons % READ_AVHRR_INSTR_CONSTANTS( trim(Instr_Const_file) ) ! - read allthe data is segment_buffer bytes_per_word = 1 word_start = meta % l1b_rec_length * ( ny_start - 1 ) Number_Of_Words = meta % l1b_rec_length * dim_seg ( 2 ) call mreadf_int (trim(file)//CHAR(0), word_start , bytes_per_word, & Number_of_Words, Number_Of_Words_Read, segment_buffer ) ! - loop over lines loop_line : do line_idx = ny_start, ny_END ! - sub-buffer for this line word_start = ( line_idx - 1 ) * meta % l1b_rec_length + 1 word_END = word_start + meta % l1b_rec_length - 1 allocate(line_buffer ( meta % l1b_rec_length )) line_buffer = segment_buffer ( word_start : word_END ) ! - unpack the data call unpack_avhrr_data_record_klm ( line_buffer , meta , data_line) !- we need the counts as a segemetn , could be also done linewise Chan_Counts_segm ( :,line_idx , : ) = data_line % chan_counts ! compute lat/lon IF ( maxval ( abs ( data_line % Lat_anchor ) ) > 85 ) then call GNOMIC_ANCHOR_INTERP ( data_line % lon_anchor, data_line % lat_anchor , meta & & , out % geo % lon (:,line_idx) , out % geo % lat (:,line_idx) ) ELSE call lagrangian_anchor_interp ( data_line % lat_anchor , meta , out % geo % lat (:,line_idx) ) call lagrangian_anchor_interp ( data_line % lon_anchor , meta , out % geo % lon (:,line_idx) ) END IF ! - for calibration reasons we need this year_since_launch = data_line % scan_year + data_line % scan_day/365. - instr_cons % launch_date allocate ( rad_line ( 4:6 , meta % num_pix) ) call thermal_calibration ( & config % use_intern_calib_thermal & , data_line & , meta & , config % chan_on_rfl & !-TODO have to make this line specIFic , year_since_launch & , rad_line ) do i = 4, 6 out % band(i ) % rad ( :, line_idx ) = rad_line (i,:) END do !IF (line_idx == 200) print*, rad_line(5,:) deallocate ( rad_line ) deallocate( line_buffer ) END do loop_line !-todo ! - reflection can be done on segement call reflectance_calibration (Chan_Counts_segm, & config % use_intern_calib_refl & , config % use_intern_calib_refl & , data_line & , year_since_launch & , out % band ( 1:3) ) out % is_set = .true. END SUBROUTINE get_avhrr_data ! determine avhrr TYPE ! should do any avhtt_TYPE configuration needed in this ! module ! ! !======================================================================================= ! ! ! ! !======================================================================================== SUBROUTINE get_avhrr_meta ( file , meta ) CHARACTER ( len = *) , INTENT (in) :: file TYPE ( avhrr_meta_TYPE) , INTENT(out) :: meta INTEGER( kind = int1) , dimension(100):: header_Buffer_temp INTEGER:: word_start INTEGER:: bytes_per_word INTEGER:: Number_Of_Words INTEGER:: Number_Of_Words_Read INTEGER(kind=int2) :: data_TYPE ! --- executable --- meta % avhrr_TYPE = ET_avhrr_t % MISSING bytes_per_word = 1 word_start = 0 Number_Of_Words = 100 print*,file//CHAR(0),'==' call mreadf_int(file//CHAR(0),word_start,bytes_per_word, & Number_Of_Words,Number_Of_Words_Read,header_Buffer_temp) print*,'m' ! - KLM IF ( header_buffer_temp(1) >= 65 ) then meta % avhrr_TYPE = ET_avhrr_t % KLM call b2_to_ui2 (header_Buffer_temp(77:78), data_TYPE) call b2_to_ui2 ( header_Buffer_temp(73:74) , meta % sensor_id ) call b2_to_ui2 ( header_buffer_temp(5:6) , meta % ver_1b ) meta % avhrr_TYPE = ET_avhrr_t % KLM ! - check wrong version_l1b IF ( meta % ver_1b > 10 .or. meta % ver_1b < 1 ) then ! swap the bytes we assume AAPP data call b2_to_ui2 ( header_buffer_temp(6:5:-1) , meta % ver_1b ) call b2_to_ui2 ( header_Buffer_temp(74:73:-1) , meta % sensor_id ) call b2_to_ui2 (header_Buffer_temp(78:77:-1), data_TYPE) meta % avhrr_TYPE = ET_avhrr_t % AAPP END IF ! fix noaa-15 bug IF ( meta % ver_1b ==1 .and. meta % sensor_id == 4 ) meta % ver_1b = 2 ELSE ! - non-klm data_TYPE = ishft ( header_buffer_temp (2) , -4 ) meta % ver_1b = 1 call b2_to_ui2 ( header_Buffer_temp(73:74) , meta % sensor_id ) meta % avhrr_TYPE = ET_avhrr_t % PRE_KLM END IF meta % is_gac_not_lac = .false. IF ( data_TYPE == 2) meta % is_gac_not_lac = .true. print*,'version level1b: ',meta % ver_1b ! populate sensor meta IF ( meta % avhrr_TYPE == ET_avhrr_t % PRE_KLM ) then IF ( meta % is_gac_not_lac ) then meta % num_pix = 409 meta % l1b_rec_length = 3220 meta % nrec_header = 2 meta % nword_clavr_start = 0 meta % nword_clavr = 0 meta % Num_data_bytes = 2728 ELSE meta % num_pix = 2048 meta % l1b_rec_length = 14800 meta % nrec_header = 1 meta % nword_clavr_start = 0 meta % nword_clavr = 0 meta % Num_data_bytes = 6952 + 6704 ! -- What ? END IF END IF IF ( meta % avhrr_TYPE == ET_avhrr_t % KLM ) then IF ( meta % is_gac_not_lac ) then meta % num_pix = 409 meta % l1b_rec_length = 4608 meta % nrec_header = 1 meta % nword_clavr_start = 4049 meta % nword_clavr = 2 * 52 meta % Num_data_bytes = 4000 - 1265 + 1 meta % Num_digtel_bytes = 4016 - 4001 + 1 meta % Num_house_bytes = 4048 - 4017 + 1 meta % Num_clavr_bytes = 4160 - 4049 + 1 ELSE meta % num_pix = 2048 meta % l1b_rec_length = 15872 meta % nrec_header = 1 meta % nword_clavr_start = 14977 meta % nword_clavr = 2 * 256 meta % Num_data_bytes = 14928 - 1265 meta % Num_digtel_bytes = 14944 - 14929 meta % Num_house_bytes = 14976 - 14971 meta % Num_clavr_bytes = 15496 - 14977 END IF END IF IF ( meta % avhrr_TYPE == ET_avhrr_t % AAPP ) then meta % num_pix = 2048 meta % l1b_rec_length = 22016 meta % nrec_header = 1 meta % nword_clavr_start = 0 meta % nword_clavr = 0 meta % first_data_byte = 1265 meta % Num_digtel_bytes = 14944 - 14929 meta % Num_house_bytes = 14976 - 14971 meta % Num_clavr_bytes = 0 END IF ! - TODO Sensor name wmo number etc... END SUBROUTINE get_avhrr_meta !======================================================================================= ! ! ! ! !======================================================================================== SUBROUTINE UNPACK_AVHRR_DATA_RECORD_KLM( rec , meta, out) INTEGER( kind = int1 ), INTENT(in), dimension(:) :: rec TYPE (avhrr_meta_TYPE) , INTENT (in) :: meta TYPE ( avhrr_line_TYPE ), INTENT(inout) :: out INTEGER(kind=int4):: i4word INTEGER(kind=int2):: i2word INTEGER(kind=int4):: start_pixel INTEGER(kind=int4):: END_pixel INTEGER:: ichan INTEGER:: ipix INTEGER:: i INTEGER:: j INTEGER:: ibyte INTEGER:: iword INTEGER:: start_word INTEGER:: END_word INTEGER:: word_clavr_start INTEGER(kind=int1):: onebyte INTEGER:: num_loc LOGICAL :: valid_therm_cal !------------------------------------------------------------------------------- ! set parameters for Gac or lac processing !------------------------------------------------------------------------------- Valid_Therm_Cal = .true. !---------------------------------------------------------------- ! unpack scan line number !---------------------------------------------------------------- call b4_to_ui4 ( rec(1:2) , out % scan_number ) !---------------------------------------------------------------------- ! unpack time_temp code !---------------------------------------------------------------------- call b2_to_ui2 (rec(3:4) , out % Scan_Year ) call b2_to_ui2 (rec(5:6) , out % Scan_day ) call b4_to_ui4 (rec(9:12), out % scan_time ) ! IF (Scan_Time(iline) > Biggest_I4) Scan_Time(iline) = Missing_Value_Int4 ! IF (Scan_Time(iline) < -1*Biggest_I4) Scan_Time(iline) = Missing_Value_Int4 !---------------------------------------------------------------------- ! unpack quality indicators - not done fully !---------------------------------------------------------------------- i2word = rec(13) IF (i2word < 0) then i2word = i2word + 256 END IF out % is_ascend = ishft(i2word, -7) == 1 i2word = rec(14) IF (i2word < 0) then i2word = i2word + 256 END IF out % Ch3a_On = ishft(ishft(i2word, 6),-6) == 1 !-TODO Ch6_On_Pixel_Mask(:,iline) = Ch3a_On(iline) !--------------------------------------------------------------------- ! unpack space view !--------------------------------------------------------------------- do j = 1,10 call b2_to_si2 ( rec ( 1161+10*(j-1):1162+10*(j-1) ) , out % Space_Count_Temp(1,j) ) call b2_to_si2 ( rec ( 1163+10*(j-1):1164+10*(j-1) ) , out % Space_Count_Temp(2,j) ) call b2_to_si2 ( rec ( 1165+10*(j-1):1166+10*(j-1) ) , out % Space_Count_Temp(3,j) ) call b2_to_si2 ( rec ( 1167+10*(j-1):1168+10*(j-1) ) , out % Space_Count_Temp(4,j) ) call b2_to_si2 ( rec ( 1169+10*(j-1):1170+10*(j-1) ) , out % Space_Count_Temp(5,j) ) END do !---------------------------------------------------------------------- ! unpack calibration coefficients !---------------------------------------------------------------------- do ichan = 1, 3 call b4_to_i4_avhrr ( rec (49+(ichan-1)*60:52+(ichan-1)*60) , IS_SIGNED, i4word) out % vis_slope_1 ( ichan) = i4word / (10.0 ** 7) call b4_to_i4_avhrr ( rec (53+(ichan-1)*60:56+(ichan-1)*60) , IS_SIGNED, i4word) out % vis_intercept_1 ( ichan) = i4word / (10.0 ** 6) call b4_to_i4_avhrr ( rec (57+(ichan-1)*60:60+(ichan-1)*60) ,IS_SIGNED, i4word ) out % vis_slope_2 ( ichan) = i4word / (10.0 ** 7) call b4_to_i4_avhrr ( rec (61+(ichan-1)*60:64+(ichan-1)*60) , IS_SIGNED, i4word) out % vis_intercept_2 ( ichan) = i4word / (10.0 ** 6) call b4_to_i4_avhrr ( rec (65+(ichan-1)*60:68+(ichan-1)*60) , IS_SIGNED, i4word) out % vis_intersection ( ichan) = i4word call b4_to_i4_avhrr ( rec ( 229+(ichan-1)*24:232+(ichan-1)*24) , IS_SIGNED, i4word) out % ir_coef_1_1b (ichan + 3) = i4word call b4_to_i4_avhrr ( rec ( 233+(ichan-1)*24:236+(ichan-1)*24) ,IS_SIGNED, i4word ) out % ir_coef_2_1b (ichan + 3) = i4word call b4_to_i4_avhrr ( rec ( 237+(ichan-1)*24:240+(ichan-1)*24) ,IS_SIGNED, i4word ) out % ir_coef_3_1b (ichan + 3) = i4word END do !----------------------------------------------------------------------------- ! handling NOAA-N change to scaling of ir_coef_3_1b for ch4 and ch5 ! the scaling changed from 10**6 to 10**7 - account for extra ! factor of 10 in NOAA-N format here !------------------------------------------------------------------------------ IF ( meta % ver_1b >= 3) then out % ir_coef_3_1b(4) = out % ir_coef_3_1b(4) / 10.0 out % ir_coef_3_1b(5) = out % ir_coef_3_1b(5) / 10.0 ENDIF !---------------------------------------------------------------------- ! unpack number of Anchor points !---------------------------------------------------------------------- ! Num_loc = Buffer_temp(53) Num_loc = 51 !---------------------------------------------------------------------- ! unpack angle data !---------------------------------------------------------------------- do j = 1, 51 ibyte = 329 + (j-1)*6 call b2_to_i2_avhrr ( rec ( ibyte:ibyte+1), IS_SIGNED, i2word ) out % solzen_anchor ( j) = i2word/(10.0 ** 2) call b2_to_i2_avhrr ( rec ( ibyte+2:ibyte+3), IS_SIGNED, i2word ) out % satzen_anchor ( j) =i2word/(10.0 ** 2) call b2_to_i2_avhrr ( rec ( ibyte+4:ibyte+5), IS_SIGNED, i2word) out % relaz_anchor ( j)=i2word/(10.0 ** 2) ENDdo !---------------------------------------------------------------------- ! unpack earth location !---------------------------------------------------------------------- j = 1 do j = 1,51 ibyte = 641 + (j-1)*8 call b4_to_i4_avhrr ( rec ( ibyte : ibyte+3), IS_SIGNED, i4word ) out % lat_anchor(j) = i4word / (10.0 ** 4) call b4_to_i4_avhrr ( rec ( ibyte+4 : ibyte+7), IS_SIGNED, i4word ) out % lon_anchor(j) = i4word / (10.0 ** 4) ENDdo !-- make sure longitudes are correct ( ?? check this) where ( out % lon_Anchor > 180) out % Lon_Anchor = out % Lon_Anchor - 360.0 END where !---------------------------------------------------------------------- ! unpack channel data !---------------------------------------------------------------------- ichan = 1 ipix = 1 do ibyte = 1265 , 1264 + meta % Num_data_bytes , 4 call b4_to_i4_avhrr ( rec ( ibyte : ibyte + 3 ) , IS_UNSIGNED, i4word ) do j = 1,3 IF (ipix > meta % Num_pix) then exit END IF out % Chan_Counts(ichan , ipix ) = ishft(ishft(i4word,2+10*(j-1)),-22) ichan = ichan + 1 !- check which channel is on , record length is only 5 channels, we have here 6 ( 3a and 3b) IF ( ichan == ETchannel_3a ) then IF ( out % Ch3a_On .eqv. .false. ) ichan = ichan + 1 END IF IF ( ichan == ETchannel_3b ) then IF ( out % Ch3a_On ) ichan = ichan + 1 END IF IF (ichan > Num_Chan) then ichan = 1 ipix = ipix + 1 END IF END do IF (ipix > meta % Num_pix) then exit END IF END do !------------------------------------------------------------------------- ! read in cloud mask !----------------------------------------------------------------------- ipix = 0 word_clavr_start = 1264 + meta % Num_data_bytes + meta % Num_digtel_bytes + meta % Num_house_bytes + 1 !------ clavr status flag onebyte = rec(word_clavr_start + 3) out % clavr_status = ishft(ishft(onebyte,7),-7) !--- cloud mask start_word = word_clavr_start + 8 END_word = start_word + meta % Num_clavr_bytes do iword = start_word, END_word onebyte = rec(iword) ipix = 1 start_pixel = (iword - start_word)*4 + 1 END_pixel = min( meta % Num_pix, start_pixel+3_int2) do ipix = start_pixel, END_pixel IF (ipix > meta % Num_pix) then exit END IF i = (ipix - start_pixel) + 1 out % Cld_Mask_Aux(ipix) = ishft(ishft(onebyte,2*(i-1)),-6) END do END do END SUBROUTINE UNPACK_AVHRR_DATA_RECORD_KLM !--------------------------------------------------------------------- ! EXPLANATION OF LAC TO GAC GEOLOCATION CONVERSION: ! ! A LAC scan line is converted to GAC line in the following way: ! ! 1. GAC processing starts from the 3rd LAC pixel; 1st and 2nd LAC pixels are ! not used. ! ! 2. The 1st GAC pixel is formed from the 3rd, 4th, 5th and 6th LAC pixel. ! ! 3. The 7th LAC pixel is skipped and the processing for the 2nd GAC pixel uses ! the 8th, 9th, 10th and 11th LAC pixel. ! ! 4. LAC pixels contributing to a GAC pixel can be described with a formula: ! ! l, l+1, l+2, l+3 for l = 5*g-2 where g = 1,..,409 ! ! `l' and `g' being LAC and GAC pixel index, respectively, in a scan line. ! ! 5. Geolocation of the third LAC pixel of the four used for a GAC pixel is ! assigned to that GAC pixel. For the 1st and 2nd GAC pixels these are the ! 5th and 10th LAX pixels. Or in short: ! ! l = 5*g (g=1,..,409) ! ! 6. The 2048th LAC pixel is not used for GAC pixel processing. The processing ! ENDs with the 2047th LAX pixels, for a total of 2045 LAX pixels used ! (2045/5=409). ! ! This SUBROUTINE performs a five-point Lagrangian interpolation for 2048 LAC ! pixels regardless of the actual length of the scan line (GAC or LAC). IF it's ! a GAC scan line, the geolocations of the 2nd and 3rd LAC pixels for each GAC ! pixel are averaged to yield the final GAC pixel's geolocation. !--------------------------------------------------------------------- SUBROUTINE LAGRANGIAN_ANCHOR_INTERP(Anchor, meta , y) REAL , dimension(:), INTENT(in):: Anchor TYPE ( avhrr_meta_TYPE) , INTENT(in) :: meta REAL , dimension(:), INTENT(out):: y REAL , dimension(2048):: temp_y REAL, dimension(:), allocatable:: temp_Anchor REAL :: x,L,xa_i,xa_j INTEGER:: m,i,j,ia,ia1,ia2,ipix,Num_Anchor INTEGER:: discon m = 2 !2 means 5 point interp Num_Anchor = size(Anchor,1) !should be 51 allocate(temp_Anchor(Num_Anchor)) !--- initialize y = 0.0 temp_y = 0.0 temp_Anchor = Anchor !-------- check for a discontinuity (ie. crossing the 180ºW meridian) discon = 0 IF ((minval(temp_Anchor) < -140.0).and.(maxval(temp_Anchor) > 140.0)) then discon = 1 where (temp_Anchor < 0.0) temp_Anchor = temp_Anchor + 360.0 END where END IF !-- perform 2m+1 point Lagrangian interpolation do ipix = 1, 2048 ia = (ipix-25)/40+1 ia1 = min(Num_Anchor,max(1,ia-m)) ia2 = max(1,min(Num_Anchor,ia+m)) x = ipix do i = ia1,ia2 L = 1.0 xa_i = (i-1)*(40)+25 do j = ia1,ia2 xa_j = (j-1)*(40)+25 IF (i /= j) L = L*(x-xa_j)/(xa_i-xa_j) END do temp_y(ipix) = temp_y(ipix)+L*temp_Anchor(i) END do END do ! IF it's a GAC scan line perform averaging of `temp_y' for the 2nd and 3rd LAC ! pixel going into every GAC pixel for the final value of `y'. IF ( meta % Num_pix == 409) then do ipix = 1, 409 y(ipix) = 0.5*(temp_y(5*ipix)+temp_y(5*ipix-1)) END do ELSE y = temp_y END IF !----- account for any discontinuities where(y > 180.0) y = y - 360.0 END where where(y < -180.0) y = 360.0 + y END where deallocate(temp_Anchor) END SUBROUTINE LAGRANGIAN_ANCHOR_INTERP ! --------------------------------------------------------------------- ! gnomic navigation Anchor interpolation routine ! note, this is only Valid for lat and lon, do not apply to angles !--------------------------------------------------------------------- SUBROUTINE GNOMIC_ANCHOR_INTERP(Lon_Anchor_deg,Lat_Anchor_deg, meta, longitude,latitude) REAL (kind=REAL4), dimension(:), INTENT(in):: Lon_Anchor_deg,Lat_Anchor_deg TYPE ( avhrr_meta_TYPE) , INTENT(in) :: meta REAL (kind=REAL4), dimension(:), INTENT(out):: longitude,latitude REAL (kind=REAL4), dimension(:), allocatable:: k,x_Anchor,y_Anchor, & Lon_Anchor,Lat_Anchor REAL (kind=REAL4), dimension(:), allocatable:: x,y REAL (kind=REAL4):: yy,xx,xa1,xa2,rho,v,v_prev,c,Lat_center,Lon_center,offset INTEGER:: ia1,ia2,ipix,Num_Anchor,Num_skip,Num_start REAL ( kind = REAL4) , parameter :: PI = 3.14159 REAL ( kind = REAL4) , parameter :: DTOR = PI/180. Num_Anchor = size(Lon_Anchor_deg) !should be 51 Num_skip = meta % Num_pix / Num_Anchor Num_start = ( meta % Num_pix - Num_skip*(Num_Anchor-1))/2 + 1 allocate(k(Num_Anchor)) allocate(x_Anchor(Num_Anchor)) allocate(y_Anchor(Num_Anchor)) allocate(Lon_Anchor(Num_Anchor)) allocate(Lat_Anchor(Num_Anchor)) allocate ( x( meta % num_pix) , y(meta% num_pix)) !--- initialize latitude = 0.0 longitude = 0.0 !--- convert Anchors to radians Lat_Anchor = Lat_Anchor_deg * dtor Lon_Anchor = Lon_Anchor_deg * dtor !--- pick center for projection (in radians) Lat_center = Lat_Anchor(Num_Anchor/2) Lon_center = Lon_Anchor(Num_Anchor/2) !--- convert to gnomic space k = sin(Lat_center) * sin(Lat_Anchor) + & cos(Lat_center) * cos(Lat_Anchor) * cos(Lon_Anchor - Lon_center) k = 1.0 / k x_Anchor = k * cos(Lat_Anchor) * sin(Lon_Anchor - Lon_center) y_Anchor = k * (cos(Lat_center) * sin(Lat_Anchor) - & sin(Lat_center) * cos(Lat_Anchor) * cos(Lon_Anchor - Lon_center)) !--- linear interp in gnomic space do ipix = 1, meta % Num_pix !--- determine two closest Anchor points ia1 = max(1,min(Num_Anchor-1,(ipix - Num_start)/Num_skip + 1)) ia2 = ia1 + 1 !--- determine weights of two closest Anchor points xa1 = (ia1-1)*(Num_skip) + Num_start xa2 = (ia2-1)*(Num_skip) + Num_start !--- determine point to interpolate to xx = ipix !--- interp x at xx yy = (xx-xa1)*(x_Anchor(ia2)-x_Anchor(ia1))/(xa2-xa1) + & x_Anchor(ia1) x(ipix) = yy !--- interp y at xx yy = (xx-xa1)*(y_Anchor(ia2)-y_Anchor(ia1))/(xa2-xa1) + & y_Anchor(ia1) y(ipix) = yy END do !--- convert other quantities v_prev = 0.0 offset = 0 do ipix = 1, meta % Num_pix !--- some needed constants rho = sqrt(x(ipix)**2 + y(ipix)**2) !--- handle center point where rho = 0.0 IF (rho == 0.0) then latitude(ipix) = Lat_center longitude(ipix) = Lon_center cycle END IF c = atan(rho) !--- v is the longitude in radians relative to the center longitude v = atan( x(ipix) * sin(c) / & (rho * cos(Lat_center) * cos(c) - y(ipix) * sin(Lat_center) * sin(c)) ) !--- look for a jump and determine offset IF (ipix == 1) then v_prev = v END IF latitude(ipix) = asin(cos(c) * sin(Lat_center) + (y(ipix)*sin(c)*cos(Lat_center) / rho)) longitude(ipix) = Lon_center + v v_prev = v END do !--- convert to degrees latitude = latitude / dtor longitude = longitude / dtor !----- account for discontinuity IF there was one where (longitude > 180.0) longitude = longitude - 360.0 END where where (longitude < -180.0) longitude = longitude + 360.0 END where deallocate(k) deallocate(x_Anchor) deallocate(y_Anchor) deallocate(Lon_Anchor) deallocate(Lat_Anchor) deallocate(x) deallocate(y) END SUBROUTINE GNOMIC_ANCHOR_INTERP SUBROUTINE alloc_out ( self , nx , ny) class ( avhrr_data_out ) :: self INTEGER, INTENT(in) :: nx , ny INTEGER:: i IF ( self % is_allocated .eqv. .false.) then allocate ( self % geo % solzen ( nx, ny)) allocate ( self % geo % satzen ( nx, ny)) allocate ( self % geo % solaz ( nx, ny)) allocate ( self % geo % sataz ( nx, ny)) allocate ( self % geo % relaz ( nx, ny)) allocate ( self % geo % lat ( nx, ny)) allocate ( self % geo % lon ( nx, ny)) allocate ( self % geo % scan_time ( ny) ) allocate ( self % geo % ascEND ( ny) ) do i = 1 , 6 allocate ( self % band(i) % ref ( nx, ny )) allocate ( self % band(i) % rad ( nx, ny )) allocate ( self % band(i) % bt ( nx, ny )) END do self % is_allocated = .true. ELSE print*,'was already allocated' print*, 'check program' print*,'mail andi.walther@ssec.wisc.edu' END IF END SUBROUTINE alloc_out SUBROUTINE dealloc_out ( self ) class ( avhrr_data_out ) :: self INTEGER:: i IF ( self % is_allocated ) then deallocate ( self % geo % solzen ) deallocate ( self % geo % satzen ) deallocate ( self % geo % solaz ) deallocate ( self % geo % sataz ) deallocate ( self % geo % relaz ) deallocate ( self % geo % lat ) deallocate ( self % geo % lon ) deallocate ( self % geo % scan_time ) deallocate ( self % geo % ascEND ) do i = 1 , 5 deallocate ( self % band(i) % ref ) deallocate ( self % band(i) % rad ) deallocate ( self % band(i) % bt ) END do self % is_allocated = .false. ELSE print*,'was not allocated' print*, 'check program' print*,'mail andi.walther@ssec.wisc.edu' END IF END SUBROUTINE dealloc_out !======================================================================================= ! ! ! ! !======================================================================================== SUBROUTINE reflectance_calibration ( chan_counts & , ref_cal_1b & , therm_cal_1b & , data_line & , time_since_launch & , out) INTEGER( kind = int2) :: chan_counts (:,:,:) LOGICAL , INTENT(in) :: Ref_Cal_1b LOGICAL , INTENT(in) :: therm_cal_1b TYPE ( avhrr_line_TYPE ), INTENT(in) :: data_line REAL(kind=REAL4), INTENT(in):: time_since_launch TYPE ( band_str) , dimension (3) :: out REAL(kind=REAL4) :: gain_low (3) REAL(kind=REAL4) :: dark_count_cal (3) REAL(kind=REAL4) :: gain_high (3) REAL(kind=REAL4) :: switch_count_cal (3) REAL(kind=REAL4) :: ref_switch ( 3) INTEGER:: i_chn , iline LOGICAL :: Valid_Ref_Cal print*,ref_cal_1b,therm_cal_1b ! - reflectance from avhrr-file or from instrument file ?! IF ( ref_cal_1b ) then ! IF so read the slope from input file do i_chn = 1 , 3 gain_low (i_chn) = data_line % vis_slope_1(i_chn) dark_count_Cal(i_chn) = -1.0*data_line % vis_intercept_1(i_chn) / data_line % vis_slope_1(i_chn) gain_high (i_chn) = data_line%vis_slope_2(i_chn) switch_count_Cal(i_chn) = data_line%vis_intersection(i_chn) ref_switch (i_chn) = data_line%vis_slope_2(i_chn) * data_line%vis_intersection(i_chn) & & + data_line%vis_intercept_2(i_chn) END do print*,gain_low (1), dark_Count_Cal(1), gain_high(1) , switch_Count_Cal (1), ref_switch (1) ELSE ! now we need the instr file !--- adjust gains to account for degradation IF (Therm_Cal_1b ) then do i_chn = 1 , 3 Dark_Count_Cal(i_chn) = instr_cons% Dark_Count (i_chn) END do ELSE ! Ch1_Dark_Count_Cal = Scan_Space_Counts_Avhrr(1,iline) ! Ch2_Dark_Count_Cal = Scan_Space_Counts_Avhrr(2,iline) ! Ch3a_Dark_Count_Cal = Scan_Space_Counts_Avhrr(3,iline) ENDIF !-- set switch counts to values in instrument files do i_chn = 1 , 3 Switch_Count_Cal(i_chn) = instr_cons% Switch_Count (i_chn) Gain_low(i_chn) = instr_cons % Gain_Low_0(i_chn) & & * (100.0 + instr_cons % Degrad_Low_1(i_chn) * time_since_launch + & instr_cons % Degrad_Low_2(i_chn)*time_since_launch**2) & & /100.0 Gain_High(i_chn) = instr_cons%Gain_High_0(i_chn)*(100.0+instr_cons%Degrad_High_1(i_chn)*time_since_launch+ & instr_cons%Degrad_High_2(i_chn)*time_since_launch**2)/100.0 ref_switch (i_chn) = data_line%vis_slope_2(i_chn) * data_line%vis_intersection(i_chn) + data_line%vis_intercept_2(i_chn) END do print*,gain_low (1), dark_Count_Cal(1), gain_high(1) , switch_Count_Cal (1), ref_switch (1) !--- test for Valid dark count IF (abs(Dark_Count_Cal(1) - instr_cons % Dark_Count ( 1 ) ) > 5) then Valid_Ref_Cal = .false. END IF END IF do i_chn = 1 , 3 where (Chan_Counts (i_chn,:,:) <= Switch_Count_Cal(i_chn) ) out ( i_chn) % Ref = Gain_low ( i_chn) *(Chan_Counts (i_chn,:,:) - Dark_Count_Cal( i_chn)) ELSE where out ( i_chn) % Ref = Ref_Switch( i_chn) + Gain_High( i_chn)*(Chan_Counts (i_chn,:,:) - Switch_Count_Cal (i_chn) ) END where END do !- TODO ! do iline = Line_Idx_Min_Segment,Line_Idx_Min_Segment + Num_Scans_Read - 1 !--- check for a bad scan, IF so then skip ! IF (Bad_Scan_Flag(iline) == 1) then ! Ref_Ch1(:,iline) = Missing_Value_REAL4 ! Ref_Ch2(:,iline) = Missing_Value_REAL4 ! Ref_Ch6(:,iline) = Missing_Value_REAL4 ! END IF !--- set Ref_Ch6 to missing for Ch3a = off ! IF (Ch3a_On(iline) == sym%NO) then ! Ref_Ch6(:,iline) = Missing_Value_REAL4 ! END IF END SUBROUTINE reflectance_calibration !--------------------------------------------------------------- ! read the avhrr constants into memory !----------------------------------------------------------------- SUBROUTINE READ_AVHRR_INSTR_CONSTANTS( self, Instr_Const_file ) class ( instr_const_TYPE ) :: self CHARACTER(len=*), INTENT(in):: Instr_Const_file INTEGER:: ios0, erstat INTEGER:: Instr_Const_lun !Instr_Const_lun = GET_LUN() Instr_Const_lun = 10 open(unit=Instr_Const_lun,file=trim(Instr_Const_file),status="old",position="rewind",action="read",iostat=ios0) erstat = 0 IF (ios0 /= 0) then erstat = 19 print *, "Error opening AVHRR constants file, ios0 = ", ios0 print*, 'file : ',trim(Instr_Const_file) stop 19 END IF read(unit=Instr_Const_lun,fmt="(a3)") self % sat_name read(unit=Instr_Const_lun,fmt=*) self % Solar_Ch20 read(unit=Instr_Const_lun,fmt=*) self % Ew_Ch20 read(unit=Instr_Const_lun,fmt=*) self % Space_Rad_3b, self % b0 (4),self % b1(4), self % b2 (4) read(unit=Instr_Const_lun,fmt=*) self % Space_Rad_4,self % b0 (5),self % b1(5),self % b2 (5) read(unit=Instr_Const_lun,fmt=*) self % Space_Rad_5, self % b0 (6) ,self % b1 (6), self % b2 (6) read(unit=Instr_Const_lun,fmt=*) self % nu_20 read(unit=Instr_Const_lun,fmt=*) self % a1_20 read(unit=Instr_Const_lun,fmt=*) self % a2_20 read(unit=Instr_Const_lun,fmt=*) self % nu_31 read(unit=Instr_Const_lun,fmt=*) self % a1_31 read(unit=Instr_Const_lun,fmt=*) self % a2_31 read(unit=Instr_Const_lun,fmt=*) self % nu_32 read(unit=Instr_Const_lun,fmt=*) self % a1_32 read(unit=Instr_Const_lun,fmt=*) self % a2_32 read(unit=Instr_Const_lun,fmt=*) self % Dark_Count(1) read(unit=Instr_Const_lun,fmt=*) self % Dark_Count(2) read(unit=Instr_Const_lun,fmt=*) self % Dark_Count(3) read(unit=Instr_Const_lun,fmt=*) self % Gain_Low_0(1),self % Degrad_Low_1(1), self % Degrad_Low_2(1) read(unit=Instr_Const_lun,fmt=*) self % Gain_High_0(1),self % Degrad_High_1(1),self % Degrad_High_2(1) read(unit=Instr_Const_lun,fmt=*) self % Gain_Low_0(2),self % Degrad_Low_1(2), self % Degrad_Low_2(2) read(unit=Instr_Const_lun,fmt=*) self % Gain_High_0(2),self % Degrad_High_1(2),self % Degrad_High_2(2) read(unit=Instr_Const_lun,fmt=*) self % Gain_Low_0(3),self % Degrad_Low_1(3),self % Degrad_Low_2(3) read(unit=Instr_Const_lun,fmt=*) self % Gain_High_0(3),self % Degrad_High_1(3),self % Degrad_High_2(3) read(unit=Instr_Const_lun,fmt=*) self % launch_date read(unit=Instr_Const_lun,fmt=*) self % Switch_Count(1) read(unit=Instr_Const_lun,fmt=*) self % Switch_Count(2) read(unit=Instr_Const_lun,fmt=*) self % Switch_Count(3) read(unit=Instr_Const_lun,fmt=*) self % Prt_Coef(:,1) read(unit=Instr_Const_lun,fmt=*) self % Prt_Coef(:,2) read(unit=Instr_Const_lun,fmt=*) self % Prt_Coef(:,3) read(unit=Instr_Const_lun,fmt=*) self % Prt_Coef(:,4) read(unit=Instr_Const_lun,fmt=*) self % prt_weight(:) read(unit=Instr_Const_lun,fmt=*) self % b1_day_mask,self % b2_day_mask,self % b3_day_mask,self % b4_day_mask close(unit=Instr_Const_lun) !-- convert solar flux in channel 3b to mean with units mW/m^2/cm^-1 self % Solar_Ch20_Nu = 1000.0 *self % Solar_Ch20 /self % Ew_Ch20 !----- define planck constants here END SUBROUTINE READ_AVHRR_INSTR_CONSTANTS !----------------------------------------------------------- ! THERMAL CALIBRATION ! ! variables ! therm_Cal_1b - use 1b cal. coefficients (true) ! klm - LOGICAL flag set to true IF this klm data ! Ch3a_On_temp - INTEGERflag set to 0 IF Ch3a is off ! count3 - vector of channel 3 counts ! count4 - vector of channel 4 counts ! count5 - vector of channel 5 counts ! d_coef - the "d" coefficient as illustrated in klm-POD ! a_coef - the "a" coefficient as illustrated in klm-POD ! b_coef - the "b" coefficient as illustrated in klm-POD ! linear_slope - slope in linear radiance computation ! linear_intercept - intercept in linear radiance computation ! nu_20, a1_20, a2_20 - radiance to temperature conversion parameters ! (see klm-PODUG for example) ! nu_31, a1_31, a2_31 - radiance to temperature conversion parameters ! nu_32, a1_32, a2_32 - radiance to temperature conversion parameters ! ! ! output ! Rad_Ch20_temp, Rad_Ch31_temp, Rad_Ch32_temp - vector of channel radiance in mW/m^2/str/cm^-1 ! Bt_Ch20_temp, Bt_Ch31_temp, Bt_Ch32_temp - vector of brightness temperature in K ! ! Some Notes ! 1- pre-klm data is first converted to a radiance assuming a linear calibration ! then the non-linear adjustment is made ! 2- Does not incorporate the correct procedure for data ! before NOAA-14. Need to review original CLAVR-1 for this. ! 3 - klm data, the non-linear correction is added to the linear radiance ! but for non-klm data is the non-linear correction is additive (see pods) ! this is reflected in the dIFference is the pre-klm and klm coefs. ! ! Revision History ! Coded January 2003 - A. Heidinger ! February 2003 - added therm_Cal_1b argument and new calibration ! ! ! some clarIFication ! the klm-POD nomenclature has changed. ! !--------------------------------------------------------------------------- SUBROUTINE thermal_calibration ( therm_cal_1b, data_line ,meta, chan_on, time_since_launch, out) LOGICAL , INTENT(in) :: therm_cal_1b TYPE ( avhrr_line_TYPE ), INTENT(in) :: data_line TYPE ( avhrr_meta_TYPE ), INTENT(in) :: meta LOGICAL ,INTENT(in) :: chan_on ( 6 ) REAL(kind=REAL4), INTENT(in):: time_since_launch REAL (kind= REAL4), dimension (4:6,meta % num_pix ) :: out INTEGER(kind=int4):: i_chan INTEGER(kind=int4):: j INTEGER, parameter:: dx = 2 INTEGER, parameter:: dy = 0 INTEGER:: nvalid INTEGER:: Elem_Idx_Min INTEGER:: Elem_Idx_Max INTEGER:: Line_Idx_Max INTEGER:: Line_Idx_Min INTEGER:: Elem_Idx INTEGER, parameter:: Ch20_Counts_Filter_Thresh = 900 REAL :: c0(4:6) , c1(4:6) , c2 (4:6) LOGICAL :: is_KLM !------------------------------------------------------------------------ ! Perform Filter on Ch20 counts for appropriate sensors ! Do this only for TIROS-N, NOAA-5,6,7,8,9,10 !------------------------------------------------------------------------ ! IF ((Sc_Id_WMO >= 706 .and. Sc_Id_WMO <= 708) .or. (Sc_Id_WMO >= 200 .and. Sc_Id_WMO <= 202)) then ! Temp_Pix_Array = Chan_Counts_Avhrr(3,:,:) ! Ch20_Counts_Filtered = Chan_Counts_Avhrr(3,:,:) ! One_Byte_Temp = 1 ! where(Temp_Pix_Array <= 100 .or. Temp_Pix_Array > 1024) ! One_Byte_Temp = 0 ! END where ! do Elem_Idx = 1, Num_Pix ! Elem_Idx_Min = min(Num_Pix,max(1,Elem_Idx - dx)) ! Elem_Idx_Max = min(Num_Pix,max(1,Elem_Idx + dx)) ! do Line_Idx = 1, Num_Scans_Read ! IF (Chan_Counts_Avhrr(3,Elem_Idx,Line_Idx) > Ch20_Counts_Filter_Thresh) then ! Line_Idx_Min = min(Num_Scans_Read,max(1,Line_Idx - dy)) ! Line_Idx_Max = min(Num_Scans_Read,max(1,Line_Idx + dy)) ! nvalid = sum(One_Byte_Temp(Elem_Idx_Min:Elem_Idx_Max,Line_Idx_Min:Line_Idx_Max)) ! IF (nvalid > 1) then ! Ch20_Counts_Filtered(Elem_Idx,Line_Idx) = & ! sum(One_Byte_Temp(Elem_Idx_Min:Elem_Idx_Max,Line_Idx_Min:Line_Idx_Max)* & ! Temp_Pix_Array(Elem_Idx_Min:Elem_Idx_Max,Line_Idx_Min:Line_Idx_Max)) / nvalid ! END IF ! END IF !END do !END do !Chan_Counts_Avhrr(3,:,:) = int(Ch20_Counts_Filtered,kind=int2) !END IF !------------------------------------------------------------------------ ! Loop through each line and calibrate !------------------------------------------------------------------------ c0 (:) = data_line % ir_coef_1_1b( ETchannel_FIRST_NIR : ETchannel_LAST_NIR ) c1 (:) = data_line % ir_coef_2_1b( 4:6 ) c2 (:) = data_line % ir_coef_3_1b( 4:6 ) ! - make it very simple! is_KLM = meta % avhrr_TYPE == ET_avhrr_t % KLM IF ( is_KLM .and. therm_Cal_1b) then !apply level-1b non-linear correction c0 (:) = data_line % ir_coef_1_1b( ETchannel_FIRST_NIR : ETchannel_LAST_NIR ) c1 (:) = data_line % ir_coef_2_1b( 4:6 ) c2 (:) = data_line % ir_coef_3_1b( 4:6 ) ELSE IF ( is_KLM .eqv. .false. .and. therm_Cal_1b) then !compute level-1b linear correction for pre-klm ! c0 (:) = ir_linear_intercept_1b( 4:6 ) ! c1 (:) = ir_linear_slope_1b( 4:6 ) ! c2 (:) = 0. ELSE IF ( therm_Cal_1b .eqv. .false.) then !generate linear radiance using internally computed coefficients ! c0 (:) = ir_linear_intercept_new( 4:6 ) ! c1 (:) = Ir_Linear_Slope_New( 4:6 ) ! c2 (:) = 0. END IF do i_chan = 4 , 6 IF ( chan_on ( i_chan ) ) then out ( i_chan, : ) = c0 (i_chan) & + c1 (i_chan) * data_line % Chan_Counts ( i_chan , : ) & + c2(i_chan) * data_line % Chan_Counts ( i_chan , : )**2 END IF END do !--- apply non-linear correction where appropriate (non-klm or any IF thermal cal redone) IF ( ( is_KLM .eqv. .false.) .or.( therm_Cal_1b .eqv. .false. )) then do i_chan = 4 , 6 out ( i_chan, : ) = out ( i_chan , : ) & + instr_cons%b0 ( i_chan) & + instr_cons%b1 (i_chan ) * data_line % Chan_Counts ( i_chan , : ) & + instr_cons%b2 ( i_chan) * data_line % Chan_Counts ( i_chan , : ) ** 2 END do END IF END SUBROUTINE thermal_calibration END module avhrr_read_mdl