subroutine main0 c c PURPOSE: c Co-locate buoy report with GOES measurements. c c Usage: c buoy_cpo_goes.k DATE=ccyyddd HOUR=hh mm GRID=nnnn REGION=xxx c (xxx can be gen, gwn, ges, gws) c Example: c buoy_cpo_goes.k DATE=2009110 HOUR=16 15 GRID=1420 REGION=gen c c INPUT: c Hourly buoy data are available 50 minutes after the hour in c buoy.hh, where hh is hour. c c Hourly GOES Imager data (all bands) and GOES SST are available c two hours after the hour. The SST data are used to indicate c clear/cloudy. c c PROCESSING: c Temporal co-location: time difference among GOES, buoy, and the c nominal hour is less than 30 minutes. c c Spatial co-location: done with standard mcIDAS navigation. c c Ship and duplicate buoy reports: excluded. c c OUTPUT: c Results for each hour is written to a file ccyyddd_hh. These are c formatted, sequential access ASCII files. Following information c was recorded for each co-location: c c INTERNAL VAR EXTERNAL VAR RANGE UNIT c ============================================================= ! write (8,'(9i7,8(11f9.2))') ! & buoy_id(i_buoy),! buoy ID 10000~99999 n/a ! & sat_id, ! satellite ID 70 & 74 n/a ! & year, ! reference CCYY 1970~2069 year ! & month, ! reference month 1~12 month ! & day, ! reference day 1~31 day ! & whole_hour, ! reference hour 0~23 hour ! & buoy_ref, ! buoy-ref time -30~30 min ! & sat_ref, ! sat.-ref time -30~30 min ! & clear, ! # clear pixels 0~9 n/a ! & lat, ! latitude, N. pos. -90~90 degree ! & lon, ! longitude, W. pos. -180~180 degree ! & angles(1), ! sat. zenith angle 0~90 degree ! & angles(2), ! solar zenith angle 0~180 degree ! & angles(3), ! rel. azimuth angle 0~180 degree ! & t_air, ! air temp 0~999 K ! & td, ! dew temp 0~999 K ! & buoy_sst, ! sea temp 0~999 K ! & wind_dir, ! wind direction 0~360 degree ! & wind_spd, ! wind speed 0~99 m/s ! & mslp, ! mean sea level pres. 0~9999 mb c These are followed by GOES data as 11-by-7 arrays. The 11 elements c correspond to the 9 pixels (#5 is the buoy location) and the average c and standard deviation of those clear pixels. This repeats 7 times c for albedo in %, Tb(2-5) in K, instantaneous SST in K, and the c archived SST in K. Note that the "average" admits only those for c which valid instantaneous SST is available. ! & ((buf(i,j),i=1,11),j=1,7) c ============================================================= c c NOTE: c Current implementation is less efficient in that the area files c may have been open-set-close many times during the execution. c c This should be insignificant since the computer is so powerful c now and there are only a few hundreds of buoy reports anyway. c c If this becomes a concern, however, one can first sort the input c data according to satellite domain. c c HISTORY: c 99/10/23 Started. c 99/10/30 Get all buoy reports even when GOES or SST are not available. c 99/11/17 Correct one error when GOES data is unavailable. c 00/06/02 Chang offset from 271 to 270 for new SST. c 00/06/15 Change buoy file from buoy.hh to buoysst.hh c 00/07/07 Fixed a bug: writing hour to buoy_file(6:7) c 00/07/21 Process one (instead of two) hour earlier buoy data c 00/08/01 In "dual processing", images close to the hour are 131-431. c 03/10/31 Get info from Pclear area file; expand buf to hold c 04/12/09 c Code extended to write both short ascii file and extended binary file c 05/01/17 c Code modifed to read variables in from new multi-band AREA files c 05/04/14 c Fixed a bug: buoy_SST value is added to output binary MDB file c 05/11/11 Bayesian outputs now correspond to Buoy location c NWP data now correspond to Buoy location c Changed NWP surface pressure to NWP 1000 mbar altitude c ADD_slot passed to no_goes function c ADD_slot(2) = image line number c ADD_slot(3) = number of lines in image c Interpolated fileds now correctly interpolated c 05/11/30 Corrected scaling of SST error c Corrected scaling of FFM outputs c Added windspeed to interpolated outputs c ADD_slot(4) = actual scan start c 09/03/11 1) Changed output file name c 2) pass GRID file number from argument to function no_goes() c 3) process buoy match seperately for each quadrant cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit none c///////////////////////////////////////////////////// cccccccccc changes made on 03-11-2009 ccccccccccccc c -- for output files c character*12 buoy_file, match_file, sst_file character*80 buoy_file, match_file, sst_file character*80 outdir, outname character*10 whoami character*11 cepoch integer len1, len2, len_trim integer iccyy integer iddd integer ihh c -- for grid files integer grid1 integer grid2 c -- for quadrant integer quad_number character*3 region_name, ckwp cccccccccc End of changes ccccccccccccccc c//////////////////////////////////////////////////// cc character*12 match_file_binary cc real buf(11,7), angles(3) cc real buf(11,8), angles(3) real buf(11,9), angles(3) logical no_goes logical first_call integer mcgetdaytime, mcdmytocyd integer mcsectodaytime, mcdaytimetosec integer ref_sec, sat_ref, buoy_ref, sat_id integer n_rpt, inside, n_ship, repeat, i_buoy integer i, j, buoy_id(3000), clear, whole_hour integer II c Variables to handle the four quadrants integer quad integer goes_img(4) integer min_lon(4) integer max_lon(4) integer min_lat(4) integer max_lat(4) c Variables to hold buoy report integer year, month, day, hour, minute, ikwp real lon,lat,buoy_sst,x,mslp,t_air,td,wind_dir,wind_spd c Variables to test for ship ID character*8 id character*1 zero, nine, cbuf(8) equivalence (cbuf(1),id) data buoy_file/'buoysst.hh'/,match_file/'ccyy_ddd.hhS'/, & sst_file/'ccyy_ddd_hh'/ data goes_img/231, 431, 131, 331/ data min_lon/ 90, 105, 30, 30/ data max_lon/180, 180, 120, 120/ data min_lat/ 0, -45, -20, -45/ data max_lat/ 45, 0, 45, -20/ data zero/'0'/ data nine/'9'/ data first_call / .true. / c GOES SST variables real ref_sat_1(49) ! Satellite % VIS real ref_sat_2(49) ! Satellite % Red reflectance real ref_sat_3(49) ! Satellite % VNIR reflectance real bt1_sat(49) ! 3.9um or equivalent (K) real bt2_sat(49) ! 6.7um WV channel (K) [GOES=CH3] real bt3_sat(49) ! 8.7um (K) real bt4_sat(49) ! 11um (K) real bt5_sat(49) ! 12um (K) real bt6_sat(49) ! 13um (K) real sst_prd(49) ! Satellite SST product (K) real sst_prd_err(49) ! Satellite SST product random error (K) real sst_prd_bias(49) ! Satellite SST product bias if cloudy (K) real pclear(49) ! Clear sky probability - first pass real pclear_spec(49) ! Clear sky probability - spectral real pclear_text(49) ! Clear sky probability - textural real pclear_revs(49) ! Clear sky probability - revised (clumpy) real prior_sst(49) ! Prior SST analysis real ref_ffm(3) ! FFM reflectance at pixel (VIS) real bt_ffm(6) ! FFM BT at pixel 3.9um real dBTdSST_ffm(6) ! FFM dBT1/dSST at pixel real dBTdTCWV_ffm(6) ! FFM dBT1/dTCWV at pixel real dBTdTair_ffm(6) ! FFM dBT1/dTATM at pixel real dBT_glint ! dBT1 due to sunglint (K) real dBT_scatt ! dBT1 due to scattering (K) real NWP_month ! NWP analaysis month real NWP_day ! NWP analysis day real NWP_hour ! NWP analysis hour real NWP_SST(4) ! NWP prior skin SSTs around pixel real NWP_SLP(4) ! NWP sea level pressure around pixel real NWP_10mV(4) ! NWP 10m wind north around pixel real NWP_10mU(4) ! NWP 10m wind east around pixel real NWP_2mRH(4) ! NWP 2m relative humidity around pixel real NWP_2mT(4) ! NWP 2m temperature around pixel real NWP_RH_prof(4,18) ! NWP realtive humidity profiles around pixel real NWP_TCWV(4) ! NWP total column water vapour around pixel real NWP_T_prof(4,18) ! NWP temperature profiles around pixel real NWP_arsl(4,6) ! NWP aerosol parameters around pixel real Diurnal_SST ! Diurnal skin SST model real NWP_WS_intrp(5) ! NWP interpolated wind speed real NWP_SSI_intrp(5) ! NWP interpolated SSI real NWP_DLI_intrp(5) ! NWP interpolated DLI real NWP_LHF_intrp(5) ! NWP interpolated LHF real NWP_SHF_intrp(5) ! NWP interpolated SHF real NWP_ST_intrp(5) ! NWP interpolated skin ST real ASST ! SST grided value at buoy location c C MDB binary write functions integer openf integer wrtrec integer closef integer c_out integer fd character*32 cCode, reg_coeff character*16 ncep_ver,ffm_ver,ref_ver logical file_open integer start_number integer proc_number integer goes_number_in c c********additional slot for 4 values real ADD_slot(4) c////////////////////////////////////////////////////////// ccccccccccccccc Added on 03-11-2009 ccccccccccccccccccc c -- to get grid files from argument grid1 = ikwp('GRID', 1, -1) if (grid1.lt.0) then write (*,'("Wrong grid", i4)') grid1 return endif grid2 = grid1 + 1 C grid2 = ikwp('GRID', 2, grid2) print*, 'grid1=',grid1,' grid2=',grid2 c -- to get quadrant info region_name = ckwp('REGION', 1, ' ') print*,"region_name = ",region_name if (region_name(1:3) .eq. 'gwn') then quad_number = 1 start_number = 200 else if (region_name(1:3) .eq. 'gws') then quad_number = 2 start_number = 400 else if (region_name(1:3) .eq. 'gen') then quad_number = 3 start_number = 100 else if (region_name(1:3) .eq. 'ges') then quad_number = 4 start_number = 300 else write (*,'("Wrong REGION NAME")') return endif print*,'quad_number = ', quad_number II = quad_number print*,'DBG- start_number = ', start_number print*,'DBG- goes_img(II) = ', goes_img(II) print*,'DBG- longitude range',min_lon(II),max_lon(II) print*,'DBG- latitude range',min_lat(II),max_lat(II) proc_number = ikwp('PROC_NUMBER',1,-1) if (proc_number.lt.0) then write (*,'("Wrong proc_number", i4)') proc_number return endif goes_number_in = start_number + proc_number ccccccccccccccc End of adding cccccccccccccccccccc c////////////////////////////////////////////////////////// c ===================================== c Determine and open input/output files c ========== =========================== c >>> Get current ccyyddd & hh0000 <<< day = ikwp('DATE', 1, -1) whole_hour = ikwp('HOUR', 1, -1)*10000 if (day.lt.0 .or. whole_hour.lt.0) then write (*,'("Wrong date or hour",2i9)') day, whole_hour return end if c if (mcgetdaytime(day,whole_hour) .ne. 0) then c write (*,'("unable to get current day/time - will stop")') c stop 'no current day/time' c end if c >>> Convert the time to second since 1970 <<< if (mcdaytimetosec(day,whole_hour, ref_sec) .ne. 0) then write (*,'("sec failed - stop",3i9)') day,whole_hour,ref_sec stop 'no second since 1970' end if print*, 'current day/hour/sec', day, whole_hour, ref_sec c >>> Set reference time <<< c Step 1: Convert to hours since 1970 c Step 2: Round to the nearest whole hour c Step 3: Subtract one hours c Step 4: Convert this reference time to second since 1970 c ref_sec = (nint(ref_sec/3600.)-1)*3600 c >>> Define input/output files by ref. time <<< c Step 1: Convert the reference time to ccyyddd & hhmmss if (mcsectodaytime(ref_sec, day,whole_hour) .ne. 0) then write (*,'("sec failed - stop",3i9)') ref_sec,day,whole_hour stop 'no day/time from second since 1970' end if print*, 'updated day/hour/sec', day, whole_hour, ref_sec c Step 2: Use ccyyddd as the first part of output file name write (match_file(1:4),'(i4.4)') day/1000 write (match_file(6:8),'(i3.3)') mod(day,1000) write ( sst_file(1:4),'(i4.4)') day/1000 write ( sst_file(6:8),'(i3.3)') mod(day,1000) c Step 3: Convert hhmmss to whole hour whole_hour = nint(whole_hour/10000.) c Step 4: Use whole hour as the second part of I/O file names write (match_file(10:11),'(i2.2)') whole_hour ! ccyyddd_hh write ( sst_file(10:11),'(i2.2)') whole_hour write ( buoy_file( 9:10),'(i2.2)') whole_hour ! buoy.hh c/////////////////////////////////////////////////////////////////// ccccccccccccccccc Added on 03-11-2009 cccccccccccccccccccc c -- to get directory call getenv('LOGNAME',whoami) len1 = len_trim(whoami) write(6,*) 'whoami=',whoami(1:len1) c outdir = '/disk2/pub/' // whoami(1:len1) // '/buoyMatch/' outdir = '/disk2/pub/SST/buoyMatch/' len1 = len_trim(outdir) c -- to get file name c cepoch iccyy = day/1000 iddd = mod(day,1000) c ihh = ikwp('HOUR', 1, -1) ihh = whole_hour write (cepoch,'(i4.4,"_",i3.3,".",i2.2)') iccyy, iddd,ihh print*,'cepoch=',cepoch c filename outname = 'match_' // region_name(1:3) // '_' // cepoch // 'S' len2 = len_trim(outname) c -- whole name for output c match_file = outdir(1:len1) // outname(1:len2) match_file = outname(1:len2) c -- buoydata name WRITE(buoy_file,'(a,a,''.'',i2.2)')'buoydata.', & region_name(1:3),proc_number ccccccccccccccccc End of adding cccccccccccccccccccc c/////////////////////////////////////////////////////////////////// print*, 'match file name: ', match_file print*, 'sst file name: ', sst_file print*, 'buoy file name: ', buoy_file c >>> Open input/output files <<< open (7, file=buoy_file, status='old', form='formatted') open (8, file=match_file, form='formatted', access='sequential') ccc open (9, file=sst_file, form='unformatted', status='old', ccc & access='direct', recl=3000/4) cCode = 'buoy_cpo.k ' ncep_ver = 'NCEP v999 ' ffm_ver = 'OPTRAN v1.3 ' ref_ver = 'STUART v1.0 ' reg_coeff = 'reg_coeff_mer.asc ' write(6,*) 'DBG-Opening binary matchup file' c Swap args fd,match_file - bob/gordana 2012jun05 cccc c_out = openf(fd,match_file,cCode,ncep_ver,ffm_ver,ref_ver, c_out = openf(match_file,fd,cCode,ncep_ver,ffm_ver,ref_ver, & reg_coeff) if (c_out .ne. 0) then file_open = .false. write(6,*) 'DBG-file_open is false' else file_open = .true. write(6,*) 'DBG-file_open is true' end if write(6,*) 'DBG-match_file ', match_file write(6,*) 'DBG-cCode ', cCode write(6,*) 'DBG-ncep_ver ', ncep_ver write(6,*) 'DBG-ffm_ver ', ffm_ver write(6,*) 'DBG-ref_ver ', ref_ver write(6,*) 'DBG-reg_coeff ', reg_coeff write(6,*) 'DBG-c_out ', c_out write(6,*) 'DBG-file_open ', file_open c ============================ c Loop over the four quadrants c ============================ c do quad=1,4 do quad=quad_number, quad_number ! change on 03-11-2009 c ---------------- c Reset everything c ---------------- rewind (7) inside = 0 repeat = 0 n_ship = 0 i_buoy = 0 do 200 n_rpt=0,9999 c ----------------------------- c Read a report, skip wave info c ----------------------------- read (7,'(a8,i5,4i3,2f9.2,13f8.2)',end=201) & id, year, month, day, hour, minute, & lon, lat, ! West pos. and North pos., in degree & buoy_sst, x, ! SST in K & C & mslp, ! mean sea level pressure in mb & t_air, td, ! air temperature and dew point in K & wind_dir, wind_spd ! wind direction in degree, speed in m/s c -------------------------------------- c Report is inside the satellite domain? c -------------------------------------- if (lat.lt.min_lat(quad) .or. lat.gt.max_lat(quad) .or. & lon.lt.min_lon(quad) .or. lon.gt.max_lon(quad)) goto 200 inside = inside + 1 c ------------------------------------------ c Report is within the required time window? c ------------------------------------------ c >>> Step 1: Convert yymmdd to ccyyddd <<< if (mcdmytocyd(day,month,year, i) .ne. 0) then write (*,'("mcdmytocyd failed for",3i9)') day,month,year goto 200 end if c >>> Step 2: Convert to second since 1970 <<< j = hour*10000 + minute*100 if (mcdaytimetosec(i,j, buoy_ref).ne.0) then write (*,'("mcdaytimetosec 3 failed",8i9)') i, j goto 200 end if c >>> Step 3: Compare buoy & reference time <<< buoy_ref = nint((buoy_ref-ref_sec)/60.) if (abs(buoy_ref) .gt. 30) then write (*,'("large buoy-ref",2i12)') ref_sec, buoy_ref goto 200 end if c ----------------------- c Report is made by ship? c ----------------------- do i=1,5 if (cbuf(i).lt.zero .or. cbuf(i).gt.nine) then ! ship: n_ship = n_ship + 1 ! ... increment ship report ... goto 200 ! ... and try the next report end if end do c ---------------------- c Report is a duplicate? c ---------------------- i_buoy = i_buoy + 1 ! # unique buoys read (id,'(i8)') buoy_id(i_buoy) ! save as the last new buoy ID do i=i_buoy-1,1,-1 if (buoy_id(i_buoy) .eq. buoy_id(i)) then ! duplicate: repeat = repeat + 1 ! ... increment repeat ... i_buoy = i_buoy - 1 ! ... decrement buoy list ... goto 200 ! ... and try the next buoy end if end do c --------------- c Zero MDB Record c --------------- dBT_glint=0.0 dBT_scatt=0.0 NWP_month=0.0 NWP_day=0.0 NWP_hour=0.0 Diurnal_SST=0.0 do i=1,49 ref_sat_1(i)=0.0 ref_sat_2(i)=0.0 ref_sat_3(i)=0.0 bt1_sat(i)=0.0 bt2_sat(i)=0.0 bt3_sat(i)=0.0 bt4_sat(i)=0.0 bt5_sat(i)=0.0 bt6_sat(i)=0.0 sst_prd(i)=0.0 sst_prd_err(i)=0.0 sst_prd_bias(i)=0.0 pclear(i)=0.0 pclear_spec(i)=0.0 pclear_text(i)=0.0 pclear_revs(i)=0.0 prior_sst(i)=0.0 if (i .lt. 7) then bt_ffm(i)=0.0 dBTdSST_ffm(i)=0.0 dBTdTCWV_ffm(i)=0.0 dBTdTair_ffm(i)=0.0 end if if (i .lt. 6) then NWP_WS_intrp(i)=0.0 NWP_SSI_intrp(i)=0.0 NWP_DLI_intrp(i)=0.0 NWP_LHF_intrp(i)=0.0 NWP_SHF_intrp(i)=0.0 NWP_ST_intrp(i)=0.0 end if if (i .lt. 5) then NWP_SST(i)=0.0 NWP_SLP(i)=0.0 NWP_10mV(i)=0.0 NWP_10mU(i)=0.0 NWP_2mRH(i)=0.0 NWP_2mT(i)=0.0 NWP_TCWV(i)=0.0 do j=1,18 NWP_RH_prof(i,j)=0.0 NWP_T_prof(i,j)=0.0 if (j .lt. 7) then NWP_arsl(i,j)=0.0 end if end do end if if (i .lt. 4) then ref_ffm(i)=0.0 end if end do c ----------------------------------------- c Get GOES data or fill with missing values c ----------------------------------------- Cc if (no_goes(goes_img(quad),ref_sec,lat,lon, Cc & sat_id,sat_ref,clear,angles,buf)) then cccc write(6,*)'Calling no_goes' c if (no_goes(first_call,goes_img(quad), ref_sec, lat, lon, if (no_goes(first_call,goes_number_in, ref_sec, lat, lon, & sat_id, sat_ref, clear, angles, buf, & ref_sat_1, ref_sat_2, ref_sat_3, & bt1_sat, bt2_sat, bt3_sat, bt4_sat, bt5_sat, bt6_sat, & sst_prd, sst_prd_err, sst_prd_bias, & pclear, pclear_spec, pclear_text, pclear_revs, prior_sst, & ref_ffm, bt_ffm, & dBTdSST_ffm, dBTdTCWV_ffm, dBTdTair_ffm, & dBT_glint, dBT_scatt, & NWP_month, NWP_day, NWP_hour, & NWP_SST, NWP_SLP, NWP_10mV, NWP_10mU, NWP_2mRH, NWP_2mT, & NWP_RH_prof, NWP_TCWV, NWP_T_prof, NWP_arsl, Diurnal_SST, & NWP_WS_intrp, NWP_SSI_intrp, NWP_DLI_intrp, c & NWP_LHF_intrp, NWP_SHF_intrp, NWP_ST_intrp, ADD_slot)) then c changes on 03-11-2009 & NWP_LHF_intrp, NWP_SHF_intrp, NWP_ST_intrp, ADD_slot, & grid1,grid2)) then sat_id = -6 sat_ref = -999 clear = -6 do i=1,3 angles(i) = -6 end do do j=1,9 do i=1,11 buf(i,j) = -6 end do end do c 99/11/17 Wrong: When there is no GOES data, this removes the correct c buoy id and uses the previous buoy id. c i_buoy = i_buoy - 1 ! ... decrement buoy list ... c ----------------------------------------- c Final check: buoy time and satellite time c ----------------------------------------- else if (abs(sat_ref-buoy_ref) .gt. 30) then write (*,'("large sat-buoy",2i9)') sat_ref,buoy_ref i_buoy = i_buoy - 1 ! ... decrement buoy list ... goto 200 end if c************pick up the gridded SST value c 11 Nov 2005 - O. Embury - ADD_slot now passed directly to no_goes c ADD_slot(1) = ASST c************************************** c -------------------- c Save the co-location c -------------------- c INTERNAL VAR EXTERNAL VAR RANGE UNIT c ============================================================= cc write (8,'(9i7,8(11f9.2))') cc write (8,'(9i7,9(11f9.2))') write (8,'(9i7,10(11f9.2))') & buoy_id(i_buoy),! buoy ID 10000~99999 n/a & sat_id, ! satellite ID 70 & 74 n/a & year, ! reference CCYY 1970~2069 year & month, ! reference month 1~12 month & day, ! reference day 1~31 day & whole_hour, ! reference hour 0~23 hour & buoy_ref, ! buoy-ref time -30~30 min & sat_ref, ! sat.-ref time -30~30 min & clear, ! # clear pixels 0~9 n/a & lat, ! latitude, N. pos. -90~90 degree & lon, ! longitude, W. pos. -180~180 degree & angles(1), ! sat. zenith angle 0~90 degree & angles(2), ! solar zenith angle 0~180 degree & angles(3), ! rel. azimuth angle 0~180 degree & t_air, ! air temp 0~999 K & td, ! dew temp 0~999 K & buoy_sst, ! sea temp 0~999 K & wind_dir, ! wind direction 0~360 degree & wind_spd, ! wind speed 0~99 m/s & mslp, ! mean sea level pres. 0~9999 mb c These are followed by GOES data as 11-by-6 arrays. The 11 elements c correspond to the 9 pixels (#5 is the buoy location) and the average c and standard deviation of those clear pixels. This repeats 6 times c for albedo in % and Tb(2-5) and SST in K. & ((buf(i,j),i=1,11),j=1,9) if (file_open) then cccc write(6,*)'Writing record to binary MDB' c_out = wrtrec(fd, buoy_id(i_buoy), sat_id,year, & month, day, whole_hour, buoy_ref, sat_ref, lat, lon, & angles(1), angles(2), angles(3), & t_air, td, buoy_SST, wind_dir, wind_spd, mslp, & ref_sat_1, ref_sat_2, ref_sat_3, & bt1_sat, bt2_sat, bt3_sat, bt4_sat, bt5_sat, bt6_sat, & sst_prd, sst_prd_err, sst_prd_bias, & pclear, pclear_spec, pclear_text, pclear_revs, prior_sst, & ref_ffm, bt_ffm, & dBTdSST_ffm, dBTdTCWV_ffm, dBTdTair_ffm, & dBT_glint, dBT_scatt, & NWP_month, NWP_day, NWP_hour, & NWP_SST, NWP_SLP, NWP_10mV, NWP_10mU, NWP_2mRH, NWP_2mT, & NWP_RH_prof, NWP_TCWV, NWP_T_prof, NWP_arsl, Diurnal_SST, & NWP_WS_intrp, NWP_SSI_intrp, NWP_DLI_intrp, & NWP_LHF_intrp, NWP_SHF_intrp, NWP_ST_intrp, ADD_slot) if (c_out .ne. 0) then write(6,*)'Unable to write to ',fd end if end if c ============================================================= 200 continue c ---------------------------------- c Report summary about this quadrant c ---------------------------------- 201 continue write (*,'(a12,i2,i4," reports",i4," inside",i4," ship", & i4," duplicate",i4," unique buoy")') & match_file, quad, n_rpt, inside, n_ship, repeat, i_buoy end do close (7) close (8) close (9) if (file_open) c_out = closef(fd) return end c **************************************** C function no_goes(image,ref_sec,lat,lon, C & sat_id,sat_ref,clear,angles,buf) c function no_goes(first_call,image, ref_sec, lat, lon, & sat_id, sat_ref, clear, angles, buf, & ref_sat_1, ref_sat_2, ref_sat_3, & bt1_sat, bt2_sat, bt3_sat, bt4_sat, bt5_sat, bt6_sat, & sst_prd, sst_prd_err, sst_prd_bias, & pclear, pclear_spec, pclear_text, pclear_revs, prior_sst, & ref_ffm, bt_ffm, & dBTdSST_ffm, dBTdTCWV_ffm, dBTdTair_ffm, & dBT_glint, dBT_scatt, & NWP_month, NWP_day, NWP_hour, & NWP_SST, NWP_SLP, NWP_10mV, NWP_10mU, NWP_2mRH, NWP_2mT, & NWP_RH_prof, NWP_TCWV, NWP_T_prof, NWP_arsl, Diurnal_SST, & NWP_WS_intrp, NWP_SSI_intrp, NWP_DLI_intrp, c & NWP_LHF_intrp, NWP_SHF_intrp, NWP_ST_intrp, ADD_slot) c changed on 03-11-2009 & NWP_LHF_intrp, NWP_SHF_intrp, NWP_ST_intrp, ADD_slot, & grid1,grid2) c **************************************** c c PURPOSE: c Get GOES measurements and the derived SST. c c HISTORY: c 99/10/28 Separated from main. c 99/11/17 old_image was updated once old_image is closed and before all c the checks. Theerfore if the image is wrong (e.g., non-exist), it'll c be caught in the first round but bypasses all checks in subsequent c rounds. For that reason, old_image now is updated at the end of all c checks. c 04/12/09 CPO c Code modified to generate extended MDB for full Bayes screening output c Buffer has been replaced with a large variable list c All bayes fields added in read loop c Code added to read GRID files for extracting NWP fieldsbuoy_cpo_LATEST.for c Unavailable fields currently zeroed c I/O LIST: implicit none INCLUDE 'gridparm.inc' ! GRID file parameters c I/O LIST: integer grid1, grid2 ! add on 03-10-2009 -- NCEP profile & surface data integer image ! IN file for GOES image integer ref_sec ! IN reference time in second real lat, lon ! IN buoy location integer sat_id ! OUT satellite ID integer sat_ref ! OUT sat. time - ref. time integer clear ! OUT # clear pixels real angles(3) ! OUT viewing geometry real buf(11,9) ! OUT GOES data & SST short list real ref_sat_1(49) ! Satellite % VIS real ref_sat_2(49) ! Satellite % Red reflectance real ref_sat_3(49) ! Satellite % VNIR reflectance real bt1_sat(49) ! 3.9um or equivalent (K) real bt2_sat(49) ! 6.7um WV channel (K) [GOES=CH3] real bt3_sat(49) ! 8.7um (K) real bt4_sat(49) ! 11um (K) real bt5_sat(49) ! 12um (K) real bt6_sat(49) ! 13um (K) real sst_prd(49) ! Satellite SST product (K) real sst_prd_err(49) ! Satellite SST product random error (K) real sst_prd_bias(49) ! Satellite SST product bias if cloudy (K) real pclear(49) ! Clear sky probability - first pass real pclear_spec(49) ! Clear sky probability - spectral real pclear_text(49) ! Clear sky probability - textural real pclear_revs(49) ! Clear sky probability - revised (clumpy) real prior_sst(49) ! Prior SST analysis real ref_ffm(3) ! FFM reflectance at pixel (VIS) real bt_ffm(6) ! FFM BT at pixel 3.9um real dBTdSST_ffm(6) ! FFM dBT1/dSST at pixel real dBTdTCWV_ffm(6) ! FFM dBT1/dTCWV at pixel real dBTdTair_ffm(6) ! FFM dBT1/dTATM at pixel real dBT_glint ! dBT1 due to sunglint (K) real dBT_scatt ! dBT1 due to scattering (K) real NWP_month ! NWP analaysis month real NWP_day ! NWP analysis day real NWP_hour ! NWP analysis hour real NWP_SST(4) ! NWP prior skin SSTs around pixel real NWP_SLP(4) ! NWP sea level pressure around pixel real NWP_10mV(4) ! NWP 10m wind north around pixel real NWP_10mU(4) ! NWP 10m wind east around pixel real NWP_2mRH(4) ! NWP 2m relative humidity around pixel real NWP_2mT(4) ! NWP 2m temperature around pixel real NWP_RH_prof(4,18) ! NWP realtive humidity profiles around pixel real NWP_TCWV(4) ! NWP total column water vapour around pixel real NWP_T_prof(4,18) ! NWP temperature profiles around pixel real NWP_arsl(4,6) ! NWP aerosol parameters around pixel real Diurnal_SST ! Diurnal skin SST model real NWP_WS_intrp(5) ! NWP interpolated wind speed real NWP_SSI_intrp(5) ! NWP interpolated SSI real NWP_DLI_intrp(5) ! NWP interpolated DLI real NWP_LHF_intrp(5) ! NWP interpolated LHF real NWP_SHF_intrp(5) ! NWP interpolated SHF real NWP_ST_intrp(5) ! NWP interpolated skin ST real NWP_monthc ! NWP analaysis month real NWP_dayc ! NWP analysis day real NWP_hourc ! NWP analysis hour real NWP_SSTc(4) ! NWP prior skin SSTs around pixel real NWP_SLPc(4) ! NWP sea level pressure around pixel real NWP_10mVc(4) ! NWP 10m wind north around pixel real NWP_10mUc(4) ! NWP 10m wind east around pixel real NWP_2mRHc(4) ! NWP 2m relative humidity around pixel real NWP_2mTc(4) ! NWP 2m temperature around pixel real NWP_RH_profc(4,18) ! NWP realtive humidity profiles around pixel real NWP_TCWVc(4) ! NWP total column water vapour around pixel real NWP_T_profc(4,18) ! NWP temperature profiles around pixel real NWP_arslc(4,6) ! NWP aerosol parameters around pixel real Diurnal_SSTc ! Diurnal skin SST model real NWP_WS_intrpc(5) ! NWP interpolated wind speed real NWP_SSI_intrpc(5) ! NWP interpolated SSI real NWP_DLI_intrpc(5) ! NWP interpolated DLI real NWP_LHF_intrpc(5) ! NWP interpolated LHF real NWP_SHF_intrpc(5) ! NWP interpolated SHF real NWP_ST_intrpc(5) ! NWP interpolated skin ST c real ADD_slot(4) logical no_goes ! OUT true if no GOES data logical first_call c c LOCAL VARIABLES: integer jmage, nvset, lit, hd(64), hd1(64), old_image c integer kmage, hd2(64) ! for pclear integer hd2(64) ! for pclear c integer kmage2, hd3(64) ! for sigma integer hd3(64) ! for sigma integer allimage(9) integer nimages, ISTAT integer opnara, nv1eas, nv1opt, mcdaytimetosec, sav_tim integer*2 tb, tb3(3), tb7(7) integer*1 sst_buf(3000) real line, elem, xin(4) real ASST ! SST grided value at buoy location integer i, j, k, kk,index data old_image/-1/ c NWP timestampe variables integer day, month, year c NWP variables in GRID1000 & GRID1001 integer lev_t(181,360,18) ! Profiles of temperature on levels (K) integer lev_rh(181,360,18) ! Profiles of relative humidity on levels (%) integer z1(181,360) ! Altitude (m) at 1000.0 mB integer surf_t(181,360) ! Surface temperature (K) integer surf_rh(181,360) ! Surface relative humidity (%) integer surf_u(181,360) ! Surface windspeed, u (m / s) integer surf_v(181,360) ! Surface windspeed, v (m / s) integer tcwv(181,360) ! TCWV integer t_at2m(181,360) ! Temperature at 2m level (K) integer surf_SHF(181,360) ! SHORT H.F. integer surf_LHF(181,360) ! LONG H.F. integer surf_P(181,360) ! Surface pressure integer rh_at2m(181,360) ! relative humidity (%) at 2m c NWP grid indices integer nlat, nlon integer latlo, lathi, lonlo, lonhi real latf, lonf integer hdgrid(64) ! for GRID files integer status ! Status flag for IGGET integer nrows, ncols ! Number of rows and columns in GRID file integer grid(181,360) ! Temporary array for storing grids integer ii ! Counter integer mccydtodmy ! McIDAS routine in gettime.c integer mciydtocyd integer IGGET integer icyd integer c_out integer ifloor save NWP_monthc ! NWP analaysis month save NWP_dayc ! NWP analysis day save NWP_hourc ! NWP analysis hour save lev_t ! Profiles of temperature on levels (K) save lev_rh ! Profiles of relative humidity on levels (%) save z1 ! Altitude (m) at 1000.0 mB save surf_t ! Surface temperature (K) save surf_rh ! Surface relative humidity (%) save surf_u ! Surface windspeed, u (m / s) save surf_v ! Surface windspeed, v (m / s) save tcwv ! TCWV save t_at2m ! Temperature at 2m level (K) save surf_SHF ! SHF save surf_LHF ! LHF save surf_P ! Surface pressure save rh_at2m ! relative humidity (%) at 2m c c ----------- c Set default c ----------- no_goes = .true. ! no GOES data if returned half way jmage = image + 11 ! ... corresponding SST c kmage = jmage + 1000 ! ... corresponding p_clear c kmage2 = jmage + 2000 ! ....corresponding sigma C CPO 13/01/2005 Include number of images variable nimages nimages=3 allimage(1) = jmage do i=2,nimages allimage(i)=jmage+1000*(i-1) end do C CPO 13/01/2005 Following has been superceded CXXc 0000 Retrieved SST CXXc 1000 First-pass clear-sky probability CXXc 2000 Revised clear-sky probability CXXc 3000 Satellite product SST random error CXXc 4000 First-pass clear-sky p spectral CXXc 5000 First-pass clear-sky p texture CXXc 6000 FFM Reflectance CXXc 7000 dBT1 from scattering estimate CXXc 8000 dBT1 from sunlight estimate C-------------------------------------------------------------------------- C C CPO 13 JAN 2005 C C SATBAND file modified to include the following definitions of the C Bayesian Cloud Screening Products for GOES satellites C C SAT = 70 72 74 76 78 (i.e. imager only) C CAL = 'PRD 'buoy_cpo_LATEST.for C C Product Band Number and Description C C 01 GOES SST, UW-CIMSS, 00/05/01 (unmasked) C 02 SST error (sigma) C 03 Clear sky probability (Bayesian) C 04 Adjusted clear sky prob. (Bayesian) C 05 Clear sky probability (spectral only) C 06 Clear sky probability (textural only) C 07 Cloud mask (Bayesian prob.) C 08 Cloud mask (adjusted Bayesian prob.) C 09 GOES SST, UW-CIMSS, 00/05/01 (masked) C 10 Modelled reflectance C 11 3.9um scattering correction (daytime) C 12 3.9um sun glint correction (daytime) C 13 dBT(3.9)/dSST C 14 dBT(3.9)/dTCWV C 15 dBT(11)/dSST C 16 dBT(11)/dTCWV C C---------------------------------------------------------------------------- C CPO 13/01/2005 Modified output from GOES Cloud Screening C Products now stored in multi-band AREA files C C 0000 Retrieved SST + error (bands 1,2) C 1000 Bayesian clear sky probability fields (bands 3,4,5,6) C 2000 Extra analysis fields (bands 10,11,12,13,14,15,16) C c---------------------------------------------------------------------------- C 04/20/2005 The bug with Bayesian clear sky probability fields c was fixed to provide the correct fields. c---------------------------------------------------------------------------- c c ================================= c Check whether this is a new image c ================================= if (image .eq. old_image) then ! always reset all output sat_id = hd(3) sat_ref= sav_tim else c ----------------------- c First, close old images c ----------------------- if (old_image .gt. 0) then ! old images were opened so ... call clsara(old_image) do i=1,nimages call clsara(allimage(i)) end do end if c --------------------- c Check 1: Files exist? c --------------------- call readd(image,hd) if (hd(1).ne.0) then write (*,'("Input AREA not present ?",2i6)') image, hd(1) return endif do i=1,nimages call readd(allimage(i),hd1) if (hd1(1).ne.0) then write (*,'("Input AREA not present ?",2i6)') allimage(i), hd1(1) return end if end do sat_id = hd(3) c ---------------------------------- c Check 2: First 10 words identical? c ---------------------------------- do j=1,nimages call readd(allimage(j),hd1) do i=1,10 if (hd(i) .ne. hd1(i)) then write (*,'("mismatch", 4i9)') allimage(j),i,hd(i),hd1(i) return end if end do end do c ------------- c Check 3: time c ------------- c >>> Step 1: Convert yyddd to ccyyddd <<< if (hd(4)/1000 .gt. 70) then i = 1900000 + hd(4) else i = 2000000 + hd(4) end if c >>> Step 2: ccyyddd & hhmmss ==> second since 1970 <<< if (mcdaytimetosec(i,hd(5), sat_ref).ne.0) then write (*,'("mcdaytimetosec 2 failed",8i9)') i, hd(5) return end if c >>> Step 3: Compare current & reference second <<< sat_ref = nint((sat_ref-ref_sec)/60.) sav_tim = sat_ref if (abs(sat_ref) .gt. 30) then write (*,'("inconsistent image time",2i11)') ref_sec,sat_ref return end if c ------------------- c Check 4: navigation c ------------------- if (nvset('AREA',image) .ne. 0) then write (*,'(" NVSET failed for AREA",i4.4,i4)') image return end if c ---------------------------- c Finally, open the new images c ---------------------------- cccc write(6,*) 'OPNARA ',image if (opnara(image).lt.0) then write (*,'("OPNARA failed for AREA",2i9)') image write(6,*) 'CLSARA ',image call clsara(image) endif c ---------------- c Update old_image c ---------------- old_image = image ! now we're sure it's a proper image end if c ---------- c Navigation c ---------- c >>> Buoy position in satellite coordinate <<< if (nv1eas(lat,lon,0., line,elem,0.) .ne. 0) then write (*,'("nv1eas failed for",2e15.5)') lat, lon return end if c >>> Buoy position in image coordinate <<< line = nint((line-hd(6))/hd(12)) elem = nint((elem-hd(7))/hd(13)) c >>> Buoy is on the edge on the image? <<< C 04/12/09 CPO increased boarder to allow for 7x7 box C if (line.lt.2 .or. line.gt.hd( 9)-1 .or. C & elem.lt.2 .or. elem.gt.hd(10)-1) then if (line.lt.4 .or. line.gt.hd( 9)-3 .or. & elem.lt.4 .or. elem.gt.hd(10)-3) then write (*,'("on image edge",2i19,2f9.0)') & hd(9),hd(10),line,elem return end if c >>> Get the angles <<< xin(1) = hd(3)*100000 + hd(4) ! ssyyddd xin(2) = hd(5)/10000 + mod(hd(5),10000)/6000. ! fractional hour xin(3) = lat ! lat (degree pos. N) xin(4) =-lon ! lon (degree pos. E) if (nv1opt('ANG ',xin,angles) .ne. 0) then write (*,'("nv1opt failed for",4e15.5)') xin return end if c ================ c Obtain GOES data c ================ c ----------------------- c Albedo: 0-1000 (0-100%) c ----------------------- call araopt(image,1,'SPAC',2) call araopt(image,1,'UNIT',lit('ALB ')) k = 0 do j=nint(line-1),nint(line+1) call redara(image,j-1,nint(elem-2),3,1,tb3) do i=1,3 k = k + 1 if (tb3(i).lt.0 .or. tb3(i).gt.1000) then buf(k,1) = -1 else buf(k,1) = tb3(i)*.1 end if end do end do k = 0 do j=nint(line-3),nint(line+3) call redara(image,j-1,nint(elem)-4,7,1,tb7) do i=1,7 k = k + 1 if (tb7(i).lt.0 .or. tb7(i).gt.1000) then ref_sat_1(k) = -1.0 else ref_sat_1(k) = tb7(i)*.1 end if end do end do do i=1,49 ref_sat_2(i) = 0.0 ref_sat_3(i) = 0.0 end do c ----------------------------- c Tb(2-5): 1500-3276 (150-327K) c ----------------------------- call araopt(image,1,'UNIT',lit('TEMP')) do kk=2,5 k = 0 do j=nint(line-1),nint(line+1) if ((sat_id .ge. 78).AND.(kk .eq. 5)) then call redara(image,j-1,nint(elem)-2,3,kk+1,tb3) else call redara(image,j-1,nint(elem)-2,3,kk,tb3) end if do i=1,3 k = k + 1 if (tb3(i).lt.1500.or.tb3(i).gt.3276) then buf(k,kk) = -1 else buf(k,kk) = tb3(i)*.1 end if end do end do end do c c********satellite BT reading******* c do kk=2,5 k = 0 do j=nint(line-3),nint(line+3) c*****1st if************************************************** if ((sat_id .ge. 78).AND.(kk .eq. 5)) then call redara(image,j-1,nint(elem)-4,7,kk+1,tb7) else call redara(image,j-1,nint(elem)-4,7,kk,tb7) end if c*****end of 1st if & end of area files reading**************** do i=1,7 k = k + 1 c c *****corrected by A.Frolov Apr.26, 2005 c if (tb7(i).lt.1500.or.tb7(i).gt.3276) then ! 2nd if c if (kk .eq. 2) then bt1_sat(k) = -1.0 end if c if (kk .eq. 3) then bt2_sat(k) = -1.0 end if c if (kk .eq. 4) then bt4_sat(k) = -1.0 end if c if (kk .eq. 5) then if (sat_id .lt. 78) then bt5_sat(k) = -1.0 else bt6_sat(k) = -1.0 end if end if c else if (kk .eq. 2) then bt1_sat(k) = tb7(i)*0.1 c call scale_back_SST(tb7(i),bt1_sat(k)) end if c if (kk .eq. 3) then bt2_sat(k) = tb7(i)*0.1 c call scale_back_SST(tb7(i),bt2_sat(k)) end if c if (kk .eq. 4) then bt4_sat(k) = tb7(i)*0.1 c call scale_back_SST(tb7(i),bt4_sat(k)) end if c if (kk .eq. 5) then if (sat_id .lt. 78) then bt5_sat(k) = tb7(i)*0.1 c call scale_back_SST(tb7(i),bt5_sat(k)) else bt6_sat(k) = tb7(i)*0.1 c call scale_back_SST(tb7(i),bt6_sat(k)) end if end if c end if !****end of the 2nd if c c*******end of A.Frolov's correction end do end do end do c --------------------- c SST: 7-255 (272-309K) c --------------------- cccc write(6,*) 'OPNARA ',jmage if (opnara(jmage).lt.0) then write (*,'("OPNARA failed for AREA",2i9)') jmage write(6,*) 'CLSARA ',jmage call clsara(jmage) endif call araopt(jmage,1,'SPAC',2) call araopt(jmage,1,'UNIT',lit('BRIT')) k = 0 do j=nint(line-1),nint(line+1) call redara(jmage,j-1,nint(elem-2),3,1,tb3) do i=1,3 k = k + 1 if (tb3(i) .lt. 0) tb3(i)=tb3(i)+255 if (tb3(i) .lt. 7) then buf(k,6) = tb3(i) else buf(k,6) = 270 + tb3(i)*.15 end if end do end do k = 0 do j=nint(line-3),nint(line+3) call redara(jmage,j-1,nint(elem)-4,7,1,tb7) do i=1,7 k = k + 1 if (tb7(i) .lt. 0) tb7(i)=tb7(i)+255 if (tb7(i) .lt. 7) then sst_prd(k) = tb7(i)*1.0 else sst_prd(k) = 270.0 + tb7(i)*0.15 end if end do end do c 11 Nov 2005 - O. Embury - Removed old "Average SST" code which was c corrupting image line/elem numbers c ----------------------------- c Average SST: 7-255 (272-309K) c ----------------------------- c line = ( 60-lat)*20 c elem = (180-lon)*20 c k = 0 c do j=-1,1 cccc read (9, rec=nint(line)+j) sst_buf c do i=nint(elem-1),nint(elem+1) c k = k + 1 c k = sst_buf(i) c if (sst_buf(i) .lt. 0) k=sst_buf(i)+255 c if (k .lt. 7) then c buf(k,7) = k c else c buf(k,7) = 270 + k*.15 c end if c end do c end do c ----------------------- c Statistics of GOES data c ----------------------- c >>> Clear the buffer <<< do j=1,7 do i=10,11 buf(i,j) = 0 end do end do c >>> Accumulation <<< clear = 0 do i=1,9 if (buf(i,6) .gt. 6) then clear = clear + 1 do j=1,7 buf(10,j) = buf(10,j) + buf(i,j) buf(11,j) = buf(11,j) + buf(i,j)*buf(i,j) end do end if end do c >>> Statistics <<< if (clear .gt. 0) then do j=1,7 buf(10,j) = buf(10,j)/clear buf(11,j) = sqrt(max(0.,buf(11,j)/clear-buf(10,j)**2)) end do end if call clsara(jmage) c ----------------------------- c all variables in a loop c --------------------------- index = 1 cccc write(6,*) 'OPNARA ', allimage(index) if ( opnara(allimage(index)).lt.0) then write (*,'("OPNARA failed for AREA",2i9)') allimage(index) write(6,*) 'CLSARA ',allimage(index) call clsara(allimage(index)) endif call araopt(allimage(index),1,'SPAC',2) call araopt(allimage(index),1,'UNIT',lit('BRIT')) c 30 Nov 2005 - O. Embury - Corrected scaling of SST error k = 0 do j=nint(line-3),nint(line+3) call redara(allimage(index),j-1,nint(elem)-4,7,2,tb7) do i=1,7 k = k + 1 sst_prd_err(k) = 0.3 + (tb7(i) * 0.0048) end do end do call clsara(allimage(index)) index = 2 do k=1,11 buf(k,8)=255 end do cccc write(6,*) 'OPNARA ', allimage(index) if ( opnara(allimage(index)).lt.0) then write (*,'("OPNARA failed for AREA",2i9)') allimage(index) write(6,*) 'CLSARA ',allimage(index) call clsara(allimage(index)) endif call araopt(allimage(index),1,'SPAC',2) call araopt(allimage(index),1,'UNIT',lit('BRIT')) k = 0 do j=nint(line-1),nint(line+1) call redara(allimage(index),j-1,nint(elem-2),3,3,tb3) do i=1,3 k = k + 1 if (tb3(i) .lt. 0) then tb3(i)=tb3(i)+255 endif buf(k,8) = tb3(i) end do end do c*********corrected by Frolov 04.19.2005 k = 0 do j=nint(line-3),nint(line+3) call redara(allimage(index),j-1,nint(elem)-4,7,3,tb7) do i=1,7 k = k + 1 pclear(k) = tb7(i) end do end do k = 0 do j=nint(line-3),nint(line+3) call redara(allimage(index),j-1,nint(elem)-4,7,4,tb7) do i=1,7 k = k + 1 pclear_revs(k) = tb7(i) end do end do k = 0 do j=nint(line-3),nint(line+3) call redara(allimage(index),j-1,nint(elem)-4,7,5,tb7) do i=1,7 k = k + 1 pclear_spec(k) = tb7(i) end do end do k = 0 do j=nint(line-3),nint(line+3) call redara(allimage(index),j-1,nint(elem)-4,7,6,tb7) do i=1,7 k = k + 1 pclear_text(k) = tb7(i) end do end do c*****convertion into probability**************************** c 30 Nov 2005 - O. Embury - removed scaling of SST error do i=1,49 call scale_back(NINT(pclear(i)),pclear(i)) call scale_back(NINT(pclear_revs(i)),pclear_revs(i)) call scale_back(NINT(pclear_spec(i)),pclear_spec(i)) call scale_back(NINT(pclear_text(i)),pclear_text(i)) c call scale_back(NINT(sst_prd_err(i)),sst_prd_err(i)) enddo c c**********end of the corrections****************************** call clsara(allimage(index)) index = 3 cccc write(6,*) 'OPNARA ', allimage(index) if ( opnara(allimage(index)).lt.0) then write (*,'("OPNARA failed for AREA",2i9)') allimage(index) write(6,*) 'CLSARA ',allimage(index) call clsara(allimage(index)) endif call araopt(allimage(index),1,'SPAC',2) call araopt(allimage(index),1,'UNIT',lit('BRIT')) c 30 Nov 2005 - O. Embury - Corrected unit scaling call redara(allimage(index),nint(line)-1,nint(elem)-1, & 1,10,tb) ref_ffm(1) = tb*1.0 call redara(allimage(index),nint(line)-1,nint(elem)-1, & 1,11,tb) dBT_scatt = tb * 0.01 call redara(allimage(index),nint(line)-1,nint(elem)-1, & 1,12,tb) dBT_glint = tb * 0.01 call redara(allimage(index),nint(line)-1,nint(elem)-1, & 1,13,tb) if ((tb .GE. 7) .AND. (tb .LE. 254)) then dBTdSST_ffm(1) = ((tb - 7.0)/(247.0 * 3.0)) + 0.6 else dBTdSST_ffm(1) = -999.99 endif call redara(allimage(index),nint(line)-1,nint(elem)-1, & 1,14,tb) if ((tb .GE. 7) .AND. (tb .LE. 254)) then dBTdTCWV_ffm(1) = (2.0*(tb - 7.0)/247.0) - 1.5 else dBTdTCWV_ffm(1) = -999.99 endif call redara(allimage(index),nint(line)-1,nint(elem)-1, & 1,15,tb) if ((tb .GE. 7) .AND. (tb .LE. 254)) then dBTdSST_ffm(4) = ((tb - 7.0)/(247.0 * 1.25)) + 0.1 else dBTdSST_ffm(4) = -999.99 endif call redara(allimage(index),nint(line)-1,nint(elem)-1, & 1,16,tb) if ((tb .GE. 7) .AND. (tb .LE. 254)) then dBTdTCWV_ffm(4) = (20.0*(tb - 7.0)/247.0) - 14.0 else dBTdTCWV_ffm(4) = -999.99 endif call redara(allimage(index),nint(line)-1,nint(elem)-1, & 1,22,tb) if ((tb .GE. 7) .AND. (tb .LE. 254)) then dBTdSST_ffm(5) = ((tb - 7.0)/(247.0 * 1.25)) + 0.1 else dBTdSST_ffm(5) = -999.99 endif call redara(allimage(index),nint(line)-1,nint(elem)-1, & 1,23,tb) if ((tb .GE. 7) .AND. (tb .LE. 254)) then dBTdTCWV_ffm(5) = (20.0*(tb - 7.0)/247.0) - 14.0 else dBTdTCWV_ffm(5) = -999.99 endif call redara(allimage(index),nint(line)-1,nint(elem)-1, & 1,25,tb) if ((tb .GE. 7) .AND. (tb .LE. 254)) then dBTdSST_ffm(2) = ((tb - 7.0)/(247.0 * 10.)) - 0.01 else dBTdSST_ffm(2) = -999.99 endif call redara(allimage(index),nint(line)-1,nint(elem)-1, & 1,26,tb) if ((tb .GE. 7) .AND. (tb .LE. 254)) then dBTdTCWV_ffm(2) = (13.0*(tb - 7.0)/247.0) - 15.0 else dBTdTCWV_ffm(2) = -999.99 endif cc call clsara(allimage(index)) c 30 Nov 2005 - O. Embury - Added FFM BTs call redara(allimage(index),nint(line)-1,nint(elem)-1, & 1,19,tb) if ((tb .GE. 1) .AND. (tb .LE. 254)) then bt_ffm(1) = 270 + tb*.15 else bt_ffm(1) = 0 endif call redara(allimage(index),nint(line)-1,nint(elem)-1, & 1,20,tb) if ((tb .GE. 1) .AND. (tb .LE. 254)) then bt_ffm(4) = 270 + tb*.15 else bt_ffm(4) = 0 endif call redara(allimage(index),nint(line)-1,nint(elem)-1, & 1,21,tb) if ((tb .GE. 1) .AND. (tb .LE. 254)) then ! bt_ffm(5) = 270 + tb*.15 pre-goes12 (12 micron) bt_ffm(5) = 247 + tb*.15 else bt_ffm(5) = 0 endif call redara(allimage(index),nint(line)-1,nint(elem)-1, & 1,24,tb) if ((tb .GE. 1) .AND. (tb .LE. 254)) then bt_ffm(2) = 230 + tb*.15 else bt_ffm(2) = 0 endif call clsara(allimage(index)) c ---------------------------- c Get stats for pclr and sigma c ---------------------------- c >>> Clear the buffer <<< do j=8,9 do i=10,11 buf(i,j) = 0 end do end do c >>> Accumulation <<< clear = 0 do i=1,9 if (buf(i,6) .gt. 6) then clear = clear + 1 do j=8,9 buf(10,j) = buf(10,j) + buf(i,j) buf(11,j) = buf(11,j) + buf(i,j)*buf(i,j) end do end if end do c >>> Statistics <<< if (clear .gt. 0) then c****this line was corrected by Frolov 09/06/2005 (responce on Owen Embury's comment) do j=8,9 buf(10,j) = buf(10,j)/clear buf(11,j) = sqrt(max(0.,buf(11,j)/clear-buf(10,j)**2)) end do end if c if (first_call) then c c Read in atmospheric profile data for temperature and relative humidity c from the file GRID1000 do ii = 1, 18 c changes on 03-11-2009 c Read in temperature grids c status = IGGET(1000,ii,maxgridpt,grid,nrows,ncols,hdgrid) status = IGGET(grid1,ii,maxgridpt,grid,nrows,ncols,hdgrid) do j = 1,360 do i = 1,181 lev_t(i,j,(19-ii)) = grid(i,j) end do end do c Read in relative humidity grids c status = IGGET(1000,(ii+20),maxgridpt,grid,nrows,ncols, status = IGGET(grid1,(ii+20),maxgridpt,grid,nrows,ncols, & hdgrid) do j = 1,360 do i = 1,181 lev_rh(i,j,(19-ii)) = grid(i,j) end do end do enddo c Write NWP timestamp c_out = mciydtocyd(hdgrid(4),icyd) if (mccydtodmy(icyd,day,month,year) .ne. 0) then write (*,'("mccydtodmy failed for",i7)') hdgrid(4) no_goes = .true. return end if NWP_day = day*1.0 NWP_month = month*1.0 C NWP_hour = nint(hdgrid(????)/10000.0)first_call c Read in the 1000mB altitude, z1, from the file GRID1000 c status = IGGET(1000,40,maxgridpt,z1,nrows,ncols,hdgrid) status = IGGET(grid1,40,maxgridpt,z1,nrows,ncols,hdgrid) c******corrected by Frolov 06/20/2005********** c Read in the surface parameters from the file GRID1001 c Read in total column water vapour, tcwv c status = IGGET(1001,1,maxgridpt,tcwv,nrows,ncols,hdgrid) status = IGGET(grid2,1,maxgridpt,tcwv,nrows,ncols,hdgrid) c Read in surface temperature, surf_t c status = IGGET(1001,2,maxgridpt,surf_t,nrows,ncols,hdgrid) status = IGGET(grid2,2,maxgridpt,surf_t,nrows,ncols,hdgrid) c Read in surface RH, surf_rh c status = IGGET(1001,3,maxgridpt,surf_rh,nrows,ncols,hdgrid) status = IGGET(grid2,3,maxgridpt,surf_rh,nrows,ncols,hdgrid) c Read in surface windspeed, surf_u c status = IGGET(1001,4,maxgridpt,surf_u,nrows,ncols,hdgrid) status = IGGET(grid2,4,maxgridpt,surf_u,nrows,ncols,hdgrid) c Read in surface windspeed, surf_v c status = IGGET(1001,5,maxgridpt,surf_v,nrows,ncols,hdgrid) status = IGGET(grid2,5,maxgridpt,surf_v,nrows,ncols,hdgrid) c******corrected by Frolov 01/05/2006********** c Read in surface SHF, surf_SHF c status = IGGET(1001,7,maxgridpt,surf_SHF,nrows,ncols,hdgrid) status = IGGET(grid2,7,maxgridpt,surf_SHF,nrows,ncols,hdgrid) c Read in surface LHF, surf_LHF c status = IGGET(1001,6,maxgridpt,surf_LHF,nrows,ncols,hdgrid) status = IGGET(grid2,6,maxgridpt,surf_LHF,nrows,ncols,hdgrid) c*****end of correction********************** c Read in surface pressure, surf_P c status = IGGET(1001,8,maxgridpt,surf_P,nrows,ncols,hdgrid) status = IGGET(grid2,8,maxgridpt,surf_P,nrows,ncols,hdgrid) c Read in temperature at 2 m, t_at2m c status = IGGET(1001,9,maxgridpt,t_at2m,nrows,ncols,hdgrid) status = IGGET(grid2,9,maxgridpt,t_at2m,nrows,ncols,hdgrid) c c ************end of the corrections******************************** c first_call = .false. NWP_monthc= NWP_month ! NWP analaysis month NWP_dayc= NWP_day ! NWP analysis day NWP_hourc= NWP_hour ! NWP analysis hour end if ! (first_call) NWP_month= NWP_monthc ! NWP analaysis month NWP_day = NWP_dayc ! NWP analysis day NWP_hour = NWP_hourc ! NWP analysis hour c ------------------------------- c Get 4 nearest model grid points c ------------------------------- nlat = 180 nlon = 360 c 11 Nov 20005 - O Embury - corrected lat/lon indices latlo = 90 - ifloor(lat) lathi = latlo + 1 lonlo = ifloor(360 - lon) + 1 c 14 May 2009 -- corrected lonlo in case of negative longitude lon lonlo = mod( lonlo-1, 360) + 1 lonhi = mod(lonlo, 360) + 1 latf = (90.0 - (latlo-1)) - lat lonf = 360.0 - (lonlo-1) - lon c ----------------------------- c Put values into MDB structure c ----------------------------- NWP_SST(1) = surf_t(latlo,lonlo)/100. NWP_SST(2) = surf_t(latlo,lonhi)/100. NWP_SST(3) = surf_t(lathi,lonhi)/100. NWP_SST(4) = surf_t(lathi,lonlo)/100. c O. Embury - 11 Nov 2005 - Bayesian code uses 1000mbar altitude c not sea level pressure. NWP_SLP(1) = z1(latlo,lonlo)/100. NWP_SLP(2) = z1(latlo,lonhi)/100. NWP_SLP(3) = z1(lathi,lonhi)/100. NWP_SLP(4) = z1(lathi,lonlo)/100. NWP_10mV(1) = surf_v(latlo,lonlo)/100. NWP_10mV(2) = surf_v(latlo,lonhi)/100. NWP_10mV(3) = surf_v(lathi,lonhi)/100. NWP_10mV(4) = surf_v(lathi,lonlo)/100. NWP_10mU(1) = surf_u(latlo,lonlo)/100. NWP_10mU(2) = surf_u(latlo,lonhi)/100. NWP_10mU(3) = surf_u(lathi,lonhi)/100. NWP_10mU(4) = surf_u(lathi,lonlo)/100. NWP_2mRH(1) = surf_rh(latlo,lonlo)/100. NWP_2mRH(2) = surf_rh(latlo,lonhi)/100. NWP_2mRH(3) = surf_rh(lathi,lonhi)/100. NWP_2mRH(4) = surf_rh(lathi,lonlo)/100. NWP_TCWV(1) = tcwv(latlo,lonlo)/100. NWP_TCWV(2) = tcwv(latlo,lonhi)/100. NWP_TCWV(3) = tcwv(lathi,lonhi)/100. NWP_TCWV(4) = tcwv(lathi,lonlo)/100. NWP_2mT(1) = t_at2m(latlo,lonlo)/100. NWP_2mT(2) = t_at2m(latlo,lonhi)/100. NWP_2mT(3) = t_at2m(lathi,lonhi)/100. NWP_2mT(4) = t_at2m(lathi,lonlo)/100. do ii = 1,18 NWP_RH_prof(1,ii) = lev_rh(latlo,lonlo,ii)/100. NWP_RH_prof(2,ii) = lev_rh(latlo,lonhi,ii)/100. NWP_RH_prof(3,ii) = lev_rh(lathi,lonhi,ii)/100. NWP_RH_prof(4,ii) = lev_rh(lathi,lonlo,ii)/100. NWP_T_prof(1,ii) = lev_t(latlo,lonlo,ii)/100. NWP_T_prof(2,ii) = lev_t(latlo,lonhi,ii)/100. NWP_T_prof(3,ii) = lev_t(lathi,lonhi,ii)/100. NWP_T_prof(4,ii) = lev_t(lathi,lonlo,ii)/100. end do c NWP_LHF_intrp(3) = surf_LHF(latlo,lonlo)/100. c NWP_SHF_intrp(3) = surf_SHF(latlo,lonlo)/100. c NWP_ST_intrp(3) = surf_t(latlo,lonlo)/100. c 11 Nov 2005 - O.Embury - Added interpolation of model fields. c 30 Nov 2005 - O.Embury - Added surface windspeed. For now use spare c fields in NWP_WS_intrp as temporary variables NWP_WS_intrp(1) = lonf * latf * surf_v(lathi,lonhi) + & lonf * (1.0-latf) * surf_v(latlo,lonhi) + & (1.0-lonf) * latf * surf_v(lathi,lonlo) + & (1.0-lonf) * (1.0-latf) * surf_v(latlo,lonlo) NWP_WS_intrp(2) = lonf * latf * surf_u(lathi,lonhi) + & lonf * (1.0-latf) * surf_u(latlo,lonhi) + & (1.0-lonf) * latf * surf_u(lathi,lonlo) + & (1.0-lonf) * (1.0-latf) * surf_u(latlo,lonlo) NWP_WS_intrp(3) = SQRT(NWP_WS_intrp(1)**2+NWP_WS_intrp(2)**2) NWP_WS_intrp(1) = 0.0 NWP_LHF_intrp(3) = lonf * latf * surf_LHF(lathi,lonhi) + & lonf * (1.0-latf) * surf_LHF(latlo,lonhi) + & (1.0-lonf) * latf * surf_LHF(lathi,lonlo) + & (1.0-lonf) * (1.0-latf) * surf_LHF(latlo,lonlo) NWP_SHF_intrp(3) = lonf * latf * surf_SHF(lathi,lonhi) + & lonf * (1.0-latf) * surf_SHF(latlo,lonhi) + & (1.0-lonf) * latf * surf_SHF(lathi,lonlo) + & (1.0-lonf) * (1.0-latf) * surf_SHF(latlo,lonlo) NWP_ST_intrp(3) = lonf * latf * surf_t(lathi,lonhi) + & lonf * (1.0-latf) * surf_t(latlo,lonhi) + & (1.0-lonf) * latf * surf_t(lathi,lonlo) + & (1.0-lonf) * (1.0-latf) * surf_t(latlo,lonlo) NWP_WS_intrp(3) = NWP_WS_intrp(3) / 100. c*******cleaning*** do ii = 1,2 NWP_WS_intrp(ii) = 0. NWP_WS_intrp(ii+3) = 0. enddo c********************* NWP_LHF_intrp(3) = NWP_LHF_intrp(3) / 100. NWP_SHF_intrp(3) = NWP_SHF_intrp(3) / 100. NWP_ST_intrp(3) = NWP_ST_intrp(3) / 100. c ---------------------------------- c Read grided SST from the area file c ----------------------------------- c RTG SST call ARSST(6320,lat,lon,ASST,ISTAT) IF( 0 .ne. ISTAT )THEN ADD_slot(1) = -999. ELSE ADD_slot(1) = ASST ENDIF c 30 Nov 2005 - O. Embury - Added extra outputs for future checking c of improved sunglint code. ADD_slot(2) = line ADD_slot(3) = elem ! ADD_slot(3) = float( hd(9)) ADD_slot(4) = float(hd(48)) c --------------------- c Return with GOES data c --------------------- no_goes = .false. return end c****************************************************** subroutine scale_back(inp1,output) c c PURPOSE: c Convert scaled count into a probability value). c c INPUT: c The output matchup file scaled cloudness probability values. c The data range is 2 - 255. c c OUTPUT: c Cloudness probability values. The data range is 0.002-0.9999. c c ============================================================= c implicit none integer inp1 ! scaled cloudness probability values real output ! cloudness probability values. c LOCAL VARIABLES: integer i real pcrl, E E = 40.546121 if ((inp1 .GE. 2) .AND.(inp1 .LE. 255)) then output = 1.002 - exp((2 - inp1)/E) else output = -1 endif return end c****************************************************** subroutine scale_back_SST(inp,output) c c PURPOSE: c Convert scaled SST count into a temperature1 value. c c INPUT: c The output matchup file scaled SST values. c The data range is 0 - 255. c c OUTPUT: c SST values. The data range is 270-350. c c ============================================================= c implicit none integer*2 inp ! scaled SST values real output ! SST values. c LOCAL VARIABLES: real E c ============================== E = 2700. if ((inp .GE. 0) .AND. (inp .LE. 255)) then output = (1.5*inp + E)/10. else output = -1 endif return end c******************************************************************** function ifloor(x) implicit none real x integer ifloor ifloor = int(x) if (x .LT. 0) then ifloor = ifloor - 1 endif return end