Skip to content

Commit

Permalink
Add additional obs (#196)
Browse files Browse the repository at this point in the history
* add decoder for AVHRR ACSPO
* fix a bug in superob when nobs=0
  • Loading branch information
gmao-cda authored Mar 22, 2024
1 parent 24d4b10 commit 96729ca
Show file tree
Hide file tree
Showing 9 changed files with 1,081 additions and 11 deletions.
3 changes: 3 additions & 0 deletions src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -91,6 +92,7 @@ set(operators
sst
sst_geostationary
sst_viirs
sst_acspo
sss
tprof
sprof
Expand All @@ -107,6 +109,7 @@ endforeach()
# different observation decoders
set(decoders
sst_viirs
sst_acspo
sss_smap
)
foreach(decoder ${decoders})
Expand Down
8 changes: 6 additions & 2 deletions src/letkf/letkf_obs.f90.new
Original file line number Diff line number Diff line change
Expand Up @@ -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'

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
5 changes: 4 additions & 1 deletion src/letkf/vars_letkf.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
!---------------------------------------------------------------------------
Expand All @@ -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..."
Expand Down
39 changes: 34 additions & 5 deletions src/model_specific/mom6/common_obs_mom6.f90.kdtree
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)))
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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_
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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'
Expand Down
161 changes: 161 additions & 0 deletions src/obs/decoder_sst_acspo.f90
Original file line number Diff line number Diff line change
@@ -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 <obsinfile> -gues <guesfile> -obsout <obsoutfile>)


!-----------------------------------------------------------------------------
! 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
Loading

0 comments on commit 96729ca

Please sign in to comment.