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: kbxwari.dlm,v 1.8 2015/06/18 18:14:31 russd Tst $ *** C$ Name: C$ kbxini - Calibration module for HIMAWARI data 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$ 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 --- HIMAWARI 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(512) ! cal block common/himawari/ & cal_input, & cal_output, & src_byte_size, & dest_byte_size, & cal_flag, & cal_block c --- set local variables to fill common block cal_input = input_cal cal_output = output_cal c --- Verify that the input calibration type if (cal_input(1:3) .ne. 'RAW') then kbxini = -1 c --- Verify that the output calibration type 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:3) .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 c --- constants integer MAXTBLSIZE parameter (MAXTBLSIZE = 65536) c --- input parameters 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*2 ibuf(*) ! I/O array containing pixels to be modified c --- external functions real temp_to_rad ! converts temperature to radiance character*12 cfi character*12 cff character*12 cfe character*4 clit c --- local variables character*12 cval integer area_number ! area number found in area directory integer ss_value integer cal_offset integer i ! do loop index integer j ! do loop index integer last_area ! last area used for araget integer last_band ! last area used for araget integer iRAW integer iRAD integer iTEMP integer iALB integer iBRIT integer init, ival integer jbuf(11000) integer*2 jbuf2(22000) real H ! Planck Constant real C ! Speed of Light real K ! Boltzmann Constant real rALB double precision dRAW double precision dRAD double precision dTEMP double precision dALB double precision RAW_SCALE double precision RAD_SCALE double precision TEMP_SCALE double precision ALB_SCALE logical ir_data ! flag indicating ir data requested logical vis_data ! flag indicating visible data requested c --- HIMAWARI 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(512) ! cal block - includes directory and tables common/himawari/ & cal_input, & cal_output, & src_byte_size, & dest_byte_size, & cal_flag, & cal_block c --- JMA Cal common integer calb_BlockLen integer calb_bandNo integer calb_bitPix integer calb_errorCount integer calb_outCount integer calb_invalidValue double precision calb_waveLen double precision calb_gainCnt2rad double precision calb_cnstCnt2rad double precision calb_rad2btpC0 double precision calb_rad2btpC1 double precision calb_rad2btpC2 double precision calb_btp2radC0 double precision calb_btp2radC1 double precision calb_btp2radC2 double precision calb_lightSpeed double precision calb_planckConst double precision calb_bolzConst double precision calb_rad2albedo common/JMA/ & calb_BlockLen, & calb_bandNo, & calb_bitPix, & calb_errorCount, & calb_outCount, & calb_invalidValue, & calb_waveLen, & calb_gainCnt2rad, & calb_cnstCnt2rad, & calb_rad2btpC0, & calb_rad2btpC1, & calb_rad2btpC2, & calb_btp2radC0, & calb_btp2radC1, & calb_btp2radC2, & calb_lightSpeed, & calb_planckConst, & calb_bolzConst, & calb_rad2albedo c --- HIMAWARI calibration table integer table(MAXTBLSIZE) common/CALTBL/ & table data last_area /-1/ data last_band /-1/ data RAW_SCALE /1.0000D0/ data RAD_SCALE /1000.0D0/ data TEMP_SCALE/100.0D0/ data ALB_SCALE /100.0D0/ data init/0/ equivalence (jbuf,jbuf2) c --- module return code kbxcal = 0 c --- set ir_data/vis_data flag if (band .le. 6) then vis_data = .TRUE. ir_data = .FALSE. else vis_data = .FALSE. ir_data = .TRUE. endif c --- read in calibration block area_number = area_dir(33) ss_value = area_dir(3) cal_offset = area_dir(63) if (last_area .ne. area_number) then call mctrace(1,'KBXWARI','CAL Function') if( cal_flag.eq.0 ) then call araget(area_number,cal_offset,512,cal_block) endif c ------ set the constants H = 6.6260755e-34 C = 2.9979246e+8 K = 1.380658e-23 c ------ values from CAL block go into JMA cal common calb_BlockLen = cal_block(2) calb_bandNo = cal_block(3) calb_bitPix = cal_block(4) calb_errorCount = cal_block(5) calb_outCount = cal_block(6) calb_invalidValue = cal_block(7) calb_waveLen = DBLE(cal_block(8)) *0.0000001D0 calb_gainCnt2rad = DBLE(cal_block(9)) *0.0000001D0 calb_cnstCnt2rad = DBLE(cal_block(10))*0.0000001D0 calb_rad2btpC0 = DBLE(cal_block(11))*0.0000001D0 calb_rad2btpC1 = DBLE(cal_block(12))*0.0000001D0 calb_rad2btpC2 = DBLE(cal_block(13))*0.000000000001D0 calb_btp2radC0 = DBLE(cal_block(14))*0.0000001D0 calb_btp2radC1 = DBLE(cal_block(15))*0.0000001D0 calb_btp2radC2 = DBLE(cal_block(16))*0.0000001D0 calb_rad2albedo = DBLE(cal_block(17))*0.0000001D0 calb_lightSpeed = DBLE(C) calb_planckConst = DBLE(H) calb_bolzConst = DBLE(K) cval = cfe( REAL(calb_waveLen), 6 ) call mctrace(1,'KBXWARI','WaveLength ='//cval) cval = cfe( REAL(calb_rad2btpC0), 6 ) call mctrace(1,'KBXWARI','RadtoBtC0 ='//cval) cval = cfe( REAL(calb_rad2btpC1), 6 ) call mctrace(1,'KBXWARI','RadtoBtC1 ='//cval) cval = cfe( REAL(calb_rad2btpC2), 6 ) call mctrace(1,'KBXWARI','RadtoBtC2 ='//cval) cval = cfe( REAL(calb_lightspeed), 6 ) call mctrace(1,'KBXWARI','LightSpeed ='//cval) cval = cfe( REAL(calb_planckConst), 6 ) call mctrace(1,'KBXWARI','Planck ='//cval) cval = cfe( REAL(calb_bolzConst), 6 ) call mctrace(1,'KBXWARI','Boltzmann ='//cval) c ------ make the lookup table call maktbl( band ) endif c --- loop through incoming RAW counts do i = 1, num_pixels c ------ use lookup table to calibrate values jbuf(i) = table(ibuf(i)+1) enddo c --- move computed values into request byte size call mpixel( & num_pixels, & 4, & dest_byte_size, & jbuf & ) c --- move data into IO array if( dest_byte_size.eq.1) then call movw( num_pixels/4, jbuf, ibuf ) elseif( dest_byte_size.eq.2) then call movw( num_pixels/2, jbuf, ibuf ) else call movw( num_pixels, jbuf, ibuf ) endif c --- save last area and band numbers last_area = area_number 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 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 logical ir_data ! flag indicating ir data requested logical vis_data ! flag indicating visible data requested c --- HIMAWARI 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(512) ! cal block - includes directory and tables common/himawari/ & cal_input, & cal_output, & src_byte_size, & dest_byte_size, & cal_flag, & cal_block kbxopt = 0 if (option .eq. 'KEYS') then c c --- param_in contains frame directory - set ir_data/vis_data flag c band = param_in(4) if (band .le. 6) then vis_data = .TRUE. ir_data = .FALSE. else vis_data = .FALSE. ir_data = .TRUE. 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') 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) if (band .le. 6) then vis_data = .TRUE. ir_data = .FALSE. else vis_data = .FALSE. ir_data = .TRUE. 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('WM**') param_out(8) = lit(' % ') param_out(9) = lit(' ') param_out(10) = 1 param_out(11) = 1000 param_out(12) = 100 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('WM**') param_out(8) = lit(' K ') param_out(9) = lit(' ') param_out(10)= 1 param_out(11)= 1000 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( 6528, param_in, cal_block ) else kbxopt = -1 endif return end subroutine RADtoTEMP( & RAD, & TEMP & ) implicit none c --- parameters double precision RAD double precision TEMP c --- external functions character*12 cfe character*12 cff c --- internal variables character*12 cval double precision radiance double precision effective_temperature double precision lambda double precision planck_c1 double precision planck_c2 integer init integer Rinit c --- JMA Cal common integer calb_BlockLen integer calb_bandNo integer calb_bitPix integer calb_errorCount integer calb_outCount integer calb_invalidValue double precision calb_waveLen double precision calb_gainCnt2rad double precision calb_cnstCnt2rad double precision calb_rad2btpC0 double precision calb_rad2btpC1 double precision calb_rad2btpC2 double precision calb_btp2radC0 double precision calb_btp2radC1 double precision calb_btp2radC2 double precision calb_lightSpeed double precision calb_planckConst double precision calb_bolzConst double precision calb_rad2albedo common/JMA/ & calb_BlockLen, & calb_bandNo, & calb_bitPix, & calb_errorCount, & calb_outCount, & calb_invalidValue, & calb_waveLen, & calb_gainCnt2rad, & calb_cnstCnt2rad, & calb_rad2btpC0, & calb_rad2btpC1, & calb_rad2btpC2, & calb_btp2radC0, & calb_btp2radC1, & calb_btp2radC2, & calb_lightSpeed, & calb_planckConst, & calb_bolzConst, & calb_rad2albedo data init/0/ data Rinit/0/ c --- initialize constants if( init.eq.0 ) then c ------ central wave length lambda = calb_waveLen/1000000.0 c ------ planck_c1 = (2 * h * c^2 / lambda^5) planck_c1 = 2.0 & * calb_planckConst & * (calb_lightSpeed**2) & / (lambda**5) c ------ planck_c2 = (h * c / k / lambda ) planck_c2 = calb_planckConst & * calb_lightSpeed & / calb_bolzConst & / lambda init = 1 endif Rinit = 1 if( RAD.ge.5.850D0 .and. & RAD.le.5.860D0 & ) Rinit = 0 c --- radiance radiance = RAD*1000000.0 c --- compute the brightness temperature TEMP = -9999.0D0; if( radiance.gt.0.0D0 ) then effective_temperature = planck_c2 / & DLOG( (planck_c1 / radiance ) + 1.0 ) TEMP = calb_rad2btpC0 & + calb_rad2btpC1 & * effective_temperature & + calb_rad2btpC2 & * (effective_temperature**2) if( Rinit.eq.0 ) then cval = cff( RAD, 6 ) call mctrace(1,'KBXWARI','RAD='//cval) cval = cff( radiance, 2 ) call mctrace(1,'KBXWARI','radiance='//cval) cval = cfe( calb_rad2btpC0, 6 ) call mctrace(1,'KBXWARI','rad2btpC0='//cval) cval = cfe( calb_rad2btpC1, 6 ) call mctrace(1,'KBXWARI','rad2btpC1='//cval) cval = cfe( calb_rad2btpC2, 6 ) call mctrace(1,'KBXWARI','rad2btpC2='//cval) cval = cfe( planck_c1, 6 ) call mctrace(1,'KBXWARI','Planck C1='//cval) cval = cfe( planck_c2, 6 ) call mctrace(1,'KBXWARI','Planck C2='//cval) cval = cff( effective_temperature, 2 ) call mctrace(1,'KBXWARI','Eff TEMP='//cval) cval = cff( TEMP, 6 ) call mctrace(1,'KBXWARI','TEMP='//cval) Rinit = 1 endif endif return end subroutine maktbl( band ) implicit none c --- constants integer MAXTBLSIZE parameter (MAXTBLSIZE = 65536) c --- parameters integer band c --- external functions c --- internal variables integer last_area ! last area used for araget integer last_band ! last area used for araget integer i integer iCNT integer iRAW integer iRAD integer iTEMP integer iALB integer iBRIT integer init, ival double precision RAW_SCALE double precision RAD_SCALE double precision TEMP_SCALE double precision ALB_SCALE double precision dRAW double precision dRAD double precision dTEMP double precision dALB logical ir_data ! flag indicating ir data requested logical vis_data ! flag indicating visible data requested real H ! Planck Constant real C ! Speed of Light real K ! Boltzmann Constant real rALB c --- JMA Cal common integer calb_BlockLen integer calb_bandNo integer calb_bitPix integer calb_errorCount integer calb_outCount integer calb_invalidValue double precision calb_waveLen double precision calb_gainCnt2rad double precision calb_cnstCnt2rad double precision calb_rad2btpC0 double precision calb_rad2btpC1 double precision calb_rad2btpC2 double precision calb_btp2radC0 double precision calb_btp2radC1 double precision calb_btp2radC2 double precision calb_lightSpeed double precision calb_planckConst double precision calb_bolzConst double precision calb_rad2albedo common/JMA/ & calb_BlockLen, & calb_bandNo, & calb_bitPix, & calb_errorCount, & calb_outCount, & calb_invalidValue, & calb_waveLen, & calb_gainCnt2rad, & calb_cnstCnt2rad, & calb_rad2btpC0, & calb_rad2btpC1, & calb_rad2btpC2, & calb_btp2radC0, & calb_btp2radC1, & calb_btp2radC2, & calb_lightSpeed, & calb_planckConst, & calb_bolzConst, & calb_rad2albedo c --- HIMAWARI 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(512) ! cal block common/himawari/ & cal_input, & cal_output, & src_byte_size, & dest_byte_size, & cal_flag, & cal_block c --- HIMAWARI calibration table integer table(MAXTBLSIZE) common/CALTBL/ & table data RAW_SCALE /1.0000D0/ data RAD_SCALE /1000.0D0/ data TEMP_SCALE/100.0D0/ data ALB_SCALE /100.0D0/ c --- set ir_data/vis_data flag if (band .le. 6) then vis_data = .TRUE. ir_data = .FALSE. else vis_data = .FALSE. ir_data = .TRUE. endif c --- loop through incoming RAW counts do i = 1, MAXTBLSIZE c --- calibrate the RAW units if (cal_output(1:3) .eq. 'RAW') then dRAW = dble(i-1) iRAW = INT( dRAW*RAW_SCALE ) c --- need to convert to another physical value else c ------ count value is table index iCNT = i-1 c ------ convert RAW to RAD if( & iCNT.eq.calb_errorCount .or. & iCNT.eq.calb_outCount .or. & iCNT.eq.0 & ) then dRAD = calb_invalidValue iRAD = 0 else dRAD = dble(iCNT)*calb_gainCnt2rad+calb_cnstCnt2rad if( dRAD.le.0.0D0 ) dRAD = calb_invalidValue iRAD = INT( dRAD*RAD_SCALE ) endif c ------ for reflectance bands if( vis_data ) then c --------- convert RAD to ALB/BRIT if( & cal_output(1:3).eq.'ALB' .or. & cal_output(1:4).eq.'BRIT' & ) then if( dRAD.eq.calb_invalidValue ) then dALB = 0 iBRIT = 0 else dALB = calb_rad2albedo*dRAD rALB = REAL(dALB) iBRIT = NINT(SQRT(rALB)*255.0) endif iALB = INT( dALB*100.0D0*ALB_SCALE ) iBRIT = MIN(MAX(iBRIT,0),255) endif c ------ for radiance bands else c --------- convert RAD to TEMP/BRIT if( & cal_output(1:4).eq.'TEMP' .or. & cal_output(1:4).eq.'BRIT' & ) then if( dRAD.eq.calb_invalidValue ) Then dTEMP = calb_invalidValue iBRIT = 255 else call RADtoTEMP(dRAD, dTEMP) call gryscl(REAL(dTEMP),iBRIT) endif iTEMP = INT( dTEMP*TEMP_SCALE ) iBRIT = MIN(MAX(iBRIT,0),255) endif endif endif c --- move the computed value into the array if( cal_output(1:3).eq.'RAW' ) table(i) = iRAW if( cal_output(1:3).eq.'RAD' ) table(i) = iRAD if( cal_output(1:3).eq.'TEM' ) table(i) = iTEMP if( cal_output(1:3).eq.'ALB' ) table(i) = iALB if( cal_output(1:3).eq.'BRI' ) table(i) = iBRIT enddo return end