Skip to content

Commit

Permalink
add SCDA for surface wind speeds
Browse files Browse the repository at this point in the history
  • Loading branch information
gmao-cda committed Feb 14, 2024
1 parent 24d4b10 commit fa66205
Show file tree
Hide file tree
Showing 4 changed files with 48 additions and 9 deletions.
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
5 changes: 4 additions & 1 deletion src/obs/params_obs.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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:
Expand Down

0 comments on commit fa66205

Please sign in to comment.