Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add additional obs #196

Merged
merged 3 commits into from
Mar 22, 2024
Merged

Add additional obs #196

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
fa66205
add SCDA for surface wind speeds
gmao-cda Feb 14, 2024
3dc1a8b
add decoder for AVHRR ACSPO
gmao-cda Mar 6, 2024
1b089a2
fix a bug in superob when nobs=0
gmao-cda Mar 8, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions src/CMakeLists.txt
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
Loading