subroutine main0 c c PURPOSE: c UW/CIMSS Regression GOES SST. c c Full navigation. c c METHOD: c Regression algorithms for GOES-8/9 are based on the GOES/BUOY c collocations from September 1996 to August 1997. Regression c algorithm for GOES-10 is based on GOES-8/GOES-10 collocations c in March-April 1998 (GOES-10 science check up), using GOES-8 c SST as truth. These are: c c G08_DAY=-12.10+1.0472T11+0.2634(T11-T12)+0.4284(T11-T12)^2+1.0984ANGLE c G08_NIT= 0.34+1.0095T39+0.1123(T39-T11)+ +0.8553ANGLE c G10_DAY=-20.13+1.0707T11+1.7041(T11-T12)+0.1407(T11-T12)^2+1.1909ANGLE c G10_NIT=-27.74+1.1010T39+0.3933(T39-T11)+ +3.3412ANGLE c c where T39, T11, and T12 are brightness temperature at 3.9, 11, and c 12 micron channels, respectively, and ANGLE is sec(satellite viewn c angle)-1. c c I/O FILES: c In the following table, # is 1-4 for GOES-E NH, GEOS-W NH, c GOES-E SH, and GOES-W SH, respectively. c #00-#23: Guess SST for the hour c #30-#32: Radiances (half-hour before, current, and half-hour after) c #34: Land mask c #35: Solar zenith angle c #40-#42: SST (half-hour before, current, and half-hour after) c c HISTORY: c 98/09/28 New SST formula. Correction for vis cal. 0K>> Regression coefficients <<< real c1(5,2) ! WU style retr coeffs real c2(8,2) ! Edinburgh style retr coeffs c >>> Size limits for GOES image <<< integer max_line, max_elem parameter (max_line=1900, max_elem=3500) c >>> Albedo thresholds for cloud <<< integer dark, bright, glint parameter (dark=40, bright=dark*2, glint=dark*4) real bd_sq parameter (bd_sq=1./((bright-dark)*(bright-dark))) c >>> More parameters for cloud detection <<< integer twilight, off_limit real small parameter (twilight=80, off_limit=2000, small=1.e-6) c >>> Code for invalid SST <<< integer space, land, cloud parameter (space=0, land=2, cloud=4) c ------------ c Major arrays c ------------ c MJW 06/25/02 integer*1 scaled_sst(max_elem) ! Out SST, (K-270)/.15 c MJW integer*2 new_guess(max_elem) ! New guess SST, K*10 integer*2 t1(max_elem,max_line) ! Albedo, %*1000 integer*2 t2(max_elem,max_line) ! Tb(3.9), K*10 integer*2 t4(max_elem,max_line) ! Tb(11), K*10 integer*2 t5(max_elem,max_line) ! Tb(12), K*10 integer*2 t_bef(max_elem,max_line) ! Earlier Tb(11), K*10 integer*2 t_aft(max_elem,max_line) ! Later Tb(11), K*10 integer*2 guess(max_elem,max_line) ! Guess SST, K*10 integer*2 retrv(max_elem,max_line) ! Retrieved SST, K*10 integer*2 incrm(max_elem,max_line) ! Increment for guess equivalence (incrm(1,1),t1(1,1)) ! Save some memory c SNM real env_dat(432) ! BT threshold envelope data real rad_bt_dat(4096) ! Radiance to BT LUT c SNM c --------------- c Other variables c --------------- c >>> For handling McIDAS AREA-file <<< integer hd(64), param(640) character*32 comment equivalence (hd(25),comment) c >>> Input, guess, and output images <<< integer image, g_sst, jmage c >>> Results from navigation <<< integer scene, sat_z, sun_z c >>> Loop variables and bounds <<< integer line, elem, ii, jj, kk c SNM integer*2 kk2 !SNM temp variable for 2.2.3 c SNM integer beg_line, end_line, beg_elem, end_elem c >>> Ealier/later images available? <<< logical has_prev_img, has_next_img c >>> Used in processing SST <<< real mean, vari, t4_t5, t2_t4, score, xx c >>> Used in updating guess <<< real w(-5:5,-5:5), weighted_sum, total_weight integer ikwp, iretr_type c SNM c >>> variables for env_check function <<< integer bt_3, bt_11, bt_12 ! Temp store for BTs c >> to allow passing to function << c >>> channel combination, results of env_check <<< integer chan, env_temp1, env_temp2 integer sat, flag ! satellite #, result of env_check integer tol ! tolerance for env_check parameter (tol=0) c >>> vairiables for Sunglint <<< integer maskglint ! C function integer jday, hour, minute c >>> sat, sun zenith in radians <<< real sat_z_rad, sun_z_rad real a2r ! deg to rad conversion factor parameter (a2r=3.14159265/180) real windspeed real glint_rad ! glint radiance real dbt ! difference BT integer dsst_res ! result of mask_glint integer aa ! counter for reading in arrays integer speed ! 0 = use wind velocity parameter (speed = 1) ! 1 = use wind speed real orglat, orglon ! satellite lat and lon real pxlat, pxlon ! pixel lat and lon c SNM c ============== c 1. PREPARATION c ============== c >>> Decide WU vs EDINBURGH retrievals <<< iretr_type = ikwp('COEF_TYPE',1,-1) if ( iretr_type .eq. 1 ) then ! Fred Wu / orignial retrieval type write (6,*) ' Using old style (WU) retrievals' else if ( iretr_type .eq. 2 ) then ! Edinburgh retrieval model write (6,*) ' Using new style (Edinburgh) retrievals' else write (6,*) ' No Retrieval Style Specified ' write (6,*) ' Must Specify COEF_TYPE= 1 or 2 for WU or New' write (6,*) ' Usage: goes_sst.k FIRST_IMAGE=<> COEF_TYPE=<>' stop 'ERROR in deciding Retrieval Type' endif c c +++++++++++++++++ c 1.1. Observations c +++++++++++++++++ c c ------------------------------ c 1.1.1. Check for Current Image c ------------------------------ call gettim(ii) print*, ii, ' SST 0005', image call readd(image,hd) ! Read the header if (hd(1) .ne. 0) then write (*,'("GOES image not available - next image")') return ! cycle end if c ----------------------------------- c 1.1.2. Get the Current Measurements c ----------------------------------- c >>> Open and set calibration <<< call opnara(image) ! Open the file call araopt(image,1,lit('SPAC'),2) ! Read as I*2 and ... call araopt(image,1,'UNIT',lit('ALB ')) ! ... calibrated as albedo c >>> Read with quality control <<< c Valid range [0.1%,99.9%], otherwise -1. do line=1,hd(9) call read_area(image,line-1,0,hd(10),1,t1(1,line),1,999,-1,-2) end do c >>> Re-set calibration <<< call araopt(image,1,'UNIT',lit('TEMP')) ! ... calibrated as Tb c >>> Read with quality control <<< c Valid range [150K,320K], otherwise -1. do line=1,hd(9) call read_area(image,line-1,0,hd(10),2,t2(1,line), & 1500,3200,-1,-2) call read_area(image,line-1,0,hd(10),4,t4(1,line), & 1500,3200,-1,-2) call read_area(image,line-1,0,hd(10),5,t5(1,line), & 1500,3200,-1,-2) end do call clsara(image) c --------------------- c 1.1.3. Previous Image c --------------------- c >>> Search for an earlier image <<< c Use one hour earlier image if available. Otherwise, try c half hour earlier. If neither available, take penulty later. do jmage=image-2,image call readd(jmage,hd) if (hd(1) .eq. 0) goto 115 ! exit - this image is available end do 115 continue if (jmage .ne. image) then c >>> An earlier image is found <<< call opnara(jmage) call araopt(jmage,1,lit('SPAC'),2) call araopt(jmage,1,'UNIT',lit('TEMP')) do line=1,hd(9) c Read with QC: valid range [150K,320K], otherwise -1. call read_area(jmage,line-1,0,hd(10),4,t_bef(1,line), & 1500,3200,-1,-2) end do call clsara(jmage) has_prev_img = .true. else c >>> An earlier image is not available <<< has_prev_img = .false. end if c ----------------- c 1.1.4. Next Image c ----------------- c >>> Search for a later image <<< c Use one hour later image if available. Otherwise, try c half hour later. If neither available, take penulty later. do jmage=image+2,image,-1 call readd(jmage,hd) if (hd(1) .eq. 0) goto 116 ! exit - this image is available end do 116 continue if (jmage .ne. image) then c >>> A later image is found <<< call opnara(jmage) call araopt(jmage,1,lit('SPAC'),2) call araopt(jmage,1,'UNIT',lit('TEMP')) do line=1,hd(9) c Read with QC: valid range [150K,320K], otherwise -1. call read_area(jmage,line-1,0,hd(10),4,t_aft(1,line), & 1500,3200,-1,-2) end do call clsara(jmage) has_next_img = .true. else c >>> An later image is not available <<< has_next_img = .false. end if c +++++++++++++++++++ c 1.2. Auxillary Data c +++++++++++++++++++ c c ---------------- c 1.2.1. Guess SST c ---------------- c Considering the diurnal variation of SST, it seems natural c to use a guess specific for the hour. c c But that has the disadvantage of guess being updated c infrequently, particularly for night hours, since SST c is retrieved each time for only 5-15% of the pixels. c c One compromise is to use one guess for all hours while c relax the SST-guess test somehow. c c Alternatively, we can increase the "influence radius" so c that retrieval is used to update a larger number of guess. c c We are experimenting with both strategies and will decide c later which one to be adopted. c c >>> Guess for each hour <<< g_sst = hd(5)/10000 + nint(mod(hd(5),10000)/6.e3) ! Hour g_sst = image - 31 + mod(g_sst,24) ! guess SST for the hour c MJW 06/26/02 cccccc g_sst = g_sst-3000 ! The New Regime - guess files <1000 c MJW c >>> One guess for all hours <<< c g_sst = image - 3 ! Guess SST for the hour c >>> Read guess SST << call readd(g_sst,hd) ! Read the header call opnara(g_sst) call araopt(g_sst,1,lit('SPAC'),2) call araopt(g_sst,1,'UNIT',lit('BRIT')) do line=1,hd(9) c Read with QC: valid range [270K,310K], otherwise stop. call read_area(g_sst,line-1,0,hd(10),1,guess(1,line), & 2700,3100,0,4) end do c >>> Update navigation of guess SST <<< c There may be a need to update the navigation of guess SST c for the sake of its displaying it in McIDAS. The following c doesn't seem to work, so no update is implemented at the c moment. This, however, should affect only the appearance c when displayed. c call clsara(g_sst) c call makara(g_sst,hd) c call opnara(g_sst) c kk = hd(34) - hd(35) c if (hd(63).ne.0) kk=hd(63)-hd(35) c call araget (image,hd(35),kk,param) c call araput (g_sst,hd(35),kk,param) c ------------------------------ c 1.2.2. Regression coefficients c ------------------------------ c c >>> Decide WU vs EDINBURGH retrievals <<< ccc iretr_type = ikwp('COEF_TYPE',1,-1) if ( iretr_type .eq. 1 ) then ! Fred Wu / orignial retrieval type c >>> Open coefficient file <<< open (66, file='../data/reg_coeff.asc', status='old', form='formatted', & access='sequential') c >>> Searching for matching sat ID <<< kk = 0 do while (hd(3) .ne. kk) read (66,*) kk c >>> Error if hit EOF <<< if (kk .gt. 99) then write (*,'("Unknown satellite",2i15)') hd(3), kk stop 'Unknown satellite - will stop' c >>> Read coefficients <<< else read (66,*) ((c1(ii,jj),ii=1,5),jj=1,2) end if end do close (66) c >>> Adjust coefficients <<< c ... because we are dealing with Tb*10 do jj=1,2 c1(1,jj) = c1(1,jj)*10 ! coeff. for const. term c1(4,jj) = c1(4,jj)*.1 ! coeff. for dt*dt c c1(5,jj) = c1(5,jj)*1.e-3 ! coeff. for [sec(theta)-1]*1000 c1(5,jj) = c1(5,jj)*1.e-2 ! SNM accounts for Tb*10 factor end do print*, ((c1(ii,jj),ii=1,5),jj=1,2) ccc SNM else if ( iretr_type .eq. 2 ) then ! Edinburgh retrieval model c >>> open new coefficients file <<< open (66, file='../data/reg_coeff_mer.asc', status='old') c >>> searching for matching sat ID <<< kk = 0 do while (hd(3) .ne. kk) read (66,*) kk c >>> Error if hit EOF <<< if (kk .gt. 99) then write (*,'("Unknown satellite",2i15)') hd(3), kk stop 'Unknown satellite - will abort' c >>> Read coefficients <<< else read (66,*) ((c2(ii,jj),ii=1,8),jj=1,2) endif end do close (66) c >>> adjust coefficients <<< c .... because we are dealing with Tb*10 do jj=1,2 c2(1,jj) = c2(1,jj)*10 ! coeff. for constant term c2(2,jj) = c2(2,jj)*1.e-2 ! coeff. for [sec(theta)-1]*1000 c2(4,jj) = c2(4,jj)*1.e-3 ! coeff. for [[sec(theta)-1]*1000]*Tb*10 c2(6,jj) = c2(6,jj)*1.e-3 ! coeff. for ["] c2(8,jj) = c2(8,jj)*1.e-3 ! coeff. for ["] end do print*, ((c2(ii,jj),ii=1,8),jj=1,2) ccc SNM else ! NO RETRIEVAL STYLE@#$%^&*( write (6,*) ' No Retrieval Style Specified ' write (6,*) ' Must Specify COEF_TYPE= 1 or 2 for WU or New' stop 'ERROR in deciding Retrieval Type' endif c +++++++++++++++++++++++++++ c 1.3. Set Up the Output File c +++++++++++++++++++++++++++ c >>> Prepare the area file <<< call readd(image,hd) ! Read the header c MJW 06/27/02 jmage = image + 11 ! retrieved SST ccccc jmage = image + 1011 ! place in 4000 range c MJW hd(11) = 1 ! bytes per pixel hd(14) = 1 ! one band only hd(19) = 1 ! data are in band 1 hd(36) = 0 comment= 'GOES SST, UW-CIMSS, 00/05/01' hd(49) = 0 hd(50) = 0 hd(51) = 0 ccc hd(52) = 'PRD ' ! calibration codde hd(52) = lit('PRD ') ! calibration codde ccc hd(53) = 'BRIT' hd(53) = lit('BRIT') call makara(jmage,hd) call opnara(jmage) kk = hd(34) - hd(35) if (hd(63).ne.0) kk=hd(63)-hd(35) call araget (image,hd(35),kk,param) call araput (jmage,hd(35),kk,param) c ============= c 2. PROCESSING c ============= c SNM >>> Read in data for Env_check and Sunglint tests <<< c >>> get julian day, hour and minute from file header <<< jday = mod(hd(4),1000) hour = int(hd(5)/10000) minute = int(hd(5)/100) minute = mod(minute,100) c >>> initialise windspeed <<< windspeed = 10.0 c >>> open BT threshold data file and read in array <<< sat = (hd(3)-70)/4 open (20, file='../data/env_thresh_dat.txt') read(20,*) (env_dat(aa),aa=1,432) c >>> open radiance to BT LUT file and read in array <<< if (hd(3) .eq. 70) then open(21, file='../data/goes8_radbt_edit.txt') read(21,*) (rad_bt_dat(aa),aa=1,4096) else if (hd(3) .gt. 70) then open(21, file='../data/goes10_radbt_edit.txt') read(21,*) (rad_bt_dat(aa),aa=1,4096) endif c SNM c --------------------------------------------- c 2.0.1. Pick a pixel and define its 3-by-3 box c --------------------------------------------- do line=1,hd(9) beg_line = max(line-1,1) end_line = min(line+1,hd(9)) do elem=1,hd(10) beg_elem = max(elem-1,1) end_elem = min(elem+1,hd(10)) c ---------------------------- c 2.0.2. Initialize test score c ---------------------------- score = 1 c SNM ------------------------------- c 2.0.3. Initialize SST retrieval c ------------------------------- retrv(elem,line) = 5 c SNM c ++++++++++++++++++++++ c 2.1. General Screening c ++++++++++++++++++++++ c c 00/06/04 ORA1 can't handle full navigation. c 00/09/29 Full navigation is now tried, with revised nav module c c The first call to navPixab2 computes only the lat/lon and c land/sea flag. This part of the computation now calls navline, c which navigates a who;e image line at once. The second call c to navPixab2 computes satellite and sun angles. c ----------------- c 2.1.1. Navigation c ----------------- c Navigating latlon is moved back here before cloudclear, c since a whole line is now computed at once. c The sun angle computations are postponed after Tb(11) test c of 2.1.3. to avoid this extra calculation for "obvious" c cloudy pixels that count almost half of all pixels. That c speeds up the processing c ccc call navpixab2(image,twilight,line,elem, scene,sat_z,sun_z) ccc call nav_pixab2(image,twilight,line,elem, scene,sat_z,sun_z) c SNM 28/01/02 c >>> Call to updated nav_pixab2 function <<< call nav_pixab3(image,twilight,line,elem, scene,sat_z,sun_z, & orglat, orglon, pxlat, pxlon) c SNM c ------------------------------------ c 2.1.2. Test on Topography (land/sea) c ------------------------------------ if (scene .lt. 4) then if (scene .eq. 0) then ! Space retrv(elem,line) = space else ! Land, Partial/Ephemeral Water retrv(elem,line) = land end if goto 10 ! cycle - next pixel end if c c ---------------------- c 2.1.3. Test for Tb(11) c ---------------------- c a) Even though SST>Tb(11) and guess could have error, c Tb(11)>(guess-10) should be sufficient evidence of cloud. c b) If Tb(11)<270K, there is little chance for c SST>271.05K, the lower limit of the accepted SST value. c if (t4(elem,line).lt.max(guess(elem,line)-100,2700)) then retrv(elem,line) = cloud ! cloud/frozen goto 10 ! cycle - next pixel end if c -------------------------------- c 2.1.4. Test for Tb(3.9) & Tb(12) c -------------------------------- c Tb(3.9) and Tb(12) are both used for cloud detection and c one of them is used for atmospheric correction. c c This test is placed after Tb(11) test so highly reflective c clouds, which cause "invalid" Tb(3.9), will not be c identified as "space" (i.e., no data). c c SNM 02/02/01 c >>> extra choice added to remove t5 check for G12 <<< if (hd(3) .le. 76) then if (t2(elem,line).lt.0 .or. t5(elem,line).lt.0) then retrv(elem,line) = cloud goto 10 end if else if (hd(3) .ge. 78) then if (t2(elem,line).lt.0) then retrv(elem,line) = cloud goto 10 end if end if c SNM c --------------------- c 2.1.4.5 Sun position c --------------------- c Now need sun angles. Postponed 'til here, because no sense c doing this calculation for space,land,or clouds c c >>> Sun position <<< ccc call navpixab2(image,twilight,line,elem, scene,sat_z,sun_z) ccc call nav_pixab2(image,twilight,line,elem, scene,sat_z,sun_z) c SNM c >>> call to updated nav_pixab2 function <<< call nav_pixab3(image,twilight,line,elem, scene, sat_z,sun_z, & orglat, orglon, pxlat, pxlon) c SNM if( mod(line,200).eq.0 .and. mod(elem,200).eq.0 )then write(6,*) 'image,line,elem,sat_z,sun_z:', > image,line,elem,sat_z,sun_z endif c ----------------------- c 2.1.5. Test for Visible c ----------------------- c During daytime, if albedo>glint [16%] anywhere or if c albedo>bright [8%] in absence of sun glint and coast, c it's cloudy. c c If albedo>> Sunglint or coastal region <<< if (sun_z .eq. twilight) then if (t1(elem,line) .gt. glint) then score = -1 ! invalid or too bright end if c >>> Elsewhere during day <<< else if (sun_z .lt. 90) then if (t1(elem,line) .gt. bright) then score = -1 ! too bright c >>> Discount the score if necessary <<< else if (t1(elem,line) .gt. dark) then score = score*(bright-t1(elem,line))**2*bd_sq end if c >>> Penulty for lack of visible <<< c ... so that all pixels go through the same # of test (discount) else score = score*0.1 end if c >>> Check the score so far <<< if (score .lt. small) then retrv(elem,line) = cloud ! failed ... goto 10 ! ... cycle for the next pixel end if c >>> Note <<< c Part of the test above (albedo>glint) could be done before c navigation to reduce the number of pixels navigated by c 20-30% during day. At night, however, errors in albedo may c incorrectly identify some pixels as cloudy. c c Navigation was done early not only for higher quality at c night but also the fact that the alternative "speeds up" the c processing only during day, which is nice but not enough. c ---------------------------- c 2.1.6. Compute Tb(11)-Tb(12) c ---------------------------- c Wu and Menzel (1998) found that Tb(11)-Tb(12) can be c particularly vulnerable to striping. Assuming spatial c variation of moisture smaller than that of SST, Legeckis c (1997) suggests to smooth Tb(11)-Tb(12) instead of Tb(11) c & Tb(12) individually, thus reducing the striping effect c while preserving the SST features. c c >>> Accumulate <<< c SNM 02/02/01 c >>> do only for G8 and G10 <<< if (hd(3) .le. 76) then t4_t5 = 0 kk = 0 do jj=beg_line,line ! Two lines only for striping do ii=beg_elem,end_elem if (t4(ii,jj).gt.0 .and. t5(ii,jj).gt.0) then t4_t5 = t4_t5 + t4(ii,jj) - t5(ii,jj) kk = kk + 1 end if end do end do c >>> Average <<< c No need to check kk>0 since at least the central pixel c has valid observations. t4_t5 = t4_t5/kk end if c SNM c ----------------------------- c 2.1.7. Compute Tb(3.9)-Tb(11) c ----------------------------- c Unlike Tb(11)-Tb(12), striping is not apparent for Tb(3.9)- c Tb(11). On the other hand, it could be harmful to smooth if c this value is also used for cloud detection. c c Note that validity of Tb(3.9) & Tb(11) for this pixel has c be confirmed at the beginning. c t2_t4 = t2(elem,line) - t4(elem,line) c ----------------------------- c 2.1.8. Atmospheric Correction c ----------------------------- c >>> Day <<< if (sun_z .lt. 90) then if ( iretr_type .eq. 1) then retrv(elem,line) = nint(c1(1,1) + c1(2,1)*t4(elem,line) + + c1(3,1)*t4_t5 + c1(4,1)*t4_t5*t4_t5 + c1(5,1)*sat_z) else if ( iretr_type .eq. 2) then retrv(elem,line) = nint(c2(1,1) + (c2(2,1)*sat_z) + + ((c2(3,1) + (c2(4,1)*sat_z))*t2(elem,line)) + + ((c2(5,1) + (c2(6,1)*sat_z))*t4(elem,line)) + + ((c2(7,1) + (c2(8,1)*sat_z))*t5(elem,line)) ) endif c >>> Night <<< else if ( iretr_type .eq. 1) then retrv(elem,line) = nint(c1(1,2) + c1(2,2)*t2(elem,line) + + c1(3,2)*t2_t4 + c1(4,2)*t2_t4*t2_t4 + c1(5,2)*sat_z) else if ( iretr_type .eq. 2) then retrv(elem,line) = nint(c2(1,2) + (c2(2,2)*sat_z) + + ((c2(3,2) + (c2(4,2)*sat_z))*t2(elem,line)) + + ((c2(5,2) + (c2(6,2)*sat_z))*t4(elem,line)) + + ((c2(7,2) + (c2(8,2)*sat_z))*t5(elem,line)) ) endif end if c ------------------------------ c 2.1.9. Test Against Recent SST c ------------------------------ c There should be little limit on how much the retrieved c SST can be warmer than guess, expecially since if cloud c contaminates SST at one time, the next true SST can appear c much warmer. c c The only worry is that if a hot SST results from errors c in observation and is used to update the guess, it may c disqualify valid SST in future. c c With quality control while reading measurements and smoothing c in updating guess, hopefully this will not be the case. The c upper limit is therefore very generous. c c Since diurnal variation of SST can reach 3K, the lower bound c should be about -3K. c c >>> Find the difference <<< kk = retrv(elem,line) - guess(elem,line) c >>> Invalid if 10K warmer or 3K cooler <<< if (kk.gt.100 .or. kk.lt.-30 .or. & retrv(elem,line).gt.3100) then score = -1 c >>> Discount if more than 1K cooler <<< else if (kk .lt. -10) then score = score*((30+kk)**2)*2.50e-3 end if c >>> Check the score so far <<< if (score .lt. small) then retrv(elem,line) = cloud ! failed ... goto 10 ! ... cycle for the next pixel end if c +++++++++++++++++++++++++ c 2.2. Spectral Consistency c +++++++++++++++++++++++++ c c ----------------------- c 2.2.1. Tb(11) & Tb(12) c ----------------------- c Simulation with 1761 TIGR profiles shows that although c Tb(11)-Tb(12)<2K counts for majority, Tb(11)-Tb(12) can c exceed 3K when very moist atmosphere (total precipitable c water 40-80 g/Kg) is viewed at large angle (60 degree). c c Ideally, one'd be more careful when Tb(11)-Tb(12)>2K. That, c however, would involve checking latitude, season, angle, c and so on, thus impractical. c c Hopefully, erroneously large Tb(11)-Tb(12) leads to large c error in SST that will be detected later. c c On the other hand, Tb(11)-Tb(12)<0 is always undesirable c (calibration error, dust, etc.) c SNM 02/02/01 c >>> do only for G8 and G10 <<< if (hd(3) .le. 76) then if (t4_t5.lt.0 .or. t4_t5.gt.40) then retrv(elem,line) = cloud ! failed ... goto 10 ! ... cycle for the next pixel end if end if c SNM c ------------------------ c 2.2.2. Tb(3.9) & Tb(11) c ------------------------ c At night, Tb(3.9) should be slightly warmer than Tb(11) for c clear sky because tau(3.9)>tau(11). c c Clouds can change this relationship in two ways. c (1) For low stratus, e(3.9) Tb(3.9)SST. c c (2) For partial clouds (semi-transparent or partial IFOV), c Tb(3.9) is much warmer than Tb(11) since dR/dTb(3.9)> c dR/dTb(11) c SNM 02/02/01 c >>> SNM threshold value changed from -3 to -6 (Tb*10) <<< c >>> do only for G8 and G10 (due to t4_t5 term) <<< if (hd(3) .le. 76) then if (t2_t4 .lt. -6) then ! Low stratotus score = -1 else if (sun_z .gt. 90) then ! At night ... xx = abs(t2_t4-2*t4_t5+20) ! predicted diff. if (xx .gt. 10) then score = -1 else if (xx .gt. 5) then score = score*(10-xx)**2*.04 end if end if end if c SNM c >>> Check the score so far <<< if (score .lt. small) then retrv(elem,line) = cloud ! failed ... goto 10 ! ... cycle for the next pixel end if c >>> Note <<< c In the tests above, signal from Tb(12) compared with that c from Tb(11) is stronger in theory but noiser in practice. c c Also, although this is a night test only, no penulty is c imposed because a similar test is done later for day only. c ------------- c 2.2.3. Albedo c ------------- c Albedo over clear ocean should be uniform, but noise and c truncation error can introduce some fluctuation. c c No penulty for night this time because: c 1) A nighttime penulty has been imposed for an albedo test; c 2) Cloud detection at night is more difficult already. c if (sun_z .le. twilight) then ! sun is high c >>> Find the min in the 3-by-3 box <<< kk2 = 999 ! SNM tmp variable to allow comparison with t1 kk = 999 do jj=beg_line,end_line do ii=beg_elem,end_elem if (t1(ii,jj).gt.0) kk=min(t1(ii,jj),kk2) end do end do c >>> Pixel brighter than neighbor? <<< kk = t1(elem,line) - kk if (kk .gt. 10) then ! cloudy if 1% higher score = -1 else if (kk .gt. 2) then ! OK if <0.2% higher score = score*((10-kk)**2)*.016 ! discount otherwise end if c >>> Check the score so far <<< if (score .lt. small) then retrv(elem,line) = cloud ! Too brighter than neighbors goto 10 ! cycle for the next pixel end if end if c --------------- c 2.2.4. Infrared c --------------- c Cloud can introduce spatial variation in Tb's, but so can c certain oceanic fronts. c c However, the impacts of clouds on Tb(3.9) and Tb(11) are c often different, whereas the impacts of oceanic thermal c features under clear sky should be similar for the two Tb's c (not the radiances, though). c c Clear tropical atmosphere, viewed at nadir during day, has c Tb(3.9)-Tb(11)<4K. Sun glint and possibly the turbidity in c the coastal area may increase this difference significantly. c c Experience also shows that, during day, probability of c cloud contamination is high when standard deviation is > 1.5K. c c Note that pixels with Tb(3.9)-Tb(11)<-0.3K should have been c excluded by now. Also, although the following test can be c used at night as well, more rigorous tests should have been c applied for the night cases. c c >>> Penulty for sunglint/coast <<< if (sun_z .eq. twilight) then score = score*.1 c >>> Cloudy if Tb(3.9)-Tb(11) > 8K <<< else if (t2_t4 .gt. 80) then score = -1 c >>> Clear if Tb(3.9)-Tb(11) < 4K <<< c This may not be true at night, but the night cases c have been handled separately. else if (t2_t4 .gt. 40) then c >>> Statistics of Tb(3.9)-Tb(11) <<< call stats(t2,t4,max_elem, & beg_line,end_line,beg_elem,end_elem, mean,vari) c >>> Score <<< if (vari .gt. 400) then score = -1 else score = score*min((80-t2_t4)**2*6.25e-4, & (400-vari)*2.5e-3) end if c >>> No penulty for night <<< c ... a similar test is done for night only end if c >>> Check the score so far <<< if (score .lt. small) then retrv(elem,line) = cloud ! Too brighter than neighbors goto 10 ! cycle for the next pixel end if c ++++++++++++++++++++++++ c 2.3. Temporal Continuity c ++++++++++++++++++++++++ c c ----------------------------------------- c 2.3.1. Compare with Previous Measurements c ----------------------------------------- c If current Tb(11) for a pixel is much cooler than that c shortly before, it is likely due to cloud. c c Sometimes the area is "clearing up", so Tb's for all pixels c are warmer. In this situation, if the Tb increase for the c pixel in question is much lower than neighboring pixels, it c suggests that this pixel has not been as clear as neighbors. c c Sometimes, the area is getting cloudy. In this case, just c because the central pixel is not much cooler than before c is inadequate arguement for clear. c c Finally, if the Tb's of the neighboring pixels do not change c uniformly, it also implies the presence of clouds. c c >>> Statistics of Tb(now)-Tb(earlier) <<< if (has_prev_img) then call stats(t4,t_bef,max_elem, & beg_line,end_line,beg_elem,end_elem, mean,vari) c >>> Temporal change of Tb(11) <<< kk = t4(elem,line) - t_bef(elem,line) kk = min(float(kk),mean,kk-mean) c >>> Cloudy if 1K cooler or stdev > 1.5K <<< if (kk.lt.-10 .or. vari.gt.230) then score = -1 c >>> Discount the score if necessary <<< else if (kk.lt.0 .or. vari.gt.30) then score = score*min((10+kk)**2*.01, (230-vari)*2.5e-5) end if c >>> Penulty for lack of earlier image <<< else score = score*.1 end if c >>> Check the score so far <<< if (score .lt. small) then retrv(elem,line) = cloud ! Too brighter than neighbors goto 10 ! cycle for the next pixel end if c ------------------------------------- c 2.3.2. Compare with Next Measurements c ------------------------------------- c >>> Statistics of Tb(now)-Tb(later) <<< if (has_next_img) then call stats(t4,t_aft,max_elem, & beg_line,end_line,beg_elem,end_elem, mean,vari) c >>> Temporal change of Tb(11) <<< kk = t4(elem,line) - t_aft(elem,line) kk = min(float(kk),mean,kk-mean) c >>> Cloudy if 1K cooler or stdev > 1.5K <<< if (kk.lt.-10 .or. vari.gt.230) then score = -1 c >>> Discount the score if necessary <<< else if (kk.lt.0 .or. vari.gt.30) then score = score*min((10+kk)**2*.01, (230-vari)*2.5e-5) end if c >>> Penulty for lack of later image <<< else score = score*.1 end if c >>> Check the score so far <<< if (score .lt. small) then retrv(elem,line) = cloud ! Too brighter than neighbors goto 10 ! cycle for the next pixel end if c SNM 02/01/24 c ----------------------- c env_check function call c ----------------------- c >> convert angles into radians sun_z_rad = sun_z*a2r sat_z_rad = (sat_z/1000.0)+1.0 sat_z_rad = ACOS(1.0/(sat_z_rad)) c >>> asign BTs in array to variables to allow passing to a function <<< bt_3 = t2(elem,line) bt_11 = t4(elem,line) bt_12 = t5(elem,line) if (hd(3) .le. 76) then ! G8 or G10 chan=0 ! channel combination number (3.9/11) c >>> Call the env_check function <<< call env_check(chan, sat, sat_z, env_dat, bt_3, bt_11, & tol, flag) env_temp1 = flag chan = 1 ! channel combination number (11/12) call env_check(chan, sat, sat_z, env_dat, bt_11, bt_12, & tol, flag) env_temp2 = flag c >>> For G12 use G10, 3.9/11 combination <<< else if (hd(3) .ge. 78) then ! G12 chan=0 call env_check(chan, sat, sat_z, env_dat, bt_3, bt_11, & tol, flag) env_temp1 = flag env_temp2 = flag endif c >>> fail the test if either 3.9/11 or 11/12 tests fail <<< if ((env_temp1 .ne. 1) .OR. (env_temp2 .ne. 1)) then retrv(elem,line) = 3 ! BT fail goto 10 endif c ------------------------- c SNM 02/01/24 >>> end of env_check <<< c SNM 02/01/28 c -------- c Sunglint c -------- c >>> call new sunglint screening for GOES-12 <<< if (hd(3) .ge. 78) then ! G12 and above call sunglint(windspeed, speed, sat_z_rad, sun_z_rad, jday, & hour, minute, orglat, orglon, pxlat, pxlon, glint_rad) c >>> call rad_to_dbt <<< call rad_to_dbt(glint_rad, sat_z_rad, bt_3, & rad_bt_dat, dbt) c >>> call mask_glint <<< dsst_res = maskglint(dbt, sat_z_rad) c >>> if glint, assign a value to retrv <<< if (dsst_res .eq. 0) then retrv(elem,line) = 6 !sunglint goto 10 endif endif c ------------------------------------------- c SNM >>> end of Sunglint <<< 10 continue end do end do c SNM >>> close the data files for env_check and sunglint <<< close(20) close(21) c SNM c ++++++++++++++++++++++ c 2.4. Spatial Coherency c ++++++++++++++++++++++ c It was difficult to detect cloud based on spatial coherency c when spatial variation of SST can be larger than that of c clouds. c c After most clouds have been removed from image, we can c assume that the spatial scale of oceanic features is c larger than residual clouds, therefore we can apply the c spatial coherency principle. c c ----------------------------------------------------- c 2.4.1. Pick a (clear) pixel and define its 9-by-9 box c ----------------------------------------------------- do line=1,hd(9) beg_line = max(line-4,1) end_line = min(line+4,hd(9)) do elem=1,hd(10) if (retrv(elem,line) .gt. cloud) then beg_elem = max(elem-4,1) end_elem = min(elem+4,hd(10)) c ------------------- c 2.4.2. Accumulation c ------------------- kk = 0 mean = 0 vari = 0 do jj=beg_line,end_line do ii=beg_elem,end_elem if (retrv(ii,jj) .gt. cloud) then xx = retrv(ii,jj) mean = mean + xx vari = vari + xx*xx kk = kk + 1 end if end do end do c ----------- c 2.4.3. Test c ----------- if (kk .lt. 3) then ! 2 SST out of 81 - too lonely retrv(elem,line) = cloud else mean = mean/kk vari = sqrt(max(1.e-20, vari/kk-mean*mean)) if (vari.gt.3. .and. mean-retrv(elem,line).gt.vari) then retrv(elem,line) = cloud end if end if end if end do end do c =============== c 3. UPDATE GUESS c =============== c c The basic strategy is to update the guess with a weighted c average of retrieval-guess differences of nearby pixels. c c ++++++++++++++++++++++++++++ c 3.1. Calculate the Increment c ++++++++++++++++++++++++++++ c Since the difference between retrieval and guess will be c used repeatedly, it's better to calculate them once to c save time. It does need some space, but by now many of c those large arrays are no longer used. c do line=1,hd(9) do elem=1,hd(10) if (retrv(elem,line) .gt. cloud) then c >>> SST is retrieved <<< incrm(elem,line) = retrv(elem,line) - guess(elem,line) else c >>> SST is not available <<< incrm(elem,line) = -32767 end if end do end do c +++++++++++++++++++++++ c 3.2. Compute the Weight c +++++++++++++++++++++++ c The weight is reversely proportional to the square c of distance to the pixel in question. c c The central pixel has a larger, but finit, weight. do jj=-4,4 do ii=-4,4 if (ii.eq.0 .and. jj.eq.0) then w(ii,jj) = 2. else w(ii,jj) = 1./(ii*ii+jj*jj) end if end do end do c +++++++++++++++++ c 3.3. Update Guess c +++++++++++++++++ c c -------------------------------------------- c 3.3.1. Pick up a pixel and set up boundaries c -------------------------------------------- do line=1,hd(9) beg_line = max(-4, 1-line) end_line = min(+4,hd(9)-line) do elem=1,hd(10) beg_elem = max(-4, 1-elem) end_elem = min(+4,hd(10)-elem) weighted_sum = 0 total_weight = 0 c -------------------------------------------- c 3.3.2. Compute weighted sum and total weight c -------------------------------------------- c In anticipation of navigation shift following satellite c position change, update is attempted for all pixles, c instead of only for those that is "water" in guess. c c For "space" and "land" pixels, there should be no valid c retrieval in the neighbor. c c Some of the "coastal" pixels will appear as "water" in c the guess. This might be confusing in viewing guess SST c but is actually desirable in processing retrieval. c do jj=beg_line,end_line do ii=beg_elem,end_elem kk = incrm(elem+ii,line+jj) if (kk .gt. -32767) then weighted_sum = weighted_sum + w(ii,jj)*kk total_weight = total_weight + w(ii,jj) end if end do end do c ----------------------------------------------- c 3.3.3. Discount if central SST is not retrieved c ----------------------------------------------- c If SST is only retrieved at a remote neighbor pixel, the c algorithm above tends to give too much weight to that value. c For example, a pixel at the corner normally carries a small c weight of 1/18, but if that the only retrieval in the c 7-by-7 box, it effectively has the total weight. c c The following logic ntends to remedy this case. c c On the other hand, if SST is available from many close c neighbors but not from the central pixel, the following c logic reduce the weight only slightly. c if (incrm(elem,line) .eq. -32767) then total_weight = total_weight + 1 end if c ----------------------- c 3.3.4. Adjust Increment c ----------------------- c Value of xx (0,1] controls how fast the guess is updated. c c xx=0 is "zero update"; the guess will remain unchanged. c xx=1 is "full update"; the guess will (almost) be the retrieval. c xx>1 is "over update"; which probably should never be used. c c A slow update, such as xx=0.8, is often desirable for two c reasons: c (1) Retrieval may have error, so we don't trust it fully. c (2) If guess follows retrieval too closely, it may have c diurnal fluctuation. c c If the retrieval is expected to be different from guess, c for example after revising the algorithm or re-initializing c guess, a faster update may be needed temporarily. c xx = 0.8 c ----------------------- c 3.3.5. Update the guess c ----------------------- c If total weight is one, guess is not updated. This happens c in two situations: c 1) Regardless of the type of pixel, there is no valid SST c retrieval at this time in the neighborhood. c 2) The pixel is not "water". c c If total weight is greater than one, it is possible that the c guess was not "water" and thus has no valid guess value. c In this case we need an estimate. Hopefully, the navigation c shift will be much less than four pixels in any direction c between valid retrievals so that the rough estimate will c be refined a couple of times before it is actually used. c if (total_weight .gt. 1) then c >>> Need update <<< if (guess(elem,line).gt.2700 .and. & guess(elem,line).lt.3100) then c >>> Have a guess for the pixel <<< kk = guess(elem,line) else c >>> Need an estimate of guess <<< c First, there must be at least one pixel in the c neightbor for which guess is available. Otherwise c total weight would be one. c c It's no good to pick the min or max because it may c cause "overflow" at the other limit. And it's not c necessary to use the average. c c It's simpler for program to use the last but faster c for computer to use the first. c do jj=beg_line,end_line do ii=beg_elem,end_elem if (guess(elem+ii,line+jj).gt.2700 .and. & guess(elem+ii,line+jj).lt.3100) then kk = guess(elem+ii,line+jj) goto 30 ! exit end if end do end do 30 continue end if c >>> Update <<< c There is a need to bound the new guess for "coastal" c pixels that can become out of bound [270,310] because c there may never be a retrieval for them. c c Mavor (2000) calculated freezing temperature of -1.9C c (271.3K) for PSU of 35. Strong (2000) confirmed that c SST rarely exceeds 32C (305K). c kk = kk + nint(xx*weighted_sum/total_weight) new_guess(elem) = min(max(2713, kk), 3050) else c >>> No need for update <<< new_guess(elem) = guess(elem,line) end if if (new_guess(elem).lt.0 .or. new_guess(elem).gt.3100 .or. & (new_guess(elem).gt.4 .and. new_guess(elem).lt.2700)) then print*, line, elem, beg_line, beg_elem, end_line, end_elem, & kk, xx, total_weight, weighted_sum print*, 't4' write (*,'(9i8)') ((t4(elem+ii,line+jj), & ii=beg_elem,end_elem),jj=beg_line,end_line) print*, 't5' write (*,'(9i8)') ((t5(elem+ii,line+jj), & ii=beg_elem,end_elem),jj=beg_line,end_line) print*, 'retrv' write (*,'(9i8)') ((retrv(elem+ii,line+jj), & ii=beg_elem,end_elem),jj=beg_line,end_line) print*, 'incrm' write (*,'(9i8)') ((incrm(elem+ii,line+jj), & ii=beg_elem,end_elem),jj=beg_line,end_line) print*, 'guess' write (*,'(9i8)') ((guess(elem+ii,line+jj), & ii=beg_elem,end_elem),jj=beg_line,end_line) print*, 'new_guess' write (*,'(9i8)') (new_guess(elem+ii),ii=beg_elem,end_elem) stop 'update' end if c ------------------------------------- c 3.3.6. Scale the SST to one-byte word c ------------------------------------- if (retrv(elem,line) .lt. 7) then scaled_sst(elem) = retrv(elem,line) else scaled_sst(elem) = min(max(7, & nint((retrv(elem,line)-2700)/1.5)), 255) end if end do c ++++++++++++++++++++++ c 3.4. Write the Results c ++++++++++++++++++++++ c >>> Retrieved SST <<< call wrtara(jmage,line-1,scaled_sst) c >>> Guess SST <<< call wrtara(g_sst,line-1,new_guess) end do call clsara(jmage) call clsara(g_sst) write (*,'("End of goes_sst_navtest3. ",/)') return end c ********************************************** subroutine stats(x,y,n,js,je,is,ie, mean,vari) c ********************************************** c c PURPOSE: c Find the basic statistics of the difference between two 2-D arrays. c c METHOD: c Compare the mean and variance of the difference between the c 3-by-3 arrays observed in two times. c c HISTORY: c 98/10/06 Started. c c I/O LIST: implicit none integer*2 x(n,*), y(n,*) ! IN The two arrays integer n ! IN The 1st dimension of the arrays integer js, je, is, ie ! IN bounds real mean, vari ! OUT mean and variance c ! vari=0 => stats unavailable c c LOCAL VARIABLES: integer i, j, diff, numb c ========== c Initialize c ========== mean = 0. vari = 0. numb = 0 c ========== c Accumulate c ========== do j=js,je do i=is,ie c x & y are quality controlled. Value is negative iff invalid. if (x(i,j).gt.0 .and. y(i,j).gt.0) then diff = x(i,j) - y(i,j) mean = mean + diff vari = vari + diff*diff numb = numb + 1 end if end do end do c ============================= c Compute the mean and variance c ============================= if (numb .gt. 0) then mean = mean/numb vari = amax1(vari/numb-mean*mean, 1.) ! .1% chance vari<1 end if return end c *************************************************************** c nav_pix is now navPixab2 is now separate file c ****************************************************** subroutine read_area(area,line,b_elem,n_elem,band,buf, & val_lo,val_hi,mis_lo,mis_hi) c ****************************************************** c c PURPOSE: c Read a line from an area-file with quality control. c c METHOD: c 1) Read McIDAS area-file 'area' into data buffer 'buf': 'n_elem' c starting with 'b_elem' (zero-based) of 'line' and 'band'. c 2) If mis_lo>mis_hi, missing values are expected: The pixels whose c values are outside of (val_lo,val_hi) are assigned value mis_lo. c 3) Othewise missing value is prohibited: If a pixel's value is outside c of (val_lo,val_hi) and (mis_lo,mis_hi), the program stops. c c NOTE: c Data buffer is of type of 2-byte integer. c c HISTORY: c 99/01/05 Accept a range of valid data and two ways to handle c missing values. c 99/01/31 Expanded to accept a range of missing values. c c I/O LIST: implicit none integer area ! IN area-file number integer line ! IN index of the line to be read integer b_elem ! IN beginning elem to be read integer n_elem ! IN number of elem to be read integer band ! IN index of the band to be read integer*2 buf(n_elem) ! OUT buffer for data integer val_lo, val_hi ! IN range for valid data integer mis_lo, mis_hi ! IN range for missing data c c LOCAL VARIABLES: integer elem c ----------------------------------------- c Read a line using standard McIDAS routine c ----------------------------------------- call redara (area,line,b_elem,n_elem,band,buf) c -------------------------------------------------- c QC Method 1: Unknown values are treated as missing c -------------------------------------------------- if (mis_lo .gt. mis_hi) then do elem=1,n_elem if (buf(elem).lt.val_lo .or. buf(elem).gt.val_hi) & buf(elem) = mis_lo end do c ---------------------------------------- c QC Method 2: Unknown value is prohibited c ---------------------------------------- else do elem=1,n_elem if ((buf(elem).lt.val_lo .or. buf(elem).gt.val_hi) .and. & (buf(elem).lt.mis_lo .or. buf(elem).gt.mis_hi)) then write (*,'(i15," is invalid for AREA",i4.4,6i7)') & buf(elem),area,val_lo,val_hi,mis_lo,mis_hi ,, line,elem stop 'ERROR in reading area-file' end if end do end if return end