C Copyright(c) 2005, Space Science and Engineering Center, UW-Madison C Refer to "McIDAS Software Acquisition and Distribution Policies" C in the file mcidas/data/license.txt C *** $Id: kbxinst.dlm,v 1.2 2015/03/31 19:15:21 rickk Tst $ *** C$ Name: C$ kbxini - Calibration module INSAT-3D C$ C$ Interface: C$ integer function C$ kbxini(character*4 input_cal,character*4 output_cal,integer io_size(*)) C$ C$ Input: C$ input_cal - Input calibration type ('RAW') C$ output_cal - Output calibration type ('TEMP','ALB','RAD','RAW','BRIT') C$ io_size - Source and destination byte depth C$ io_size(1) - Source pixel size (2) C$ io_size(2) - Destination pixel size (1,2,4) C$ C$ Output: C$ Error messages if return value is not 0. (To be backward compatible C$ with the existing API, the return value does not indicate which C$ error was encountered.) C$ C$ Return values: C$ 0 - success C$ -1 - error C$ C$ Remarks: C$ Note also that for ADDE we need to comment out all SDEST and DDEST C$ calls, since messages cannot be sent to "stdout" pipe mixed with data. C$ C$ EDEST messages will be sent to "stderr" and can be left in, but will C$ not be returned to a remote client. C$ C$ C$ Categories: C$ image C$ display C$ calibration C$ met/science integer function kbxini(input_cal,output_cal,io_size) implicit none c c --- input parametes c character*4 input_cal ! input calibration type (RAW) character*4 output_cal ! output cal type (TEMP,ALB,RAD,RAW,BRIT) integer io_size(*) ! source and destination byte depth c c --- local variables c c --- INSAT-3D calibration block character*4 cal_input ! input calibration type (RAW) character*4 cal_output ! output cal type (TEMP, ALB, RAD, RAW, BRIT) integer src_byte_size ! number of bytes for source pixel integer dest_byte_size ! number of bytes for destination pixel integer cal_flag ! 0= from AREA, 1= from server integer cal_block(684) ! cal block - includes directory and tables common/inss/ & cal_input, & cal_output, & src_byte_size, & dest_byte_size, & cal_flag, & cal_block c c --- set local variables to fill common block c cal_input = input_cal cal_output = output_cal c c --- Verify that the input calibration type c if (cal_input(1:3) .ne. 'RAW') then kbxini = -1 c c --- Verify that the output calibration type c elseif ((cal_output(1:4) .ne. 'TEMP') .and. & (cal_output(1:3) .ne. 'RAD') .and. & (cal_output(1:3) .ne. 'RAW') .and. & (cal_output(1:4) .ne. 'ALB') .and. & (cal_output(1:4) .ne. 'BRIT')) then kbxini = -1 else src_byte_size = io_size(1) dest_byte_size = io_size(2) cal_flag = 0 kbxini = 0 endif return end C$ Name: C$ kbxcal - Converts input DN to output of correct byte size & dimensions C$ C$ Interface: C$ integer function C$ kbxcal(integer prefix(*), integer area_dir(*), integer num_pixels, C$ integer band, integer ibuf(*) ) C$ C$ Input: C$ prefix - Line prefix calibration information C$ area_dir - Area directory associated with the IBUF data C$ num_pixels - Number of pixels in IBUF array C$ band - Band number for data in IBUF array C$ ibuf - array of data to be converted C$ C$ Output: C$ ibuf - array of converted data C$ C$ Return values: C$ 0 - success C$ -1 - failure to build a calibration lookup table C$ (output will always be returned in IBUF, even with an C$ invalid or unbuilt table) C$ C$ Remarks: C$ KBXINI determines if the conversion is permitted, by putting C$ constraints on the conversion types. Here, the type is assumed C$ already validated, and the only decision is whether to build a C$ table for a new valid type or use the table already existing. C$ Since data is always returned, it is the user's responsibility to C$ test return value KBXCAL for 0 to guarantee returned data is valid. C$ C$ Categories: C$ image C$ display C$ calibration C$ met/science integer function kbxcal(prefix,area_dir,num_pixels,band,ibuf) implicit none include 'areaparm.inc' ! defines NUMAREAOPTIONS integer MAX_SNDR_IR_VAL ! Maximum raw value for SNDR IR integer MAX_SNDR_VIS_VAL ! Maximum raw value for SNDR VIS integer MAX_IMGR_VAL ! Maximum raw value for IMGR IR integer MAX_VAL c c --- Currently Sounder has a maximum of 16383 values c --- and imager 1024 values c parameter (MAX_SNDR_IR_VAL=16383) parameter (MAX_SNDR_VIS_VAL=4096) parameter (MAX_IMGR_VAL=1024) parameter (MAX_VAL=16383) c c --- input parameters c integer prefix(*) ! line prefix information to determine detector integer area_dir(*) ! area directory integer num_pixels ! number of pixels in ibuf array integer band ! band number integer ibuf(*) ! I/O array containing pixels to be modified c c --- external functions c real rad_to_temp ! converts radiance to temperature character*12 cfi character*12 cfe character*12 cfd character*4 clit c c --- local variables c character*4 vis_type character*72 debug character*12 cval character*12 cval2 character*12 cval3 character*12 cval4 character*12 cval5 integer area_number ! area number found in area directory integer cal_offset ! offset into calibration block integer i ! do loop index integer j ! do loop index integer invert_flag integer offset ! band based offset to get to coeficients integer last_area ! last area used for araget integer last_band ! last area used for araget integer NUM_VAL integer alb_brit(MAX_VAL) ! brightness table determined from albedo to send to MPIXTB integer ssnum ! satellite ID (word 3 in area directory) integer temp_brit(MAX_VAL) ! brightness table determined from temperatures integer temp_table(MAX_VAL) ! temperature table determined from raw values integer x ! do loop index integer alb_table(MAX_VAL) ! albedo table determined from raw values integer vis_rad_table(MAX_VAL) ! albedo table determined from raw values integer rad_table(MAX_VAL) ! radiance table determined from temperatures real albedo ! albedo calculated emperically real radiance ! radiance for rad_to_temp real vis_radiance ! visible radiance double precision squared_coef ! squared coefficient for quadratic equation to determin radiance of albedo double precision linear_coef ! linear coefficient for quadratic equation to determin radiance of albedo double precision constant_coef ! constant coefficient for quadratic equation to determin radiance of albedo integer squared_factor ! factor of 10 to determine squared coefficient integer linear_factor ! factor of 10 to determine linear coefficient integer constant_factor ! factor of 10 to determine constant coefficient integer temp_count real temperature ! temperature read in from calibration table real wavelength ! central wavelength for band integer wavelength_factor ! central wavelength for band logical ir_data ! flag indicating ir data requested logical vis_data ! flag indicating visible data requested c --- INSAT-3D calibration block character*4 cal_input ! input calibration type RAW character*4 cal_output ! output cal type (TEMP, ALB, RAD, RAW, BRIT) integer src_byte_size ! number of bytes for source pixel integer dest_byte_size ! number of bytes for destination pixel integer cal_flag ! 0= from AREA, 1=from server integer cal_block(684) ! cal block - includes directory common/inss/ & cal_input, & cal_output, & src_byte_size, & dest_byte_size, & cal_flag, & cal_block data last_area /-1/ data last_band /-1/ kbxcal = 0 ssnum = area_dir(3) c c --- set ir_data/vis_data flag, odd numbered SS values are Imager c --- even numbered SS values are Sounder c if (mod(ssnum,2) .eq. 0) then if (band .lt. 19) then vis_data = .FALSE. ir_data = .TRUE. NUM_VAL = MAX_SNDR_IR_VAL else vis_data = .TRUE. ir_data = .FALSE. NUM_VAL = MAX_SNDR_VIS_VAL endif else NUM_VAL = MAX_IMGR_VAL if (band .gt. 1) then vis_data = .FALSE. ir_data = .TRUE. else vis_data = .TRUE. ir_data = .FALSE. endif endif c c --- read in calibration block c area_number = area_dir(33) cal_offset = area_dir(63) if (last_area .ne. area_number) then if( cal_flag .eq. 0 ) then call araget(area_number,cal_offset,684,cal_block) endif c c --- INSAT-3D uses a quadratic equation to convert raw c --- values to physical quantities. Set each of the values c offset=(band-1)*9 do 999 i=offset,offset+9 cval = cfi(cal_block(i)) call mctrace(1,'KBXINST','cal_block('//cfi(i)//')='//cval) 999 continue squared_factor = cal_block(offset+2) squared_coef = float(cal_block(offset+1))/ & 10.0**squared_factor linear_factor = cal_block(offset+4) linear_coef = float(cal_block(offset+3))/ & 10.0**linear_factor constant_factor = cal_block(offset+6) constant_coef = float(cal_block(offset+5))/ & 10.0**constant_factor invert_flag = cal_block(offset+7) wavelength_factor = cal_block(offset+9) wavelength = float(cal_block(offset+8)) / & 10.0**wavelength_factor if (ir_data) then call mctrace(1,'KBXINST',' ir data') else call mctrace(1,'KBXINST',' vis data') endif cval = cfi(band) call mctrace(1,'KBXINST','band number ='//cval) cval = cfi(cal_offset) call mctrace(1,'KBXINST','calibration offset ='//cval) cval = cfi(cal_block(offset+1)) call mctrace(1,'KBXINST','squared_coefficient ='//cval) cval = cfi(squared_factor) call mctrace(1,'KBXINST','squared_factor ='//cval) cval = cfi(cal_block(offset+3)) call mctrace(1,'KBXINST','linear_coefficient ='//cval) cval = cfi(linear_factor) call mctrace(1,'KBXINST','linear_factor ='//cval) cval = cfi(cal_block(offset+5)) call mctrace(1,'KBXINST','constant_coefficient ='//cval) cval = cfi(constant_factor) call mctrace(1,'KBXINST','constant_factor ='//cval) cval = cfi(invert_flag) call mctrace(1,'KBXINST','invert_flag ='//cval) cval = cfe(wavelength, 6) call mctrace(1,'KBXINST','wavelength ='//cval) cval = cfd(squared_coef, 20) call mctrace(1,'KBXINST','calc squared term ='//cval) cval = cfd(linear_coef, 20) call mctrace(1,'KBXINST','calc linear term ='//cval) cval = cfd(constant_coef, 20) call mctrace(1,'KBXINST','calc constant term ='//cval) endif last_area = area_number if ((vis_data) .and. (last_band .ne. band)) then call mctrace(1,'KBXINST','making vis table') do 30 i=1,NUM_VAL if (invert_flag .eq. 1) then temp_count = (NUM_VAL-i) else temp_count = i endif vis_radiance = squared_coef*(float(temp_count-1))**2+ & linear_coef*float(temp_count-1) + & constant_coef vis_rad_table(i) = nint(vis_radiance * 10000.0) albedo=(i*100.0)/1024.0 alb_table(i) = nint(albedo*10.0) alb_brit(i) = nint(25.5*sqrt(albedo)) c c --- Make sure that brightness values stay between 0 and 255 c if (alb_brit(i) .gt. 255) then alb_brit(i) = 255 elseif (alb_brit(i) .lt. 0) then alb_brit(i) = 0 endif c if (mod(i,100) .eq. 0) then c cval = cfi(alb_brit(i)) c cval2 = cfe(vis_radiance, 6) c cval3 = cfd(squared_coef,6) c cval4 = cfd(linear_coef,6) c cval5 = cfd(constant_coef,6) c call mctrace(1,'KBXINST',' '// c & ' '//cfi(i)// c & ' '//cval // c & ' '//cval2 // c & ' '//cval3 // c & ' '//cval4 // c & ' '//cval5 // c & ' ') c endif 30 continue endif c c --- if data is from ir band c --- create temp, radiance and brit lookup tables c if ((ir_data) .and. (last_band .ne. band)) then c c --- create temperature, radiance and brightness tables c do 110 i=1,NUM_VAL if (invert_flag .eq. 1) then temp_count = (NUM_VAL-i) else temp_count = i endif radiance = squared_coef*(float(temp_count-1))**2 + & linear_coef*float(temp_count-1) + & constant_coef temperature = rad_to_temp(radiance,wavelength) c if (mod(i,100) .eq. 0) then cval = cfe(radiance, 6) cval2 = cfe(temperature, 6) call mctrace(1,'KBXINST',cfi(i) // & ' '//cval//' '//cval2) c endif temp_table(i) = nint(temperature * 100.) rad_table(i) = nint(radiance * 1000.) call gryscl(temperature,temp_brit(i)) 110 continue endif c c --- Convert ibuf to appropriate quantity c if (vis_data) then if (cal_output(1:3) .eq. 'RAW') then call mpixel(num_pixels,src_byte_size, & dest_byte_size,ibuf) elseif (cal_output(1:4) .eq. 'RAD') then call mpixtb(num_pixels,src_byte_size, & dest_byte_size,ibuf,vis_rad_table) elseif (cal_output(1:4) .eq. 'ALB') then call mpixtb(num_pixels,src_byte_size, & dest_byte_size,ibuf,alb_table) elseif (cal_output(1:4) .eq. 'BRIT') then call mpixtb(num_pixels,src_byte_size, & dest_byte_size,ibuf,alb_brit) endif endif if (ir_data) then if (cal_output(1:3) .eq. 'RAW') then call mpixel(num_pixels,src_byte_size, & dest_byte_size,ibuf) elseif (cal_output(1:3) .eq. 'RAD') then call mpixtb(num_pixels,src_byte_size, & dest_byte_size,ibuf,rad_table) elseif (cal_output(1:4) .eq. 'TEMP') then call mpixtb(num_pixels,src_byte_size, & dest_byte_size,ibuf,temp_table) elseif (cal_output(1:4) .eq. 'BRIT') then call mpixtb(num_pixels,src_byte_size, & dest_byte_size,ibuf,temp_brit) endif endif last_band = band return end C$ Name: C$ kbxopt - Returns auxiliary parameters for setup or sets internal state C$ C$ Interface: C$ integer function C$ kbxopt(character*4 option, integer param_in(*), integer param_out(*)) C$ C$ Input: C$ option - Option ('KEYS', 'INFO', 'CALB' ) C$ param_in - Array of input parameters C$ C$ Output: C$ param_out - Array of output parameters C$ (The number of input and output parameters is determined C$ by the option or function executed. The calling routine C$ is responsible for having arrays large enough to send C$ and receive output for the option or function requested.) C$ C$ Return values: C$ 0 - success C$ -1 - Invalid function requested C$ -3 - Error in breakpoint table (table cannot be set up) C$ C$ Remarks: C$ Check the d.pgm code to see how the various functions are used. C$ C$ Categories: C$ image C$ display C$ calibration C$ met/science integer function kbxopt(option,param_in,param_out) implicit none c c --- input parameters c character*4 option ! option integer param_in(*) ! input parameters c c --- output parameters c integer param_out(*) ! output parameters c c --- external functions c character*12 cfi integer brkset ! checks SU table integer ischar ! checks for character string integer lit ! four byte integer representation of char*4 c c --- local variables c character*8 su_file ! stretch table file integer band ! band number integer ssnum ! satellite sensor number logical ir_data ! flag indicating ir data requested logical vis_data ! flag indicating visible data requested c c --- INSAT-3D calibration block c character*4 cal_input ! input calibration type (RAW) character*4 cal_output ! output cal type (TEMP, ALB, RAD, RAW, BRIT) integer src_byte_size ! number of bytes for source pixel integer dest_byte_size ! number of bytes for destination pixel integer cal_flag ! 0=from AREA, 1=from server integer cal_block(684) ! cal block - includes directory and tables integer i common/inss/ & cal_input, & cal_output, & src_byte_size, & dest_byte_size, & cal_flag, & cal_block kbxopt = 0 c c --- set ir_data/vis_data flag c if (option .eq. 'KEYS') then c c --- param_in contains frame directory - set ir_data/vis_data flag c band = param_in(4) call mctrace(1,'KBXINST','keys param-1 '//cfi(param_in(1))) call mctrace(1,'KBXINST','keys param-2 '//cfi(param_in(2))) call mctrace(1,'KBXINST','keys param-3 '//cfi(param_in(3))) call mctrace(1,'KBXINST','keys param-4 '//cfi(param_in(4))) if (mod(ssnum,2) .eq. 0) then if (band .lt. 19) then vis_data = .FALSE. ir_data = .TRUE. else vis_data = .TRUE. ir_data = .FALSE. endif else if (band .gt. 1) then vis_data = .FALSE. ir_data = .TRUE. else vis_data = .TRUE. ir_data = .FALSE. endif endif call mctrace(1,'KBXINST','check band '//cfi(band)) if (vis_data) then param_out(1) = 4 param_out(2) = lit('RAW ') param_out(3) = lit('RAD ') param_out(4) = lit('ALB') param_out(5) = lit('BRIT') elseif (ir_data) then param_out(1) = 4 param_out(2) = lit('RAW ') param_out(3) = lit('RAD ') param_out(4) = lit('TEMP') param_out(5) = lit('BRIT') endif c c --- Check for a stretch table (SU) c if (ischar(param_in(38)) .eq. 1) then call movwc(param_in(38),su_file) if (brkset(su_file,cal_input) .ne. 0) then kbxopt = -3 endif endif elseif (option .eq. 'BRKP') then c c --- param_in(1) is the su table name c call movwc(param_in(1),su_file) if (brkset(su_file,cal_input) .ne. 0) then kbxopt = -3 endif elseif (option .eq. 'INFO') then c c --- param_in contains frame directory - set ir_data/vis_data flag c band = param_in(4) call mctrace(1,'KBXINST','info param-1 '//cfi(param_in(1))) call mctrace(1,'KBXINST','info param-2 '//cfi(param_in(2))) call mctrace(1,'KBXINST','info param-3 '//cfi(param_in(3))) call mctrace(1,'KBXINST','info param-4 '//cfi(param_in(4))) if (mod(ssnum,2) .eq. 0) then if (band .lt. 19) then vis_data = .FALSE. ir_data = .TRUE. else vis_data = .TRUE. ir_data = .FALSE. endif else if (band .gt. 1) then vis_data = .FALSE. ir_data = .TRUE. else vis_data = .TRUE. ir_data = .FALSE. endif endif if (vis_data) then param_out(1) = 4 param_out(2) = lit('RAW ') param_out(3) = lit('RAD ') param_out(4) = lit('ALB') param_out(5) = lit('BRIT') param_out(6) = lit(' ') param_out(7) = lit('1/SR') param_out(8) = lit('* ') param_out(9) = lit(' ') param_out(10) = 1 param_out(11) = 10000 param_out(12) = 10 param_out(13) = 1 elseif (ir_data) then param_out(1) = 4 param_out(2) = lit('RAW ') param_out(3) = lit('RAD ') param_out(4) = lit('TEMP') param_out(5) = lit('BRIT') param_out(6) = lit(' ') param_out(7) = lit('MW**') param_out(8) = lit(' K ') param_out(9) = lit(' ') param_out(10)= 1 param_out(11)= 10000 param_out(12)= 100 param_out(13)= 1 endif c --- param_in contains calibration block from server elseif (option .eq. 'CALB') then cal_flag = 1 call movw( 171, param_in, cal_block ) else kbxopt = -1 endif return end C$ Name: C$ rad_to_temp - Converts radiance to temperature for INSAT-3D calibration C$ C$ Interface: C$ real function C$ rad_to_temp(real radiance, integer band, isat) C$ rad_to_temp(real radiance, real wavelength) C$ C$ Input: C$ radiance - radiance C$ wavelength - central wavelength C$ C$ Return values: C$ C$ temperature - temperature - units of Kelvin C$ C$ Remarks: C$ C$ Categories: C$ calibration real function rad_to_temp(radiance,central_wavelength) implicit none c c - symbolic constants & shared data c integer NUM_CONST parameter (NUM_CONST=2) c c --- input parameters c real radiance ! radiance real central_wavelength ! radiance integer band ! band number integer isat ! satellite ss number c c --- local variables c real wavelength ! central wavelength real temperature ! temperature value returned real adjusted_rad ! radiance adjusted by derived constants real boltzman_constant real central_wavenumber real cwl real cwn real constant1 ! constant using plancks constant and speed of light real constant2 ! constant using plancks constant, speed of light and boltzman constant real speed_of_light real planck_constant data boltzman_constant /1.380658e-23/ data planck_constant /6.6260755e-34/ data speed_of_light /2.9979246e+8/ c c --- Converting from radiance to temperature base on INSAT-3D documentation c cwl = central_wavelength/1000000.0 central_wavenumber = central_wavelength/1000000.0 cwn = central_wavenumber constant1 = 2.0 * planck_constant * speed_of_light**2 constant2 = (planck_constant * speed_of_light) / boltzman_constant temperature = constant2 / (cwl * & log(constant1 / (1.0e6 * & radiance * 10.0 * & cwn**5.0) + 1.0)) rad_to_temp = temperature return end