diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 9060fa0..eadd00a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -59,6 +59,7 @@ set (SRCS obs/read_avhrr_pathfinder.f90 obs/read_geostationary.f90 obs/read_viirs.f90 + obs/read_sst_acspo.f90 obs/read_smap.f90 obs/read_aviso_adt.f90 obs/read_ice_txt.f90 @@ -91,6 +92,7 @@ set(operators sst sst_geostationary sst_viirs + sst_acspo sss tprof sprof @@ -107,6 +109,7 @@ endforeach() # different observation decoders set(decoders sst_viirs + sst_acspo sss_smap ) foreach(decoder ${decoders}) diff --git a/src/letkf/letkf_obs.f90.new b/src/letkf/letkf_obs.f90.new index d3015ad..cf3b81e 100644 --- a/src/letkf/letkf_obs.f90.new +++ b/src/letkf/letkf_obs.f90.new @@ -102,8 +102,8 @@ SUBROUTINE set_letkf_obs !CDA: related to AOEI REAL(r_size) :: oerr_aoei REAL(r_size) :: varhdxf ! ensemble spread for h(x) - INTEGER :: cnt_aoei_sst, cnt_aoei_sss - REAL(r_size) :: max_infl_sst, max_infl_sss + INTEGER :: cnt_aoei_sst, cnt_aoei_sss, cnt_aoei_sfcwnd + REAL(r_size) :: max_infl_sst, max_infl_sss, max_infl_sfcwnd WRITE(6,'(A)') 'Hello from set_letkf_obs' @@ -329,6 +329,9 @@ quality_control : if (.true.) then elseif (NINT(obselm(n)) .eq. id_sss_obs) then cnt_aoei_sss = cnt_aoei_sss + 1 if (max_infl_sss < oerr_aoei / obserr(n)) max_infl_sss = oerr_aoei / obserr(n) + elseif (NINT(obselm(n)) .eq. id_sfcwnd_obs) then + cnt_aoei_sfcwnd = cnt_aoei_sfcwnd + 1 + if (max_infl_sfcwnd < oerr_aoei / obserr(n)) max_infl_sfcwnd = oerr_aoei / obserr(n) endif obserr(n) = oerr_aoei endif @@ -458,6 +461,7 @@ endif quality_control WRITE(6,*) "cnt_aoei_sst = ", cnt_aoei_sst, "max_infl_sst=", max_infl_sst WRITE(6,*) "cnt_aoei_sss = ", cnt_aoei_sss, "max_infl_sss=", max_infl_sss + WRITE(6,*) "cnt_aoei_sfcwnd = ", cnt_aoei_sfcwnd, "max_infl_sfcwnd=", max_infl_sfcwnd WRITE(6,*) "QC check:" WRITE(6,*) "cnt_qc_sstval = ", cnt_qc_sstval diff --git a/src/letkf/vars_letkf.f90 b/src/letkf/vars_letkf.f90 index 2d4b75a..81df2fd 100644 --- a/src/letkf/vars_letkf.f90 +++ b/src/letkf/vars_letkf.f90 @@ -41,7 +41,8 @@ MODULE vars_letkf !STEVE: I think this would be better in "vars_obs.f90", maybe. (ISSUE) SUBROUTINE get_iobs(oelm,iobs) USE params_obs, ONLY: id_u_obs, id_v_obs, id_t_obs, id_s_obs, & - id_ssh_obs, id_ssh_obs, id_sst_obs, id_sss_obs, id_eta_obs + id_ssh_obs, id_ssh_obs, id_sst_obs, id_sss_obs, id_eta_obs, & + id_sfcwnd_obs REAL(r_size), INTENT(IN) :: oelm INTEGER, INTENT(OUT) :: iobs !--------------------------------------------------------------------------- @@ -64,6 +65,8 @@ SUBROUTINE get_iobs(oelm,iobs) iobs=7 CASE(id_eta_obs) !(OCEAN) iobs=8 + CASE(id_sfcwnd_obs) ! (ATM) + iobs=9 CASE DEFAULT WRITE(6,*) "vars_letkf.f90 :: there is no variable localization for obs-type :: ", oelm WRITE(6,*) "vars_letkf.f90 :: FATAL ERROR, exiting..." diff --git a/src/model_specific/mom6/common_obs_mom6.f90.kdtree b/src/model_specific/mom6/common_obs_mom6.f90.kdtree index e811640..a9e4046 100644 --- a/src/model_specific/mom6/common_obs_mom6.f90.kdtree +++ b/src/model_specific/mom6/common_obs_mom6.f90.kdtree @@ -487,6 +487,9 @@ SUBROUTINE monit_dep(nn,elm,dep,qc) REAL(r_size) :: rmse_u,rmse_v,rmse_t,rmse_s,rmse_ssh,rmse_eta,rmse_sst,rmse_sss !(OCEAN) REAL(r_size) :: bias_u,bias_v,bias_t,bias_s,bias_ssh,bias_eta,bias_sst,bias_sss !(OCEAN) INTEGER :: n,iu,iv,it,is,issh,ieta,isst,isss !(OCEAN) + REAL(r_size) :: rmse_sfcwnd ! (ATM) + REAL(r_size) :: bias_sfcwnd ! (ATM) + INTEGER :: isfcwnd ! (ATM) rmse_u = 0.0d0 rmse_v = 0.0d0 @@ -496,6 +499,9 @@ SUBROUTINE monit_dep(nn,elm,dep,qc) rmse_eta = 0.0d0 !(OCEAN) rmse_sst = 0.0d0 !(OCEAN) rmse_sss = 0.0d0 !(OCEAN) + + rmse_sfcwnd = 0.d0 ! (ATM) + bias_u = 0.0d0 bias_v = 0.0d0 bias_t = 0.0d0 @@ -504,6 +510,9 @@ SUBROUTINE monit_dep(nn,elm,dep,qc) bias_eta = 0.0d0 !(OCEAN) bias_sst = 0.0d0 !(OCEAN) bias_sss = 0.0d0 !(OCEAN) + + bias_sfcwnd = 0.d0 ! (ATM) + iu = 0 iv = 0 it = 0 @@ -512,6 +521,9 @@ SUBROUTINE monit_dep(nn,elm,dep,qc) ieta = 0 !(OCEAN) isst = 0 !(OCEAN) isss = 0 !(OCEAN) + + isfcwnd = 0 !(ATM) + DO n=1,nn if (qc(n) /= 1) CYCLE SELECT CASE(NINT(elm(n))) @@ -547,6 +559,10 @@ SUBROUTINE monit_dep(nn,elm,dep,qc) rmse_sss = rmse_sss + dep(n)**2 !(OCEAN) bias_sss = bias_sss + dep(n) !(OCEAN) isss = isss + 1 !(OCEAN) + CASE(id_sfcwnd_obs) ! (ATM) + rmse_sfcwnd = rmse_sfcwnd + dep(n)**2 ! (ATM) + bias_sfcwnd = bias_sfcwnd + dep(n) ! (ATM) + isfcwnd = isfcwnd + 1 ! (ATM) END SELECT enddo if (iu == 0) then @@ -605,14 +621,21 @@ SUBROUTINE monit_dep(nn,elm,dep,qc) rmse_sss = SQRT(rmse_sss / REAL(isss,r_size)) !(OCEAN) bias_sss = bias_sss / REAL(isss,r_size) !(OCEAN) endif + if (isfcwnd == 0) then !(ATM) + rmse_sfcwnd = undef !(ATM) + bias_sfcwnd = undef !(ATM) + else !(ATM) + rmse_sfcwnd = SQRT(rmse_sfcwnd / REAL(isfcwnd,r_size)) !(ATM) + bias_sfcwnd = bias_sfcwnd / REAL(isfcwnd,r_size) !(ATM) + endif !(ATM) WRITE(6,'(A)') '== OBSERVATIONAL DEPARTURE =============================================' - WRITE(6,'(7A12)') 'U','V','T','S','SSH','eta','SST','SSS' !(OCEAN) - WRITE(6,'(7ES12.3)') bias_u,bias_v,bias_t,bias_s,bias_ssh,bias_eta,bias_sst,bias_sss !(OCEAN) - WRITE(6,'(7ES12.3)') rmse_u,rmse_v,rmse_t,rmse_s,rmse_ssh,rmse_eta,rmse_sst,rmse_sss !(OCEAN) + WRITE(6,'(9A12)') 'U','V','T','S','SSH','eta','SST','SSS','SFCWND' !(OCEAN) + WRITE(6,'(9ES12.3)') bias_u,bias_v,bias_t,bias_s,bias_ssh,bias_eta,bias_sst,bias_sss,bias_sfcwnd !(OCEAN) + WRITE(6,'(9ES12.3)') rmse_u,rmse_v,rmse_t,rmse_s,rmse_ssh,rmse_eta,rmse_sst,rmse_sss,rmse_sfcwnd !(OCEAN) WRITE(6,'(A)') '== NUMBER OF OBSERVATIONS TO BE ASSIMILATED ============================' - WRITE(6,'(7A12)') 'U','V','T','S','SSH','eta','SST','SSS' !(OCEAN) - WRITE(6,'(7I12)') iu,iv,it,is,issh,ieta,isst,isss !(OCEAN) + WRITE(6,'(9A12)') 'U','V','T','S','SSH','eta','SST','SSS','SFCWND' !(OCEAN) + WRITE(6,'(9I12)') iu,iv,it,is,issh,ieta,isst,isss,isfcwnd !(OCEAN) WRITE(6,'(A)') '========================================================================' END SUBROUTINE monit_dep @@ -630,6 +653,7 @@ SUBROUTINE get_nobs(cfile,nrec,nn,errIfNoObs) INTEGER :: ios INTEGER :: iu,iv,it,is,issh,ieta,isst,isss,ix,iy,iz !(OCEAN) INTEGER :: nprof !(OCEAN) + INTEGER :: isfcwnd !(ATM) REAL(r_sngl) :: lon_m1, lat_m1 INTEGER :: iunit LOGICAL :: ex, errIfNoObs_ @@ -654,6 +678,8 @@ SUBROUTINE get_nobs(cfile,nrec,nn,errIfNoObs) nprof = 0 !(OCEAN) lon_m1 = 0 !(OCEAN) lat_m1 = 0 !(OCEAN) + + isfcwnd = 0 !(ATM) iunit=91 if (dodebug) print *, "get_nobs::" INQUIRE(FILE=cfile,EXIST=ex) @@ -692,6 +718,8 @@ SUBROUTINE get_nobs(cfile,nrec,nn,errIfNoObs) iy = iy + 1 !(OCEAN) CASE(id_z_obs) !(OCEAN) iz = iz + 1 !(OCEAN) + CASE(id_sfcwnd_obs) !(ATM) + isfcwnd = isfcwnd + 1 !(ATM) END SELECT if ((wk(2) .ne. lon_m1) .or. (wk(3) .ne. lat_m1)) then nprof=nprof+1 @@ -713,6 +741,7 @@ SUBROUTINE get_nobs(cfile,nrec,nn,errIfNoObs) WRITE(6,'(A12,I10)') ' X:',ix !(OCEAN) WRITE(6,'(A12,I10)') ' Y:',iy !(OCEAN) WRITE(6,'(A12,I10)') ' Z:',iz !(OCEAN) + WRITE(6,'(A12,I10)') ' SFCWND:',isfcwnd !(ATM) CLOSE(iunit) else WRITE(6,'(2A)') cfile,' does not exist -- skipped' diff --git a/src/obs/decoder_sst_acspo.f90 b/src/obs/decoder_sst_acspo.f90 new file mode 100644 index 0000000..c5605b8 --- /dev/null +++ b/src/obs/decoder_sst_acspo.f90 @@ -0,0 +1,161 @@ +PROGRAM decoder_sst_acspo + USE common + USE params_model + USE vars_model + USE common_oceanmodel + !USE vars_obs + USE common_obs_oceanmodel + USE read_sst_acspo, ONLY: read_acspo_nc, sst_acspo_data + + IMPLICIT NONE + + !----------------------------------------------------------------------------- + ! Command line inputs: + !----------------------------------------------------------------------------- + CHARACTER(slen) :: obsinfile='obsin.nc' !IN (default) + CHARACTER(slen) :: obsoutfile='obsout.dat' !OUT(default) + + !----------------------------------------------------------------------------- + ! Obs data arrays + !----------------------------------------------------------------------------- + TYPE(sst_acspo_data), ALLOCATABLE :: obs_data(:) + INTEGER :: nobs + CHARACTER(10) :: Syyyymmddhh = "YYYYMMDDHH" + ! 1234567890 + INTEGER :: delta_seconds = -999 ! max difference from Syyymmddhh + CHARACTER(256) :: valid_id = "NONE" + + 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 + !----------------------------------------------------------------------------- + if (delta_seconds>0) then + CALL read_acspo_nc(trim(obsinfile), min_quality_level, obs_data, nobs, valid_id, & + Syyyymmddhh, delta_seconds) + else + CALL read_acspo_nc(trim(obsinfile), min_quality_level, obs_data, nobs, valid_id) + end if + + 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_sst.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('-validid') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + !read (arg2,*) min_quality_level + valid_id = trim(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), " )" + stop 1 + end select +enddo + +END SUBROUTINE process_command_line + +END PROGRAM decoder_sst_acspo diff --git a/src/obs/obsop_sst_acspo.f90 b/src/obs/obsop_sst_acspo.f90 new file mode 100644 index 0000000..2f90f4e --- /dev/null +++ b/src/obs/obsop_sst_acspo.f90 @@ -0,0 +1,537 @@ +PROGRAM obsop_sst_acspo + USE common + USE params_model + USE vars_model + USE common_oceanmodel + USE params_obs, ONLY: nobs, id_sst_obs + USE params_obs, ONLY: DO_REMOVE_65N + USE vars_obs + USE common_obs_oceanmodel + USE read_sst_acspo, ONLY: read_acspo_nc, sst_acspo_data +#ifdef DYNAMIC + USE input_nml_oceanmodel, ONLY: read_input_namelist +#endif + + IMPLICIT NONE + + !----------------------------------------------------------------------------- + ! Command line inputs: + !----------------------------------------------------------------------------- + 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) + REAL(r_size) :: obserr_scaling=1.0d0 !STEVE: use this to scale the input observations + REAL(r_size) :: obs_randselect=1.0d0 !STEVE: use this to select random input observations + + !----------------------------------------------------------------------------- + ! Obs data arrays + !----------------------------------------------------------------------------- + LOGICAL :: BINARY_INPUT = .FALSE. + CHARACTER(10) :: Syyyymmddhh = "YYYYMMDDHH" + INTEGER :: delta_seconds = -999 + CHARACTER(256) :: valid_id = "NONE" + TYPE(sst_acspo_data), ALLOCATABLE :: obs_data(:) + + 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(:) + + !----------------------------------------------------------------------------- + ! Model background data arrays + !----------------------------------------------------------------------------- + REAL(r_size), ALLOCATABLE :: v3d(:,:,:,:) + REAL(r_size), ALLOCATABLE :: v2d(:,:,:) + + !----------------------------------------------------------------------------- + ! Miscellaneous + !----------------------------------------------------------------------------- + REAL(r_size) :: dk,tg,qg + REAL(r_size) :: ri,rj,rk + INTEGER :: n, i,j,k + REAL(r_size), DIMENSION(1) :: rand + + LOGICAL :: remap_obs_coords = .true. + + INTEGER :: bdyobs=2 !STEVE: use of boundary obs. + ! 1 := less restrictive, remove obs inside boundary + ! 2 := remove all observations touching a boundary + + !----------------------------------------------------------------------------- + ! Debugging parameters + !----------------------------------------------------------------------------- + LOGICAL :: debug_obsfilter = .false. + !STEVE: to adjust writing to output file + LOGICAL :: verbose = .false. + LOGICAL :: dodebug1 = .false. + LOGICAL :: print1st = .true. + + INTEGER :: cnt_obs_u=0, cnt_obs_v=0, cnt_obs_t=0, cnt_obs_s=0 + INTEGER :: cnt_obs_ssh=0, cnt_obs_sst=0, cnt_obs_sss=0, cnt_obs_eta=0 + INTEGER :: cnt_obs_x=0, cnt_obs_y=0, cnt_obs_z=0 + INTEGER, DIMENSION(nv3d+nv2d), SAVE :: cnt_obs = 0 + + !STEVE: for debugging observation culling: + INTEGER :: cnt_yout=0, cnt_xout=0, cnt_zout=0, cnt_triout=0 + INTEGER :: cnt_rigtnlon=0, cnt_nearland=0, cnt_oerlt0=0, cnt_altlatrng=0 + + !----------------------------------------------------------------------------- + ! Instantiations specific to this observation type: + !----------------------------------------------------------------------------- + INTEGER :: min_quality_level=5 ! CDA + INTEGER :: typ = id_sst_obs + LOGICAL :: DO_SUPEROBS = .false. + REAL(r_size), DIMENSION(:,:), ALLOCATABLE :: superobs, delta, M2 ! for online computation of the mean and variance + INTEGER, DIMENSION(:,:), ALLOCATABLE :: supercnt + INTEGER :: idx + INTEGER :: cnt_obs_thinning = 0 + REAL(r_size) :: min_oerr = 0.2 !(K) + + !----------------------------------------------------------------------------- + ! Initialize the common_oceanmodel module, and process command line options + !----------------------------------------------------------------------------- +#ifdef DYNAMIC + CALL read_input_namelist +#endif + CALL set_common_oceanmodel + CALL process_command_line !(get: -obsin -gues -obsout ) + + + !----------------------------------------------------------------------------- + ! Read observations from file + !----------------------------------------------------------------------------- + call create_empty_obsfile(trim(obsoutfile)) + + if (BINARY_INPUT) then ![processed binary input] + print *, "obsop_sst_acspo.f90:: reading binary file=",trim(obsinfile) + CALL get_nobs(trim(obsinfile),8,nobs,errIfNoObs=.false.) + + if (nobs<=0) then + print*, "obsop_sst_acspo.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] + print *, "obsop_sst_acspo.f90:: reading nc file=",trim(obsinfile) + if (delta_seconds>0) then + CALL read_acspo_nc(trim(obsinfile), min_quality_level, obs_data, nobs, valid_id, & + Syyyymmddhh, delta_seconds) + else + CALL read_acspo_nc(trim(obsinfile), min_quality_level, obs_data, nobs, valid_id) + end if + + if (nobs<0) then + print*, "obsop_sst_acspo.f90: 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) ) + + print *, "obsop_sst_acspo.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.d0 + oqc(i) = 0 + obhr(i) = obs_data(i)%hour + enddo + DEALLOCATE(obs_data) + end if + + 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) + print *, "odat(1:40) = ", odat(1:40) + print *, "oerr(1:40) = ", oerr(1:40) + endif + + !----------------------------------------------------------------------------- + ! Update the coordinate to match the model grid + !----------------------------------------------------------------------------- + if (remap_obs_coords) then + CALL center_obs_coords(rlon,oerr,nobs) + endif + + !----------------------------------------------------------------------------- + ! 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..." + + !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 + + 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) + 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)) +! if (dodebug1) print *, "supercnt(",i,",",j,") = ", supercnt(i,j) + enddo + !"superobs" contains the mean + !M2 contains the variance: + WHERE(supercnt > 1) M2 = M2 / (supercnt - 1) + + idx=0 + do j=1,nlat-1 + do i=1,nlon-1 + 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) !(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_sst_obs + if (dodebug1) print *, "odat(idx) = ", odat(idx) + if (dodebug1) print *, "oerr(idx) = ", oerr(idx) + if (dodebug1) print *, "rlon(idx) = ", rlon(idx) + if (dodebug1) print *, "rlat(idx) = ", rlat(idx) + if (dodebug1) print *, "ocnt(idx) = ", supercnt(i,j) + endif + enddo + enddo + 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 + endif + + !----------------------------------------------------------------------------- + ! Read model forecast for this member + !----------------------------------------------------------------------------- + ALLOCATE( v3d(nlon,nlat,nlev,nv3d) ) + ALLOCATE( v2d(nlon,nlat,nv2d) ) + !STEVE: (ISSUE) this can be improved by only reading in the file for needed for this data. + CALL read_diag(guesfile,v3d,v2d) + WRITE(6,*) '****************' + + !----------------------------------------------------------------------------- + ! Cycle through all observations + !----------------------------------------------------------------------------- + WRITE(6,*) "Cycling through observations..." + + ohx=0.0d0 + oqc=0 + idx=0 + DO n=1,nobs + ! Thin observations + ! Select random subset of observations to speed up processing + !STEVE: I know this would be more efficient up above when obs are first introduced. + thin : if (obs_randselect < 1.0d0) then + CALL com_rand(1,rand) + if (rand(1) > obs_randselect) then + cnt_obs_thinning = cnt_obs_thinning + 1 + CYCLE + endif + endif thin + + !--------------------------------------------------------------------------- + ! Count bad obs errors associated with observations, and skip the ob + !--------------------------------------------------------------------------- + if (oerr(n) <= 0) then + !STEVE: this occurred for a few synthetic obs, perhaps due to Dave's code generating obs errors + cnt_oerlt0 = cnt_oerlt0 + 1 + CYCLE + endif + + !--------------------------------------------------------------------------- + ! Convert the physical coordinate to model grid coordinate (note: real, not integer) + !--------------------------------------------------------------------------- + CALL phys2ijk(elem(n),rlon(n),rlat(n),rlev(n),ri,rj,rk) !(OCEAN) + + !--------------------------------------------------------------------------- + ! Filter in the tripolar region until localization is examined in the arctic !(ISSUE) + !--------------------------------------------------------------------------- + if (DO_REMOVE_65N .and. rlat(n) > 65) then + if (verbose) WRITE(6,'(A)') "Latitude above 65.0N, in tripolar region. Removing observation..." + cnt_triout = cnt_triout + 1 + CYCLE + endif + + !--------------------------------------------------------------------------- + ! Filter out observations that are out of range for the grid + !--------------------------------------------------------------------------- + if (CEILING(ri) < 2 .OR. nlon+1 < CEILING(ri)) then + if (verbose) WRITE(6,'(A)') '* X-coordinate out of range' + if (verbose) WRITE(6,'(A,F6.2,A,F6.2)') '* ri=',ri,', olon=', rlon(n) + cnt_xout = cnt_xout + 1 + CYCLE + endif + if (CEILING(rj) < 2 .OR. nlat < CEILING(rj)) then + if (verbose) WRITE(6,'(A)') '* Y-coordinate out of range' + if (verbose) WRITE(6,'(A,F6.2,A,F6.2)') '* rj=',rj,', olat=',rlat(n) + cnt_yout = cnt_yout + 1 + CYCLE + endif + !STEVE: check against kmt, not nlev (OCEAN) + if (CEILING(rk) > nlev) then + CALL itpl_2d(kmt0,ri,rj,dk) + WRITE(6,'(A)') '* Z-coordinate out of range' + WRITE(6,'(A,F6.2,A,F10.2,A,F6.2,A,F6.2,A,F10.2)') & + & '* rk=',rk,', olev=',rlev(n),& + & ', (lon,lat)=(',rlon(n),',',rlat(n),'), kmt0=',dk + cnt_zout = cnt_zout + 1 + CYCLE + endif + if (CEILING(rk) < 2 .AND. rk < 1.00001d0) then !(OCEAN) + rk = 1.00001d0 !(OCEAN) + endif !(OCEAN) + + !--------------------------------------------------------------------------- + ! Check the observation against boundaries + !--------------------------------------------------------------------------- + !STEVE: Check to make sure it's in the ocean, as determined (OCEAN) + ! by mom4's topography map. + ! (note: must do it after coordinate checks, or the coordinate + ! could be outside of the range of the kmt array) + boundary_points : if (ri > nlon) then + !STEVE: I have to check what it does for this case... + ! but it causes an error in the next line if ri > nlon + if (verbose) WRITE(6,'(A)') '* coordinate is not on mom4 model grid: ri > nlon' + cnt_rigtnlon = cnt_rigtnlon + 1 + if (cnt_rigtnlon <= 1) then + WRITE(6,*) "STEVE: ri > nlon (cnt_rigtnlon)" + WRITE(6,*) "ri = ", ri + WRITE(6,*) "nlon = ", nlon + WRITE(6,*) "rj = ", rj + WRITE(6,*) "elem(n) = ", elem(n) + WRITE(6,*) "rlon(n) = ", rlon(n) + WRITE(6,*) "rlat(n) = ", rlat(n) + WRITE(6,*) "rlev(n) = ", rlev(n) + WRITE(6,*) "rk = ", rk + endif + CYCLE + else + !STEVE: check this, case 1 allows more observations, case 2 is more restrictive + select case (bdyobs) + case(1) + if (kmt(NINT(ri),NINT(rj)) .lt. 1) then + if (debug_obsfilter) then + WRITE(6,'(A)') '* coordinate is on or too close to land, according to kmt' + WRITE(6,'(A,F6.2,A,F6.2)') '* ri=',ri,', rj=',rj + WRITE(6,*) "kmt cell = ", kmt(NINT(ri),NINT(rj)) + endif + cnt_nearland = cnt_nearland + 1 + CYCLE + elseif (kmt(NINT(ri),NINT(rj)) .lt. rk) then + if (debug_obsfilter) then + WRITE(6,'(A)') '* coordinate is on or too close to underwater topography, according to kmt' + WRITE(6,'(A,F6.2,A,F6.2,A,F6.2)') '* ri=',ri,', rj=',rj, ', rk=',rk + WRITE(6,*) "kmt cell = ", kmt(NINT(ri),NINT(rj)) + endif + cnt_nearland = cnt_nearland + 1 + CYCLE + endif + case(2) + if(kmt(CEILING(ri),CEILING(rj)) .lt. 1 .or. & + kmt(CEILING(ri),FLOOR(rj)) .lt. 1 .or. & + kmt(FLOOR(ri),CEILING(rj)) .lt. 1 .or. & + kmt(FLOOR(ri),FLOOR(rj)) .lt. 1) THEN + + if (debug_obsfilter) then + WRITE(6,'(A)') '* coordinate is too close to land, according to kmt' + WRITE(6,'(A,F6.2,A,F6.2)') '* ri=',ri,', rj=',rj + WRITE(6,*) "kmt cell = ", kmt(FLOOR(ri):CEILING(ri),FLOOR(rj):CEILING(rj)) + endif + cnt_nearland = cnt_nearland + 1 + CYCLE + elseif(kmt(CEILING(ri),CEILING(rj)) .lt. rk .or. & + kmt(CEILING(ri),FLOOR(rj)) .lt. rk .or. & + kmt(FLOOR(ri),CEILING(rj)) .lt. rk .or. & + kmt(FLOOR(ri),FLOOR(rj)) .lt. rk) THEN + + if (debug_obsfilter) then + WRITE(6,'(A)') '* coordinate is too close to underwater topography, according to kmt' + WRITE(6,'(A,F6.2,A,F6.2,A,F6.2)') '* ri=',ri,', rj=',rj, ', rk=',rk + WRITE(6,*) "kmt cell = ", kmt(FLOOR(ri):CEILING(ri),FLOOR(rj):CEILING(rj)) + endif + cnt_nearland = cnt_nearland + 1 + CYCLE + endif + end select + endif boundary_points + + !--------------------------------------------------------------------------- + ! observation operator (computes H(x)) for specified member + !--------------------------------------------------------------------------- + CALL Trans_XtoY(elem(n),ri,rj,rk,v3d,v2d,ohx(n)) + idx=idx+1 + oqc(n) = 1 + enddo !1:nobs + + !----------------------------------------------------------------------------- + ! Print out the counts of observations removed for various reasons + !----------------------------------------------------------------------------- + WRITE(6,*) "In obsop_sst_acspo.f90::" + WRITE(6,*) "observations at start = ", nobs + WRITE(6,*) "== observations removed for: ==" + WRITE(6,*) "cnt_obs_thinning = ", cnt_obs_thinning + WRITE(6,*) "cnt_oerlt0 = ", cnt_oerlt0 + WRITE(6,*) "cnt_xout = ", cnt_xout + WRITE(6,*) "cnt_yout = ", cnt_yout + WRITE(6,*) "cnt_zout = ", cnt_zout + WRITE(6,*) "cnt_triout = ", cnt_triout + WRITE(6,*) "cnt_rigtnlon = ", cnt_rigtnlon + WRITE(6,*) "cnt_nearland = ", cnt_nearland + WRITE(6,*) "cnt_altlatrng = ", cnt_altlatrng + WRITE(6,*) "===============================" + WRITE(6,*) "observations kept = ", idx + WRITE(6,*) "===============================" + + !----------------------------------------------------------------------------- + ! Write the observations and their associated innovations to file + !----------------------------------------------------------------------------- + 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 ) + +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_acspo.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('-gues') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + guesfile = arg2 + case('-obsout') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + obsoutfile = arg2 + case('-rm65N') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + read (arg2,*) DO_REMOVE_65N + case('-superob') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + read (arg2,*) DO_SUPEROBS + case('-binary') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + read (arg2,*) BINARY_INPUT + case('-thin') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + read (arg2,*) obs_randselect + case('-scale') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + read (arg2,*) obserr_scaling + case('-minerr') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + read (arg2,*) min_oerr + 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('-validid') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + valid_id = trim(arg2) + case('-maxdt') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + read (arg2,*) delta_seconds + case('-debug') + CALL GET_COMMAND_ARGUMENT(i+1,arg2) + PRINT *, "Argument ", i+1, " = ",TRIM(arg2) + read (arg2,*) dodebug1 + 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 obsop_sst_acspo diff --git a/src/obs/params_obs.f90 b/src/obs/params_obs.f90 index 28e1261..aeea767 100644 --- a/src/obs/params_obs.f90 +++ b/src/obs/params_obs.f90 @@ -19,7 +19,8 @@ MODULE params_obs PUBLIC INTEGER,SAVE :: nobs -INTEGER,PARAMETER :: nid_obs=8 !STEVE: sets the dimension of the obs arrays - must be updated depending on choice of obs used. +INTEGER,PARAMETER :: nid_obs=9 !STEVE: sets the dimension of the obs arrays - must be updated depending on choice of obs used. + !CDA: controls which type of obs can affect anal vars, used by var_local(nv3d+nv2d+nv4d,nid_obs) in vars_letkf.90 INTEGER,PARAMETER :: id_u_obs=2819 INTEGER,PARAMETER :: id_v_obs=2820 INTEGER,PARAMETER :: id_t_obs=3073 @@ -41,6 +42,8 @@ MODULE params_obs INTEGER,PARAMETER :: id_ui_obs=8337 !(SIS) ice drift u INTEGER,PARAMETER :: id_vi_obs=8334 !(SIS) ice drift v +INTEGER,PARAMETER :: id_sfcwnd_obs = 100001 ! (ATM) 10-meter winds + !!------------------------------------------------------------ !! unique ID's for observations in COUPLED SYSTEM !! STEVE: the following will replace what is above: diff --git a/src/obs/read_sst_acspo.f90 b/src/obs/read_sst_acspo.f90 new file mode 100644 index 0000000..86683b3 --- /dev/null +++ b/src/obs/read_sst_acspo.f90 @@ -0,0 +1,330 @@ +MODULE read_sst_acspo + USE common, ONLY: r_sngl, r_dble, r_size, slen + USE params_obs, ONLY: id_sst_obs + IMPLICIT NONE + + PRIVATE + + PUBLIC :: sst_acspo_data + PUBLIC :: read_acspo_nc + PRIVATE :: inspect_obs_data + + TYPE sst_acspo_data + REAL(r_size) :: x_grd(2) ! longitude, latitude + REAL(r_size) :: value ! actual physical value of the parameter measured at this grid point + REAL(r_size) :: oerr ! observation standard error + REAL(r_size) :: hour ! Hour of observation + INTEGER :: qkey ! Quality key + INTEGER :: typ ! type of observation (elem) + LOGICAL :: kept ! tells letkf whether this obs is kept for assimilation + END TYPE sst_acspo_data + +CONTAINS + +SUBROUTINE read_acspo_nc(obsinfile, min_quality_level, obs_data, nobs, validId, Syyyymmddhh, delta_seconds) + USE m_ncio, ONLY: nc_get_fid, nc_close_fid, nc_rddim, nc_rdvar2d, nc_rdvar1d, & + nc_rdatt, nc_fndvar + IMPLICIT NONE + + CHARACTER(*),INTENT(IN) :: obsinfile + INTEGER, INTENT(IN) :: min_quality_level + TYPE(sst_acspo_data),ALLOCATABLE,INTENT(INOUT) :: obs_data(:) + INTEGER, INTENT(OUT) :: nobs + CHARACTER(*), INTENT(IN) :: validId ! validation id from global attr e.g."AVHRRF_MA-STAR-L2P-v2.8". default is "NONE" for no check, otherwise check if file id equals validId. Throw errors if unmatch + CHARACTER(10),INTENT(IN), OPTIONAL :: Syyyymmddhh + !1234567890 + INTEGER, INTENT(IN), OPTIONAL :: delta_seconds ! qc'ed if abs(obstime-Syyyymmddhh)>delta_seconds + + INTEGER :: fid, ni, nj, i, j, n, qlev + + REAL(r_size),ALLOCATABLE :: alon2d(:,:), alat2d(:,:) + REAL(r_size),ALLOCATABLE :: sea_surface_temperature(:,:), stde(:,:), sst_dtime(:,:) + INTEGER,ALLOCATABLE :: sst_time_in_seconds_since19780101(:,:) + INTEGER(4), ALLOCATABLE :: quality_level(:,:), l2p_flags_acc(:,:) + INTEGER(1), ALLOCATABLE :: i1buf2d(:,:) + INTEGER(2), ALLOCATABLE :: i2buf2d(:,:) + LOGICAL, ALLOCATABLE :: valid(:,:) ! .true. if no missing info at this grid by checking different FillValues in nc file + INTEGER :: FillValue ! missing value indicated by vars in nc file + REAL(r_size) :: offset, scale_factor + LOGICAL :: dodebug = .true. + INTEGER :: obsdate(8), sdate(8), sdate_in_seconds_since19780101, sdate_in_min_since19780101 + INTEGER :: qcdate(5), qcdate_in_seconds_since19780101, qcdate_in_min_since19780101 + REAL :: tdelta(5) + INTEGER :: time(1), sst_time + CHARACTER(256) :: acspoId + +!------------------------------------------------------------------------------- +! Read dimension & geo info: +!------------------------------------------------------------------------------- + CALL nc_get_fid(trim(obsinfile),fid) + + CALL nc_rdatt(fid, "","id", acspoId) + WRITE(6,*) "[msg] read_acspo_nc::id from file =", trim(acspoId) + WRITE(6,*) "[msg] read_acspo_nc::id from user =", trim(validId) + if (trim(validId) .ne. "NONE") then + if (trim(validId) .ne. trim(acspoId)) then + WRITE(6,*) "[err] validId and file Id does not match." + STOP (25) + endif + endif + + CALL nc_rddim(fid,"ni", ni) + CALL nc_rddim(fid,"nj", nj) + WRITE(6,*) "[msg] read_acspo_nc::ni, nj =", ni, nj + + ALLOCATE(alon2d(ni,nj),alat2d(ni,nj)) + ALLOCATE(valid(ni,nj)) + valid = .true. + + CALL nc_rdvar2d(fid,"lon",alon2d) + WRITE(6,*) "[msg] read_acspo_nc::alon2d: min, max=", & + minval(alon2d),maxval(alon2d) + CALL nc_rdvar2d(fid,"lat",alat2d) + WRITE(6,*) "[msg] read_acspo_nc::alat2d: min, max=", & + minval(alat2d),maxval(alat2d) + +!------------------------------------------------------------------------------- +! Read the observed SST data & other info from obs file (obsinfile) +! [NOTE]: +! The SST data from file uses K as its units, while MOM6 uses degC +! We convert SST from K to degC here. +!------------------------------------------------------------------------------- + + ALLOCATE(i1buf2d(ni,nj)) + ALLOCATE(i2buf2d(ni,nj)) + ALLOCATE(l2p_flags_acc(ni,nj)) + + + CALL nc_rdvar2d(fid, "sea_surface_temperature", i2buf2d) + CALL nc_rdatt(fid, "sea_surface_temperature","scale_factor", scale_factor) + CALL nc_rdatt(fid, "sea_surface_temperature","add_offset", offset) + CALL nc_rdatt(fid, "sea_surface_temperature","_FillValue", FillValue) + if (dodebug) WRITE(6,*) "_FillValue, scale_factor, add_offset=", & + FillValue, scale_factor, offset + where (i2buf2d == FillValue) + valid = .false. + end where + ALLOCATE(sea_surface_temperature(ni,nj)) + sea_surface_temperature = REAL(i2buf2d,r_size)*scale_factor + offset + WRITE(6,*) "[msg] read_acspo_nc::sst: min, max=", & + minval(sea_surface_temperature), maxval(sea_surface_temperature), & + minval(sea_surface_temperature, mask=valid),& + maxval(sea_surface_temperature, mask=valid) + + where (valid) + sea_surface_temperature = cvt_temp_K2C(sea_surface_temperature) + end where + WRITE(6,*) "[msg] read_acspo_nc::sst: Converting from K to degC" + WRITE(6,*) "[msg] read_acspo_nc::sst: min, max=", & + minval(sea_surface_temperature), maxval(sea_surface_temperature), & + minval(sea_surface_temperature, mask=valid),& + maxval(sea_surface_temperature, mask=valid) + +! SST error + CALL nc_rdvar2d(fid, "sses_standard_deviation", i1buf2d) + CALL nc_rdatt(fid, "sses_standard_deviation", "scale_factor", scale_factor) + CALL nc_rdatt(fid, "sses_standard_deviation", "add_offset", offset) + CALL nc_rdatt(fid, "sses_standard_deviation", "_FillValue", FillValue) + if (dodebug) WRITE(6,*) "_FillValue, scale_factor, add_offset=", & + FillValue, scale_factor, offset + where (i1buf2d == FillValue) + valid = .false. + end where + ALLOCATE(stde(ni,nj)) + stde = REAL(i1buf2d,r_size)*scale_factor + offset + WRITE(6,*) "[msg] read_acspo_nc::stde: min, max=", & + minval(stde), maxval(stde), & + minval(stde, mask=valid),maxval(stde, mask=valid) + +!------------------------------------------------------------------------------- +! Read the base time and time offset +!------------------------------------------------------------------------------- + CALL nc_rdvar1d(fid, "time", time) + WRITE(6,*) "[msg] read_acspo_nc::time=",time(1) + sst_time = time(1) + sdate = [1981, 1, 1, 0, 0, 0, 0, 0] ! YYYY/MON/DAY/TZONE/HR/MIN/SEC/MSEC + tdelta = [0.0,0.0,0.0,REAL(time(1)),0.0] ! DAY/HR/MIN/SEC/MSEC + CALL w3movdat(tdelta, sdate, obsdate) + WRITE(6,*) "[msg] read_acspo_nc::obstime=", obsdate(1:3),obsdate(5:7) + call W3FS21(sdate(1:5), sdate_in_min_since19780101) + sdate_in_seconds_since19780101 = 60 * sdate_in_min_since19780101 + WRITE(6,*) "[msg] read_acspo_nc::obstime_since_19780101=", sdate_in_seconds_since19780101 + +! sst_dtime + CALL nc_rdvar2d(fid, "sst_dtime", i2buf2d) + CALL nc_rdatt(fid, "sst_dtime", "scale_factor", scale_factor) + CALL nc_rdatt(fid, "sst_dtime", "add_offset", offset) + CALL nc_rdatt(fid, "sst_dtime", "_FillValue", FillValue) + if (dodebug) WRITE(6,*) "_FillValue, scale_factor, add_offset=", & + FillValue, scale_factor, offset + where (i2buf2d == FillValue) + valid = .false. + end where + ALLOCATE(sst_dtime(ni,nj)) + sst_dtime = REAL(i2buf2d,r_size)*scale_factor + offset + WRITE(6,*) "[msg] read_acspo_nc::sst_dtime: min, max=", & + minval(sst_dtime), maxval(sst_dtime), & + minval(sst_dtime, mask=valid),maxval(sst_dtime, mask=valid) + + ALLOCATE(sst_time_in_seconds_since19780101(ni,nj)) + sst_time_in_seconds_since19780101 = NINT(sst_dtime) + sst_time + & + sdate_in_seconds_since19780101 + WRITE(6,*) "[msg] read_acspo_nc::sst_dtime_since19780101: min, max=", & + minval(sst_time_in_seconds_since19780101, mask=valid), & + maxval(sst_time_in_seconds_since19780101, mask=valid) + + 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_acspo_nc::Syyyymmddhh=", Syyyymmddhh, & + "qcdate(:)=",qcdate(:) + call W3FS21(qcdate,qcdate_in_min_since19780101) + qcdate_in_seconds_since19780101 = 60 * qcdate_in_min_since19780101 + WRITE(6,*) "[msg] read_acspo_nc::qcdate_in_seconds_since19780101, delta_seconds=", & + qcdate_in_seconds_since19780101, delta_seconds + + ! calculate time difference + where (ABS(sst_time_in_seconds_since19780101 - & + qcdate_in_seconds_since19780101) > delta_seconds) + valid = .false. + end where + + end if + +!------------------------------------------------------------------------------- +! Read the l2p_flags (ASCPO writes additional QC flags into it) +!L2P common flags in bits 1-6 and data provider flags (from ACSPO mask) in bits 9-16: +!0: bit01 (0=IR: 1=microwave); +!1: Y bit02 (0=ocean; 1=land); +!2: Y bit03 (0=no ice; 1=ice); +!3-7: bits04-08 (reserved,set to 0); +!8: Y bit09 (0=radiance valid; 1=invalid); +!9: bit10 (0=night; 1=day); +!10: Y bit11 (0=ocean; 1=land); +!11: Y bit12 (0=good quality data; 1=degraded quality data due to \"twilight\" region); +!12: Y bit13 (0=no glint; 1=glint); +!13: Y bit14 (0=no snow/ice; 1=snow/ice); +!14: Y bits15-16 (00=clear; 01=probably clear; 10=cloudy; 11=clear-sky mask undefined)" ; +!------------------------------------------------------------------------------- + CALL nc_rdvar2d(fid, "l2p_flags", i2buf2d) + l2p_flags_acc = 0 + l2p_flags_acc = l2p_flags_acc + ibits(i2buf2d,1,1) + & ! skip land + ibits(i2buf2d,2,1) + & ! skip ice + ibits(i2buf2d,8,1) + & ! skip invalid rad + ibits(i2buf2d,10,1) + & ! skip land + ibits(i2buf2d,11,1) + & ! skip twlight + ibits(i2buf2d,12,1) + & ! skip glint + ibits(i2buf2d,13,1) + & ! skip snow/ice + abs(ibits(i2buf2d,14,2)) ! skip non-clear + where (l2p_flags_acc /= 0) + valid = .false. + end where + +!------------------------------------------------------------------------------- +! Read the Quality Key (5 is best quality; 0 missing data) +! 5 == best_quality +! 4-2 == not used +! 1 == bad data +! 0 == missing data +!------------------------------------------------------------------------------- + ALLOCATE(quality_level(ni,nj)) + CALL nc_rdvar2d(fid, "quality_level", quality_level) + CALL nc_rdatt(fid, "quality_level", "_FillValue", FillValue) + if (dodebug) then + WRITE(6,*) "_FillValue=", FillValue + ! check num of obs for each quality level + do i = -1, 5 + i1buf2d = 0 + if (i==-1) then + qlev = FillValue + else + qlev = i + end if + where (quality_level == qlev) + i1buf2d = 1 + end where + n = sum(int(i1buf2d,4)) + WRITE(6,*) "[msg] read_acspo_nc::quality_level, # of obs, %=", i, n, n*100.0/(ni*nj) + end do + end if + where (quality_level == FillValue) + valid = .false. + end where + !WRITE(6,*) "[msg] read_acspo_nc::quality_level: min, max=", & + ! minval(quality_level, mask=valid), maxval(quality_level, mask=valid) + +!------------------------------------------------------------------------------- +! Close the netcdf file +!------------------------------------------------------------------------------- + CALL nc_close_fid(fid) + +!------------------------------------------------------------------------------- +! CONVERT netcdf data to argo_data format: +!------------------------------------------------------------------------------- + WRITE(6,*) "Finished reading netCDF file, formatting data..." + +! determine num of valid obs + WRITE(6,*) "[msg] read_acspo_nc::use obs with min_quality_level=", min_quality_level + nobs = 0 + do j = 1, nj; do i = 1, ni + if (quality_level(i,j) >= min_quality_level.and.valid(i,j)) then + nobs = nobs + 1 + end if + end do; end do + WRITE(6,*) "[msg] read_acspo_nc::nobs_retained, %=", nobs, nobs*100.0/(ni*nj) + +! fill into data struct + ALLOCATE(obs_data(nobs)) + n = 0 + do j = 1, nj + do i = 1, ni + if (quality_level(i,j) >= min_quality_level .and. valid(i,j)) then + n = n + 1 + obs_data(n)%typ = id_sst_obs + obs_data(n)%x_grd(1) = alon2d(i,j) + obs_data(n)%x_grd(2) = alat2d(i,j) + obs_data(n)%hour = sst_time_in_seconds_since19780101(i,j)/3600. + obs_data(n)%value = sea_surface_temperature(i,j) + obs_data(n)%oerr = stde(i,j) + obs_data(n)%qkey = quality_level(i,j) + end if + end do + end do + CALL inspect_obs_data(obs_data) + if (n/=nobs) then + WRITE(6,*) "[err] read_acspo_nc::n/=nobs: n, nobs=", n, nobs + STOP (26) + end if + +END SUBROUTINE read_acspo_nc + +SUBROUTINE inspect_obs_data(obs_data) + IMPLICIT NONE + TYPE(sst_acspo_data),INTENT(IN) :: obs_data(:) + + WRITE(6,*) "[msg] read_acspo_nc::info" + WRITE(6,*) " nobs=", size(obs_data) + WRITE(6,*) " x_grd(1): min, max=", minval(obs_data(:)%x_grd(1)), maxval(obs_data(:)%x_grd(1)) + WRITE(6,*) " x_grd(2): min, max=", minval(obs_data(:)%x_grd(2)), maxval(obs_data(:)%x_grd(2)) + WRITE(6,*) " hour: min, max=", minval(obs_data(:)%hour), maxval(obs_data(:)%hour) + WRITE(6,*) " value: min, max=", minval(obs_data(:)%value), maxval(obs_data(:)%value) + WRITE(6,*) " oerr: min, max=", minval(obs_data(:)%oerr), maxval(obs_data(:)%oerr) + WRITE(6,*) " qkey: min, max=", minval(obs_data(:)%qkey), maxval(obs_data(:)%qkey) + +END SUBROUTINE + + +ELEMENTAL FUNCTION cvt_temp_K2C(tf) RESULT (tc) + IMPLICIT NONE + + REAL(r_size),INTENT(IN) :: tf ! temperature in Kelvin + REAL(r_size) :: tc ! temperature in degree Celsius + + REAL(r_size),PARAMETER :: ADD_OFFSET_K2C = -273.15_r_size + + tc = tf + ADD_OFFSET_K2C + +END FUNCTION + +END MODULE read_sst_acspo diff --git a/src/obs/superob.f90 b/src/obs/superob.f90 index 6d75835..a21d8b0 100644 --- a/src/obs/superob.f90 +++ b/src/obs/superob.f90 @@ -81,9 +81,9 @@ PROGRAM superob ! Read observations from file !----------------------------------------------------------------------------- nobs = 0 - CALL get_nobs(trim(obsinfile),8,nobs) + CALL get_nobs(trim(obsinfile),8,nobs, errIfNoObs=.false.) - if (nobs<0) then ! [NOBS <= 0] + if (nobs<=0) then ! [NOBS <= 0] print*, "[warning] superob :: nobs = 0. will output an empty file" call create_empty_obsfile(trim(obsoutfile))