diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 74e871e..b50c332 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -107,6 +107,7 @@ endforeach() # different observation decoders set(decoders sst_viirs + sss_smap ) foreach(decoder ${decoders}) add_executable(OCN.decoder.${decoder}.x obs/decoder_${decoder}.f90) diff --git a/src/obs/decoder_sss_smap.f90 b/src/obs/decoder_sss_smap.f90 new file mode 100644 index 0000000..b8610d4 --- /dev/null +++ b/src/obs/decoder_sss_smap.f90 @@ -0,0 +1,171 @@ +PROGRAM decoder_sss_smap + USE common + USE params_model + USE vars_model + USE common_oceanmodel + !USE vars_obs + USE common_obs_oceanmodel + USE read_smap, ONLY: sss_data, read_jpl_smap_l2_sss_h5, & + read_jpl_smap_l3_sss_nc + + IMPLICIT NONE + + !----------------------------------------------------------------------------- + ! Command line inputs: + !----------------------------------------------------------------------------- + CHARACTER(slen) :: obsinfile='obsin.nc' !IN (default) + CHARACTER(slen) :: obsoutfile='obsout.dat' !OUT(default) + INTEGER :: obs_level = -1 !CDA: read level-2 smap if obs_level=2, + !CDA: level-3 smap if obs_level=3 + + !----------------------------------------------------------------------------- + ! Obs data arrays + !----------------------------------------------------------------------------- + TYPE(sss_data), ALLOCATABLE :: obs_data(:) + INTEGER :: nobs + CHARACTER(10) :: Syyyymmddhh = "YYYYMMDDHH" + ! 1234567890 + INTEGER :: delta_seconds = -999 ! max difference from Syyymmddhh + + REAL(r_size), ALLOCATABLE :: elem(:) + REAL(r_size), ALLOCATABLE :: rlon(:) + REAL(r_size), ALLOCATABLE :: rlat(:) + REAL(r_size), ALLOCATABLE :: rlev(:) + REAL(r_size), ALLOCATABLE :: odat(:) + REAL(r_size), ALLOCATABLE :: oerr(:) + REAL(r_size), ALLOCATABLE :: ohx(:) + REAL(r_size), ALLOCATABLE :: obhr(:) + INTEGER , ALLOCATABLE :: oqc(:) + + !----------------------------------------------------------------------------- + ! Miscellaneous + !----------------------------------------------------------------------------- + INTEGER :: n, i + + !----------------------------------------------------------------------------- + ! Debugging parameters + !----------------------------------------------------------------------------- + !STEVE: to adjust writing to output file + LOGICAL :: print1st = .true. + + !----------------------------------------------------------------------------- + ! Instantiations specific to this observation type: + !----------------------------------------------------------------------------- + INTEGER :: min_quality_level=5 ! CDA + + !----------------------------------------------------------------------------- + ! Initialize the common_oceanmodel module, and process command line options + !----------------------------------------------------------------------------- + CALL process_command_line !(get: -obsin -gues -obsout ) + + + !----------------------------------------------------------------------------- + ! Read observations from file + !----------------------------------------------------------------------------- + SELECT CASE(obs_level) + CASE(2) + if (delta_seconds>0) then + CALL read_jpl_smap_l2_sss_h5(trim(obsinfile), obs_data, nobs, & + Syyyymmddhh, delta_seconds) + else + CALL read_jpl_smap_l2_sss_h5(trim(obsinfile), obs_data, nobs) + end if + CASE(3) + CALL read_jpl_smap_l3_sss_nc(trim(obsinfile), obs_data, nobs) + CASE DEFAULT + WRITE(6,*) "[err] decoder_sss_smap.f90::Unsupported obs_level=", obs_level, & + ". obs_level supported should either be 2 or 3" + STOP (10) + END SELECT + + ALLOCATE( elem(nobs) ) + ALLOCATE( rlon(nobs) ) + ALLOCATE( rlat(nobs) ) + ALLOCATE( rlev(nobs) ) + ALLOCATE( odat(nobs) ) + ALLOCATE( oerr(nobs) ) + ALLOCATE( ohx(nobs) ) + ALLOCATE( oqc(nobs) ) + ALLOCATE( obhr(nobs) ) + + print *, "decoder_sss_smap.f90:: starting nobs = ", nobs + do i=1,nobs + elem(i) = obs_data(i)%typ + rlon(i) = obs_data(i)%x_grd(1) + rlat(i) = obs_data(i)%x_grd(2) + rlev(i) = 0.d0 + odat(i) = obs_data(i)%value + oerr(i) = obs_data(i)%oerr + ohx(i) = 0 + oqc(i) = 1 + obhr(i) = obs_data(i)%hour + enddo + DEALLOCATE(obs_data) + + if (print1st) then + i=1 + print *, "i, elem,rlon,rlat,rlev,odat,oerr,ohx,oqc,obhr = ", i, elem(i),rlon(i),rlat(i),rlev(i),odat(i),oerr(i),ohx(i),oqc(i),obhr(i) + i=nobs + print *, "i, elem,rlon,rlat,rlev,odat,oerr,ohx,oqc,obhr = ", i, elem(i),rlon(i),rlat(i),rlev(i),odat(i),oerr(i),ohx(i),oqc(i),obhr(i) + endif + + call write_obs3(trim(obsoutfile),nobs,elem,rlon,rlat,rlev, & + odat,oerr,obhr,oqc,qcflag_in=.true.) + + DEALLOCATE( elem,rlon,rlat,rlev,odat,oerr,ohx,oqc,obhr ) + +CONTAINS + +SUBROUTINE process_command_line +!=============================================================================== +! Process command line arguments +!=============================================================================== +IMPLICIT NONE +INTEGER, PARAMETER :: slen2=1024 +CHARACTER(slen2) :: arg1,arg2 +INTEGER :: i, ierr +INTEGER, DIMENSION(3) :: values + +! STEVE: add input error handling! +! inputs are in the format "-x xxx" +do i=1,COMMAND_ARGUMENT_COUNT(),2 + CALL GET_COMMAND_ARGUMENT(i,arg1) + PRINT *, "In obsop_sst_podaac.f90::" + PRINT *, "Argument ", i, " = ",TRIM(arg1) + + select case (arg1) + case('-obsin') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + obsinfile = arg2 + case('-obsout') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + obsoutfile = arg2 + case('-minqc') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + read (arg2,*) min_quality_level + case('-qcyyyymmddhh') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + !read (arg2,*) min_quality_level + Syyyymmddhh = arg2 + case('-maxdt') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + read (arg2,*) delta_seconds + case('-olevel') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + read (arg2,*) obs_level + case default + PRINT *, "ERROR: option is not supported: ", arg1 + PRINT *, "(with value : ", trim(arg2), " )" + stop 1 + end select +enddo + +END SUBROUTINE process_command_line + +END PROGRAM decoder_sss_smap diff --git a/src/obs/obsop_sss.f90 b/src/obs/obsop_sss.f90 index 5e808b3..2872f69 100644 --- a/src/obs/obsop_sss.f90 +++ b/src/obs/obsop_sss.f90 @@ -18,7 +18,6 @@ PROGRAM obsop_sss !----------------------------------------------------------------------------- ! Command line inputs: !----------------------------------------------------------------------------- - CHARACTER(slen) :: navinfile='navin.nc' !IN (default) CHARACTER(slen) :: obsinfile='obsin.nc' !IN (default) CHARACTER(slen) :: guesfile='gues' !IN (default) i.e. prefix to '.ocean_temp_salt.res.nc' CHARACTER(slen) :: obsoutfile='obsout.dat' !OUT(default) @@ -30,6 +29,9 @@ PROGRAM obsop_sss !----------------------------------------------------------------------------- ! Obs data arrays !----------------------------------------------------------------------------- + LOGICAL :: BINARY_INPUT = .FALSE. + CHARACTER(10) :: Syyyymmddhh = "YYYYMMDDHH" + INTEGER :: delta_seconds = -999 TYPE(sss_data), ALLOCATABLE :: obs_data(:) REAL(r_size), ALLOCATABLE :: elem(:) @@ -100,45 +102,74 @@ PROGRAM obsop_sss CALL set_common_oceanmodel CALL process_command_line !(get: -obsin -gues -obsout ) - ALLOCATE(superobs(nlon,nlat),delta(nlon,nlat),M2(nlon,nlat),supercnt(nlon,nlat)) - !----------------------------------------------------------------------------- ! Read observations from file !----------------------------------------------------------------------------- - SELECT CASE(obs_level) - CASE(2) - CALL read_jpl_smap_l2_sss_h5(trim(obsinfile), obs_data, nobs) - CASE(3) - CALL read_jpl_smap_l3_sss_nc(trim(obsinfile), obs_data, nobs) - CASE DEFAULT - WRITE(6,*) "[err] obsop_sss.f90::Unsupported obs_level=", obs_level, & - ". obs_level supported should either be 2 or 3" - STOP (10) - END SELECT - - ALLOCATE( elem(nobs) ) - ALLOCATE( rlon(nobs) ) - ALLOCATE( rlat(nobs) ) - ALLOCATE( rlev(nobs) ) - ALLOCATE( odat(nobs) ) - ALLOCATE( oerr(nobs) ) - ALLOCATE( ohx(nobs) ) - ALLOCATE( oqc(nobs) ) - ALLOCATE( obhr(nobs) ) - - print *, "obsop_sss.f90:: starting nobs = ", nobs - do i=1,nobs - elem(i) = obs_data(i)%typ - rlon(i) = obs_data(i)%x_grd(1) - rlat(i) = obs_data(i)%x_grd(2) -! rlev(i) = obs_data(i)%x_grd(3) !not applicable for SSS - odat(i) = obs_data(i)%value - oerr(i) = obs_data(i)%oerr - ohx(i) = 0 - oqc(i) = 0 - obhr(i) = obs_data(i)%hour - enddo - DEALLOCATE(obs_data) + call create_empty_obsfile(trim(obsoutfile)) + + if (BINARY_INPUT) then ![processed binary input] + print *, "obsop_sss.f90:: reading binary file=",trim(obsinfile) + CALL get_nobs(trim(obsinfile),8,nobs,errIfNoObs=.false.) + + if (nobs<=0) then + print*, "obsop_sss.f90: nobs=",nobs,"<=0. Exit now..." + STOP (0) + endif + + ALLOCATE( elem(nobs) ) + ALLOCATE( rlon(nobs) ) + ALLOCATE( rlat(nobs) ) + ALLOCATE( rlev(nobs) ) + ALLOCATE( odat(nobs) ) + ALLOCATE( oerr(nobs) ) + ALLOCATE( ohx(nobs) ) + ALLOCATE( oqc(nobs) ) + ALLOCATE( obhr(nobs) ) + + CALL read_obs3(trim(obsinfile), nobs, elem, rlon, rlat, rlev, odat, oerr, obhr) + oqc(:) = 0 ! now set all obs as bad. + ohx(:) = 0.d0 + else ! [nc input] + SELECT CASE(obs_level) + CASE(2) + if (delta_seconds>0) then + CALL read_jpl_smap_l2_sss_h5(trim(obsinfile), obs_data, nobs, & + Syyyymmddhh, delta_seconds) + else + CALL read_jpl_smap_l2_sss_h5(trim(obsinfile), obs_data, nobs) + end if + CASE(3) + CALL read_jpl_smap_l3_sss_nc(trim(obsinfile), obs_data, nobs) + CASE DEFAULT + WRITE(6,*) "[err] obsop_sss.f90::Unsupported obs_level=", obs_level, & + ". obs_level supported should either be 2 or 3" + STOP (10) + END SELECT + + ALLOCATE( elem(nobs) ) + ALLOCATE( rlon(nobs) ) + ALLOCATE( rlat(nobs) ) + ALLOCATE( rlev(nobs) ) + ALLOCATE( odat(nobs) ) + ALLOCATE( oerr(nobs) ) + ALLOCATE( ohx(nobs) ) + ALLOCATE( oqc(nobs) ) + ALLOCATE( obhr(nobs) ) + + print *, "obsop_sss.f90:: starting nobs = ", nobs + do i=1,nobs + elem(i) = obs_data(i)%typ + rlon(i) = obs_data(i)%x_grd(1) + rlat(i) = obs_data(i)%x_grd(2) + rlev(i) = 0.d0 + odat(i) = obs_data(i)%value + oerr(i) = obs_data(i)%oerr + ohx(i) = 0 + oqc(i) = 0 + obhr(i) = obs_data(i)%hour + enddo + DEALLOCATE(obs_data) + end if if (print1st) then print *, "elem(1),rlon(1),rlat(1),obhr(1) = ", elem(1),rlon(1),rlat(1),obhr(1) @@ -157,27 +188,36 @@ PROGRAM obsop_sss ! Bin the obs and estimate the obs error based on bin standard deviations !----------------------------------------------------------------------------- if (DO_SUPEROBS) then + ALLOCATE(superobs(nlon,nlat),delta(nlon,nlat),M2(nlon,nlat),supercnt(nlon,nlat)) + superobs=0.0; delta = 0.0; M2 = 0.0; supercnt = 0 + print *, "Computing superobs..." - supercnt = 0 - superobs = 0.0d0 !STEVE: (ISSUE) this may have problems in the arctic for the tripolar grid, ! or for any irregular grid in general. do n=1,nobs ! for each ob, ! if (dodebug1) print *, "n = ", n - ! Scan the longitudes - do i=1,nlon-1 - if (lon(i+1) > rlon(n)) exit - enddo - ! Scan the latitudes - do j=1,nlat-1 - if (lat(j+1) > rlat(n)) exit - enddo + if (DO_REMOVE_65N .and. rlat(n) > 65) CYCLE + !! Scan the longitudes + !do i=1,nlon-1 + ! if (lon(i+1) > rlon(n)) exit + !enddo + !! Scan the latitudes + !do j=1,nlat-1 + ! if (lat(j+1) > rlat(n)) exit + !enddo + !CALL phys2ijk(elem(n),rlon(n),rlat(n),rlev(n),ri,rj,rk) !(OCEAN) + CALL phys2ij_nearest(rlon(n),rlat(n),ri,rj) ! OCEAN + + if (CEILING(ri) < 1 .OR. nlon < CEILING(ri)) CYCLE + if (CEILING(rj) < 1 .OR. nlat < CEILING(rj)) CYCLE + + i = NINT(ri); j = NINT(rj) supercnt(i,j) = supercnt(i,j) + 1 - delta(i,j) = odat(n) - superobs(i,j) + delta(i,j) = odat(n) - superobs(i,j) superobs(i,j) = superobs(i,j) + delta(i,j)/supercnt(i,j) - M2(i,j) = M2(i,j) + delta(i,j)*(odat(n) - superobs(i,j)) + M2(i,j) = M2(i,j) + delta(i,j)*(odat(n) - superobs(i,j)) ! if (dodebug1) print *, "supercnt(",i,",",j,") = ", supercnt(i,j) enddo !"superobs" contains the mean @@ -187,14 +227,14 @@ PROGRAM obsop_sss idx=0 do j=1,nlat-1 do i=1,nlon-1 - if (supercnt(i,j) > 1) then + if (supercnt(i,j) > 1 .and. wet(i,j)>0.5) then ! only retain ocn pts defined by MOM idx = idx+1 if (dodebug1) print *, "idx = ", idx odat(idx) = superobs(i,j) oerr(idx) = obserr_scaling*(min_oerr + SQRT(M2(i,j))) - rlon(idx) = (lon(i+1)+lon(i))/2.0d0 - rlat(idx) = (lat(j+1)+lat(j))/2.0d0 - rlev(idx) = 0 + rlon(idx) = lon(i) !(lon(i+1)+lon(i))/2.0d0 + rlat(idx) = lat(j) !(lat(j+1)+lat(j))/2.0d0 + rlev(idx) = 0.d0 elem(idx) = id_sss_obs if (dodebug1) print *, "odat(idx) = ", odat(idx) if (dodebug1) print *, "oerr(idx) = ", oerr(idx) @@ -204,7 +244,7 @@ PROGRAM obsop_sss endif enddo enddo - print *, "DO_SUPEROBS:: superobs reducing from ",nobs," to ",idx, " observations." + print *, "DO_SUPEROBS:: superobs reducing from ",nobs," to ",idx, " observations." print *, "with min obs error = ", MINVAL(oerr(1:idx)) print *, "with max obs error = ", MAXVAL(oerr(1:idx)) nobs = idx @@ -395,7 +435,16 @@ PROGRAM obsop_sss !----------------------------------------------------------------------------- ! Write the observations and their associated innovations to file !----------------------------------------------------------------------------- - CALL write_obs2(obsoutfile,nobs,elem,rlon,rlat,rlev,odat,oerr,ohx,oqc,obhr) + CALL write_obs2(obsoutfile,nobs,elem(1:nobs), & + rlon(1:nobs), & + rlat(1:nobs), & + rlev(1:nobs), & + odat(1:nobs), & + oerr(1:nobs), & + ohx(1:nobs), & + oqc(1:nobs), & + obhr(1:nobs), & + lappend_in = .true.) DEALLOCATE( elem,rlon,rlat,rlev,odat,oerr,ohx,oqc,obhr,v3d,v2d ) @@ -423,10 +472,6 @@ SUBROUTINE process_command_line CALL GET_COMMAND_ARGUMENT(i+1,arg2) PRINT *, "Argument ", i+1, " = ",TRIM(arg2) obsinfile = arg2 - case('-navin') - CALL GET_COMMAND_ARGUMENT(i+1,arg2) - PRINT *, "Argument ", i+1, " = ",TRIM(arg2) - navinfile = arg2 case('-gues') CALL GET_COMMAND_ARGUMENT(i+1,arg2) PRINT *, "Argument ", i+1, " = ",TRIM(arg2) @@ -463,6 +508,15 @@ SUBROUTINE process_command_line CALL GET_COMMAND_ARGUMENT(i+1,arg2) PRINT *, "Argument ", i+1, " = ",TRIM(arg2) read (arg2,*) obs_level + case('-qcyyyymmddhh') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + !read (arg2,*) min_quality_level + Syyyymmddhh = arg2 + case('-maxdt') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + read (arg2,*) delta_seconds case default PRINT *, "ERROR: option is not supported: ", arg1 PRINT *, "(with value : ", trim(arg2), " )" diff --git a/src/obs/read_smap.f90 b/src/obs/read_smap.f90 index ac721ba..ea1ad47 100644 --- a/src/obs/read_smap.f90 +++ b/src/obs/read_smap.f90 @@ -205,7 +205,7 @@ SUBROUTINE read_jpl_smap_l3_sss_nc(obsinfile, obs_data, nobs) end if end do end do - CALL inspect_obs_data(obs_data,subname="read_jpl_smap_l2_sss_h5") + CALL inspect_obs_data(obs_data,subname="read_jpl_smap_l3_sss_h5") if (n/=nobs) then WRITE(6,*) "[err] read_jpl_smap_l2_sss_h5::n/=nobs: n, nobs=", n, nobs STOP (27) @@ -220,7 +220,7 @@ SUBROUTINE read_jpl_smap_l3_sss_nc(obsinfile, obs_data, nobs) ! - Data documentation webpage: ! https://podaac.jpl.nasa.gov/dataset/SMAP_JPL_L2B_SSS_CAP_V5 !------------------------------------------------------------------------------- -SUBROUTINE read_jpl_smap_l2_sss_h5(obsinfile, obs_data, nobs) +SUBROUTINE read_jpl_smap_l2_sss_h5(obsinfile, obs_data, nobs, Syyyymmddhh, delta_seconds) USE m_h5io, ONLY: h5_get_fid, h5_close_fid, h5_rdvarshp, h5_rdvar2d, & h5_rdvar1d, h5_rdatt, & HID_T, HSIZE_T, i4, r4 @@ -228,22 +228,31 @@ SUBROUTINE read_jpl_smap_l2_sss_h5(obsinfile, obs_data, nobs) CHARACTER(*),INTENT(IN) :: obsinfile TYPE(sss_data),ALLOCATABLE,INTENT(INOUT) :: obs_data(:) INTEGER, INTENT(OUT) :: nobs + CHARACTER(10),INTENT(IN), OPTIONAL :: Syyyymmddhh + !1234567890 + INTEGER, INTENT(IN), OPTIONAL :: delta_seconds ! qc'ed if abs(obstime-Syyyymmddhh)>delta_seconds INTEGER(HID_T) :: fid INTEGER(HSIZE_T),allocatable :: varsizes(:) ! (2) REAL(r_size),ALLOCATABLE :: alon2d(:,:), alat2d(:,:) ! (nlines,npixels) REAL(r_size),ALLOCATABLE :: sea_surface_salinity(:,:), stde(:,:) ! (nlines,npixels) REAL(r_size),ALLOCATABLE :: sss_time(:) ! (nlines) + INTEGER,ALLOCATABLE :: sss_time_in_seconds_since19780101(:) ! (nlines) INTEGER(i4), ALLOCATABLE :: quality_flag(:,:) ! (nlines,npixels), need to ensure bits>16. Use 32bits here. LOGICAL, ALLOCATABLE :: valid(:,:) ! .true. if no missing info at this grid by checking different FillValues in ncfile REAL(r4) :: r4FillValue ! need to ensure bits=32 INTEGER(i4) :: i4FillValue ! need to ensure bits>16. Use 32bits here - INTEGER :: sss_sdate(8), sss_edate(8), sdate(8) + INTEGER :: sss_sdate(8), sss_edate(8) + INTEGER :: sdate(8), sdate_in_seconds_since19780101, sdate_in_min_since19780101 + INTEGER :: qc_sdate(8), qc_edate(8) + INTEGER :: qcdate(5), qcdate_in_seconds_since19780101, qcdate_in_min_since19780101 REAL :: tdelta(5) INTEGER :: nlines, npixels, i, j, n, istat LOGICAL :: dodebug = .true. INTEGER,PARAMETER :: QUAL_FLAG_SSS_USABLE_GOOD = 0 INTEGER,PARAMETER :: QUAL_FLAG_SSS_USABLE_BAD = 1 + INTEGER,PARAMETER :: QUAL_FLAG_SSS_HAS_LAND = 1 + INTEGER,PARAMETER :: QUAL_FLAG_SSS_HAS_ICE = 1 !------------------------------------------------------------------------------- ! Open the hdf5 file @@ -258,26 +267,74 @@ SUBROUTINE read_jpl_smap_l2_sss_h5(obsinfile, obs_data, nobs) nlines = varsizes(1); npixels = varsizes(2) WRITE(6,*) "nlines, npixels=", nlines, npixels + ALLOCATE(alon2d(nlines,npixels), alat2d(nlines,npixels)) + ALLOCATE(valid(nlines,npixels)) + valid = .true. !------------------------------------------------------------------------------- ! Read row time !------------------------------------------------------------------------------- sdate = [2015, 1, 1, 0, 0, 0, 0, 0] ! YYYY/MON/DAY/TZONE/HR/MIN/SEC/MSEC + CALL W3FS21(sdate(1:5), sdate_in_min_since19780101) + sdate_in_seconds_since19780101 = 60 * sdate_in_min_since19780101 + ALLOCATE(sss_time(nlines)) - CALL h5_rdvar1d(fid, "/row_time", sss_time) + ALLOCATE(sss_time_in_seconds_since19780101(nlines)) + + CALL h5_rdvar1d(fid, "/row_time", sss_time) ! elapse seconds since sdate tdelta = [0.0,0.0,0.0,REAL(sss_time(1)),0.0] ! DAY/HR/MIN/SEC/MSEC - call w3movdat(tdelta, sdate, sss_sdate) + CALL w3movdat(tdelta, sdate, sss_sdate) WRITE(6,*) "[msg] read_jpl_smap_l2_sss_h5::start date=", sss_sdate(1:3),sss_sdate(5:7) tdelta = [0.0,0.0,0.0,REAL(sss_time(nlines)),0.0] ! DAY/HR/MIN/SEC/MSEC - call w3movdat(tdelta, sdate, sss_edate) + CALL w3movdat(tdelta, sdate, sss_edate) WRITE(6,*) "[msg] read_jpl_smap_l2_sss_h5::end date=", sss_edate(1:3),sss_edate(5:7) + sss_time_in_seconds_since19780101(:) = sdate_in_seconds_since19780101 + & + CEILING(sss_time(:)) + WRITE(6,*) "[msg] read_jpl_smap_l2_sss_h5::sss_dtime_since19780101: min, max=", & + minval(sss_time_in_seconds_since19780101), & + maxval(sss_time_in_seconds_since19780101) + + ! time filtering + if (PRESENT(Syyyymmddhh) .and. PRESENT(delta_seconds)) then + ! convert qc center time + qcdate = 0 + READ(Syyyymmddhh,"(I4,I2,I2,I2)") qcdate(1),qcdate(2),qcdate(3),qcdate(4) + WRITE(6,*) "[msg] read_jpl_smap_l2_sss_h5::Syyyymmddhh=", Syyyymmddhh, & + "qcdate(:)=",qcdate(:) + call W3FS21(qcdate,qcdate_in_min_since19780101) + qcdate_in_seconds_since19780101 = 60 * qcdate_in_min_since19780101 + WRITE(6,*) "[msg] read_jpl_smap_l2_sss_h5::qcdate_in_seconds_since19780101, delta_seconds=", & + qcdate_in_seconds_since19780101, delta_seconds + + sdate = [1978, 1, 1, 0, 0, 0, 0, 0] ! YYYY/MON/DAY/TZONE/HR/MIN/SEC/MSEC + tdelta = [0.0, & + REAL((qcdate_in_seconds_since19780101-delta_seconds)/3600), & + 0.0, & + REAL(mod(qcdate_in_seconds_since19780101-delta_seconds,3600)), & + 0.0] + CALL w3movdat(tdelta, sdate, qc_sdate) + tdelta = [0.0,& + REAL((qcdate_in_seconds_since19780101+delta_seconds)/3600), & + 0.0,& + REAL(mod(qcdate_in_seconds_since19780101+delta_seconds,3600)),& + 0.0] + CALL w3movdat(tdelta, sdate, qc_edate) + WRITE(6,*) "[DEBUG]: qc_sdate=", qc_sdate(1:3), qc_sdate(5:7) + WRITE(6,*) "[DEBUG]: qc_edate=", qc_edate(1:3), qc_edate(5:7) + + ! calculate time difference + do i = 1, nlines + if (ABS(sss_time_in_seconds_since19780101(i) - & + qcdate_in_seconds_since19780101) > delta_seconds) then + valid(i,:) = .false. + end if + end do + end if + + !------------------------------------------------------------------------------- ! Read lat & lon info !------------------------------------------------------------------------------- - ALLOCATE(alon2d(nlines,npixels), alat2d(nlines,npixels)) - ALLOCATE(valid(nlines,npixels)) - valid = .true. - CALL h5_rdvar2d(fid, "/lon", alon2d) CALL h5_rdatt(fid, "/lon", "_FillValue", r4FillValue) if (dodebug) WRITE(6,*) "_FillValue=", NINT(r4FillValue) @@ -306,13 +363,34 @@ SUBROUTINE read_jpl_smap_l2_sss_h5(obsinfile, obs_data, nobs) where (NINT(sea_surface_salinity)==NINT(r4FillValue)) valid = .false. end where + CALL h5_rdatt(fid, "/smap_sss", "valid_max", r4FillValue) + WRITE(6,*) "valid_max_yo=", r4FillValue + where(sea_surface_salinity > r4FillValue) + valid = .false. + end where + CALL h5_rdatt(fid, "/smap_sss", "valid_min", r4FillValue) + WRITE(6,*) "valid_min_yo=", r4FillValue + where( sea_surface_salinity < r4FillValue) + valid = .false. + endwhere + WRITE(6,*) "[msg] read_jpl_smap_l2_sss_h5::smap_sss: min, max=", & minval(sea_surface_salinity, mask=valid), & maxval(sea_surface_salinity, mask=valid) CALL h5_rdvar2d(fid, "/smap_sss_uncertainty", stde) CALL h5_rdatt(fid, "/smap_sss_uncertainty", "_FillValue", r4FillValue) - where (NINT(stde)==NINT(r4FillValue)) + where (NINT(stde)==NINT(r4FillValue)) ! + valid = .false. + end where + CALL h5_rdatt(fid, "/smap_sss_uncertainty", "valid_max", r4FillValue) + WRITE(6,*) "valid_max_error=", r4FillValue + where(stde > r4FillValue) + valid = .false. + end where + CALL h5_rdatt(fid, "/smap_sss_uncertainty", "valid_min", r4FillValue) + WRITE(6,*) "valid_min_error=", r4FillValue + where(stde < r4FillValue) valid = .false. end where WRITE(6,*) "[msg] read_jpl_smap_l2_sss_h5::smap_sss_uncertainty: min, max=", & @@ -364,7 +442,10 @@ SUBROUTINE read_jpl_smap_l2_sss_h5(obsinfile, obs_data, nobs) CALL h5_rdvar2d(fid, "/quality_flag", quality_flag) CALL h5_rdatt(fid, "/quality_flag", "_FillValue", i4FillValue) if (dodebug) WRITE(6,*) "i4FillValue", i4FillValue - where (quality_flag == i4FillValue .or. IBITS(quality_flag,0,1)==QUAL_FLAG_SSS_USABLE_BAD) ! or BIT 0 QC FLAG is bad + where (quality_flag == i4FillValue .or. & + IBITS(quality_flag,0,1)==QUAL_FLAG_SSS_USABLE_BAD .or. & ! overall bad + IBITS(quality_flag,7,1)==QUAL_FLAG_SSS_HAS_LAND .or. & ! has land + IBITS(quality_flag,8,1)==QUAL_FLAG_SSS_HAS_ICE ) ! ! has ice valid = .false. end where WRITE(6,*) "[msg] read_jpl_smap_l2_sss_h5::quality_flag: min, max=", & @@ -398,7 +479,7 @@ SUBROUTINE read_jpl_smap_l2_sss_h5(obsinfile, obs_data, nobs) obs_data(n)%typ = id_sss_obs obs_data(n)%x_grd(1) = alon2d(i,j) obs_data(n)%x_grd(2) = alat2d(i,j) - obs_data(n)%hour = 0.0 ! [FIXME]: need to determine later what hours to use here + obs_data(n)%hour = sss_time_in_seconds_since19780101(i)/3600. obs_data(n)%value = sea_surface_salinity(i,j) obs_data(n)%oerr = stde(i,j) end if