diff --git a/src/enkf/enkf.f90 b/src/enkf/enkf.f90 index d35613b585..479f60c019 100644 --- a/src/enkf/enkf.f90 +++ b/src/enkf/enkf.f90 @@ -51,7 +51,7 @@ module enkf ! NH, tropics and SH, and in the horizontal, vertical and time dimensions, ! using the namelist parameters corrlengthnh, corrlengthtr, corrlengthsh, ! lnsigcutoffnh, lnsigcutofftr, lnsigcutoffsh (lnsigcutoffsatnh, -! lnsigcutoffsattr, lnsigcutoffsatsh for satellite obs, similar for ps obs) +! lnsigcutoffsattr, lnsigcutoffsatsh for satellite obs, similar for ps and fed obs) ! obtimelnh, obtimeltr, obtimelsh. The length scales should be given in km for the ! horizontal, hours for time, and 'scale heights' (units of -log(p/pref)) in the ! vertical. The function used for localization (function taper) diff --git a/src/enkf/enkf_obs_sensitivity.f90 b/src/enkf/enkf_obs_sensitivity.f90 index 6c37936f31..72296d5934 100644 --- a/src/enkf/enkf_obs_sensitivity.f90 +++ b/src/enkf/enkf_obs_sensitivity.f90 @@ -36,6 +36,7 @@ module enkf_obs_sensitivity use params, only: efsoi_flag,latbound,nlevs,nanals,datestring, & lnsigcutoffsatnh,lnsigcutoffsattr,lnsigcutoffsatsh, & lnsigcutoffpsnh,lnsigcutoffpstr,lnsigcutoffpssh, & + lnsigcutofffednh,lnsigcutofffedtr,lnsigcutofffedsh, & lnsigcutoffnh,lnsigcutofftr,lnsigcutoffsh, & corrlengthnh,corrlengthtr,corrlengthsh, & obtimelnh,obtimeltr,obtimelsh,letkf_flag, & @@ -292,6 +293,8 @@ subroutine read_ob_sens lnsigl(nob) = latval(deglat,lnsigcutoffsatnh,lnsigcutoffsattr,lnsigcutoffsatsh) else if (obtype(nob)(1:3) == ' ps') then lnsigl(nob) = latval(deglat,lnsigcutoffpsnh,lnsigcutoffpstr,lnsigcutoffpssh) + else if (obtype(nob)(1:3) == 'fed') then + lnsigl(nob) = latval(deglat,lnsigcutofffednh,lnsigcutofffedtr,lnsigcutofffedsh) else lnsigl(nob)=latval(deglat,lnsigcutoffnh,lnsigcutofftr,lnsigcutoffsh) end if diff --git a/src/enkf/enkf_obsmod.f90 b/src/enkf/enkf_obsmod.f90 index ea8f6446fb..eb4f9c8e58 100644 --- a/src/enkf/enkf_obsmod.f90 +++ b/src/enkf/enkf_obsmod.f90 @@ -109,6 +109,8 @@ module enkf_obsmod lnsigcutoffnh, lnsigcutoffsh, lnsigcutofftr, corrlengthnh,& corrlengthtr, corrlengthsh, obtimelnh, obtimeltr, obtimelsh,& lnsigcutoffsatnh, lnsigcutoffsatsh, lnsigcutoffsattr,& + lnsigcutofffednh, lnsigcutofffedsh, lnsigcutofffedtr,& + corrlengthfednh, corrlengthfedtr, corrlengthfedsh, & varqc, huber, zhuberleft, zhuberright, modelspace_vloc, & lnsigcutoffpsnh, lnsigcutoffpssh, lnsigcutoffpstr, neigv, & lnsigcutoffrdrnh, lnsigcutoffrdrsh, lnsigcutoffrdrtr,& @@ -276,6 +278,8 @@ subroutine readobs() lnsigl(nob) = latval(deglat,lnsigcutoffsatnh,lnsigcutoffsattr,lnsigcutoffsatsh) else if (obtype(nob)(1:3) == ' ps') then lnsigl(nob) = latval(deglat,lnsigcutoffpsnh,lnsigcutoffpstr,lnsigcutoffpssh) + else if (obtype(nob)(1:3) == 'fed') then + lnsigl(nob) = latval(deglat,lnsigcutofffednh,lnsigcutofffedtr,lnsigcutofffedsh) else if ( (obtype(nob)(1:3) == 'dbz' .or. obtype(nob)(1:3) == ' rw') .and. l_use_enkf_directZDA ) then lnsigl(nob) = latval(deglat,lnsigcutoffrdrnh,lnsigcutoffrdrtr,lnsigcutoffrdrsh) else @@ -293,6 +297,9 @@ subroutine readobs() if ( (obtype(nob)(1:3) == 'dbz' .or. obtype(nob)(1:3) == ' rw') .and. l_use_enkf_directZDA ) then corrlengthsq(nob)=latval(deglat,corrlengthrdrnh,corrlengthrdrtr,corrlengthrdrsh)**2 end if + if (obtype(nob)(1:3) == 'fed') then + corrlengthsq(nob)=latval(deglat,corrlengthfednh,corrlengthfedtr,corrlengthfedsh)**2 + end if obtimel(nob)=latval(deglat,obtimelnh,obtimeltr,obtimelsh) end do diff --git a/src/enkf/innovstats.f90 b/src/enkf/innovstats.f90 index 68668218fc..853532c9b9 100644 --- a/src/enkf/innovstats.f90 +++ b/src/enkf/innovstats.f90 @@ -45,6 +45,7 @@ subroutine print_innovstats(obfit,obsprd) nobsspd_nh,nobsspd_sh,nobsspd_tr,& nobsgps_nh,nobsgps_sh,nobsgps_tr,& nobsdbz_nh,nobsdbz_sh,nobsdbz_tr,& + nobsfed_nh,nobsfed_sh,nobsfed_tr,& nobsrw_nh,nobsrw_sh,nobsrw_tr,& nobsq_nh,nobsq_sh,nobsq_tr,nobswnd_nh,nobswnd_sh,nobswnd_tr,& nobsoz_nh,nobsoz_sh,nobsoz_tr,nobsps_sh,nobsps_nh,nobsps_tr,nob @@ -67,6 +68,9 @@ subroutine print_innovstats(obfit,obsprd) sumdbz_nh,biasdbz_nh,sumdbz_spread_nh,sumdbz_oberr_nh,& sumdbz_sh,biasdbz_sh,sumdbz_spread_sh,sumdbz_oberr_sh,& sumdbz_tr,biasdbz_tr,sumdbz_spread_tr,sumdbz_oberr_tr,& + sumfed_nh,biasfed_nh,sumfed_spread_nh,sumfed_oberr_nh,& + sumfed_sh,biasfed_sh,sumfed_spread_sh,sumfed_oberr_sh,& + sumfed_tr,biasfed_tr,sumfed_spread_tr,sumfed_oberr_tr,& sumrw_nh,biasrw_nh,sumrw_spread_nh,sumrw_oberr_nh,& sumrw_sh,biasrw_sh,sumrw_spread_sh,sumrw_oberr_sh,& sumrw_tr,biasrw_tr,sumrw_spread_tr,sumrw_oberr_tr,& @@ -112,6 +116,9 @@ subroutine print_innovstats(obfit,obsprd) nobsdbz_nh = 0 nobsdbz_sh = 0 nobsdbz_tr = 0 + nobsfed_nh = 0 + nobsfed_sh = 0 + nobsfed_tr = 0 nobsrw_nh = 0 nobsrw_sh = 0 nobsrw_tr = 0 @@ -168,6 +175,12 @@ subroutine print_innovstats(obfit,obsprd) sumdbz_nh,biasdbz_nh,sumdbz_spread_nh,sumdbz_oberr_nh,nobsdbz_nh,& sumdbz_sh,biasdbz_sh,sumdbz_spread_sh,sumdbz_oberr_sh,nobsdbz_sh,& sumdbz_tr,biasdbz_tr,sumdbz_spread_tr,sumdbz_oberr_tr,nobsdbz_tr) + else if (obtype(nob)(1:3) == 'fed') then + call obstats(obfit(nob),oberrvar_orig(nob),& + obsprd(nob),obloclat(nob),& + sumfed_nh,biasfed_nh,sumfed_spread_nh,sumfed_oberr_nh,nobsfed_nh,& + sumfed_sh,biasfed_sh,sumfed_spread_sh,sumfed_oberr_sh,nobsfed_sh,& + sumfed_tr,biasfed_tr,sumfed_spread_tr,sumfed_oberr_tr,nobsfed_tr) else if (obtype(nob)(1:3) == ' rw') then call obstats(obfit(nob),oberrvar_orig(nob),& obsprd(nob),obloclat(nob),& @@ -216,6 +229,9 @@ subroutine print_innovstats(obfit,obsprd) call printstats(' all dbz',sumdbz_nh,biasdbz_nh,sumdbz_spread_nh,sumdbz_oberr_nh,nobsdbz_nh,& sumdbz_sh,biasdbz_sh,sumdbz_spread_sh,sumdbz_oberr_sh,nobsdbz_sh,& sumdbz_tr,biasdbz_tr,sumdbz_spread_tr,sumdbz_oberr_tr,nobsdbz_tr) + call printstats(' all fed',sumfed_nh,biasfed_nh,sumfed_spread_nh,sumfed_oberr_nh,nobsfed_nh,& + sumfed_sh,biasfed_sh,sumfed_spread_sh,sumfed_oberr_sh,nobsfed_sh,& + sumfed_tr,biasfed_tr,sumfed_spread_tr,sumfed_oberr_tr,nobsfed_tr) call printstats(' all rw',sumrw_nh,biasq_nh,sumrw_spread_nh,sumrw_oberr_nh,nobsrw_nh,& sumrw_sh,biasrw_sh,sumrw_spread_sh,sumrw_oberr_sh,nobsrw_sh,& sumrw_tr,biasrw_tr,sumrw_spread_tr,sumrw_oberr_tr,nobsrw_tr) diff --git a/src/enkf/params.f90 b/src/enkf/params.f90 index 62701c24a7..36b0c9c207 100644 --- a/src/enkf/params.f90 +++ b/src/enkf/params.f90 @@ -124,6 +124,8 @@ module params real(r_single),public :: lnsigcutoffnh,lnsigcutofftr,lnsigcutoffsh,& lnsigcutoffsatnh,lnsigcutoffsattr,lnsigcutoffsatsh,& lnsigcutoffpsnh,lnsigcutoffpstr,lnsigcutoffpssh +real(r_single),public :: corrlengthfednh,corrlengthfedtr,corrlengthfedsh, & + lnsigcutofffednh,lnsigcutofffedtr,lnsigcutofffedsh real(r_single),public :: corrlengthrdrnh,corrlengthrdrtr,corrlengthrdrsh, & lnsigcutoffrdrnh,lnsigcutoffrdrtr,lnsigcutoffrdrsh real(r_single),public :: analpertwtnh,analpertwtsh,analpertwttr,sprd_tol,saterrfact @@ -261,6 +263,8 @@ module params lnsigcutoffnh,lnsigcutofftr,lnsigcutoffsh,& lnsigcutoffsatnh,lnsigcutoffsattr,lnsigcutoffsatsh,& lnsigcutoffpsnh,lnsigcutoffpstr,lnsigcutoffpssh,& + corrlengthfednh,corrlengthfedsh,corrlengthfedtr,& + lnsigcutofffednh,lnsigcutofffedsh,lnsigcutofffedtr,& fgfileprefixes,fgsfcfileprefixes,anlfileprefixes, & anlsfcfileprefixes,incfileprefixes,incsfcfileprefixes,& statefileprefixes,statesfcfileprefixes, & @@ -317,6 +321,10 @@ subroutine read_namelist() corrlengthrdrnh = 10 corrlengthrdrtr = 10 corrlengthrdrsh = 10 +! corrlength (km) for GLM flash extent density +corrlengthfednh = 30_r_single +corrlengthfedtr = 30_r_single +corrlengthfedsh = 30_r_single ! read in localization length scales from an external file. readin_localization = .false. ! min and max inflation. @@ -341,6 +349,9 @@ subroutine read_namelist() lnsigcutoffrdrnh = 0.2_r_single ! value for radar lnsigcutoffrdrtr = 0.2_r_single ! value for radar lnsigcutoffrdrsh = 0.2_r_single ! value for radar +lnsigcutofffednh = 2._r_single ! value for GLM flash extent density +lnsigcutofffedtr = 2._r_single ! value for GLM flash extent density +lnsigcutofffedsh = 2._r_single ! value for GLM flash extent density ! ob time localization obtimelnh = 1.e10_r_single obtimeltr = 1.e10_r_single @@ -813,6 +824,10 @@ subroutine read_namelist() corrlengthrdrnh = corrlengthrdrnh * 1.e3_r_single/rearth corrlengthrdrtr = corrlengthrdrtr * 1.e3_r_single/rearth corrlengthrdrsh = corrlengthrdrsh * 1.e3_r_single/rearth +! rescale covariance localization length for GLM FED +corrlengthfednh = corrlengthfednh * 1.e3_r_single/rearth +corrlengthfedtr = corrlengthfedtr * 1.e3_r_single/rearth +corrlengthfedsh = corrlengthfedsh * 1.e3_r_single/rearth ! convert targe area boundary into radians tar_minlat = tar_minlat * deg2rad diff --git a/src/enkf/readconvobs.f90 b/src/enkf/readconvobs.f90 index d1f4ec3ff8..a5383069a1 100644 --- a/src/enkf/readconvobs.f90 +++ b/src/enkf/readconvobs.f90 @@ -42,9 +42,9 @@ module readconvobs !> observation types to read from netcdf files -integer(i_kind), parameter :: nobtype = 11 +integer(i_kind), parameter :: nobtype = 12 character(len=3), dimension(nobtype), parameter :: obtypes = (/' t', ' q', ' ps', ' uv', 'tcp', & - 'gps', 'spd', ' pw', ' dw', ' rw', 'dbz' /) + 'gps', 'spd', ' pw', ' dw', ' rw', 'dbz', 'fed' /) contains @@ -79,7 +79,7 @@ subroutine get_num_convobs_bin(obspath,datestring,num_obs_tot,num_obs_totdiag,id integer(i_kind) :: iunit, nchar, nreal, ii, mype, ios, idate, i, ipe, ioff0 integer(i_kind),dimension(2) :: nn,nobst, nobsps, nobsq, nobsuv, nobsgps, & nobstcp,nobstcx,nobstcy,nobstcz,nobssst, nobsspd, nobsdw, nobsrw, nobspw, & - nobsdbz + nobsdbz, nobsfed character(8),allocatable,dimension(:):: cdiagbuf real(r_single),allocatable,dimension(:,:)::rdiagbuf real(r_kind) :: errorlimit,errorlimit2,error,pres,obmax @@ -104,6 +104,7 @@ subroutine get_num_convobs_bin(obspath,datestring,num_obs_tot,num_obs_totdiag,id nobspw = 0 nobsgps = 0 nobsdbz = 0 + nobsfed = 0 nobstcp = 0; nobstcx = 0; nobstcy = 0; nobstcz = 0 init_pass = .true. peloop: do ipe=0,npefiles @@ -187,6 +188,9 @@ subroutine get_num_convobs_bin(obspath,datestring,num_obs_tot,num_obs_totdiag,id else if (obtype == 'dbz') then nobsdbz = nobsdbz + nn num_obs_tot = num_obs_tot + nn(2) + else if (obtype == 'fed') then + nobsfed = nobsfed + nn + num_obs_tot = num_obs_tot + nn(2) else if (obtype == 'gps') then nobsgps = nobsgps + nn num_obs_tot = num_obs_tot + nn(2) @@ -231,6 +235,7 @@ subroutine get_num_convobs_bin(obspath,datestring,num_obs_tot,num_obs_totdiag,id write(6,100) 'dw',nobsdw(1),nobsdw(2) write(6,100) 'rw',nobsrw(1),nobsrw(2) write(6,100) 'dbz',nobsdbz(1),nobsdbz(2) + write(6,100) 'fed',nobsfed(1),nobsfed(2) write(6,100) 'tcp',nobstcp(1),nobstcp(2) if (nobstcx(2) .gt. 0) then write(6,100) 'tcx',nobstcx(1),nobstcx(2) @@ -1075,6 +1080,7 @@ subroutine get_convobs_data_bin(obspath, datestring, nobs_max, nobs_maxdiag, & if (obtype == ' t' .or. obtype == ' uv' .or. obtype == ' ps' .or. & obtype == 'tcp' .or. obtype == ' q' .or. obtype == 'spd' .or. & obtype == 'sst' .or. obtype == ' rw' .or. obtype == 'dbz' .or. & + obtype == 'fed' .or. & obtype == 'gps' .or. obtype == ' dw' .or. obtype == ' pw') then ! direct reflectivitiy DA has a different routine for dbz obs.