diff --git a/fix b/fix index 5722cd4d25..6a42a29dbb 160000 --- a/fix +++ b/fix @@ -1 +1 @@ -Subproject commit 5722cd4d2519222137c5b356bdbc01bb34c5f1f4 +Subproject commit 6a42a29dbbc9fca3453cc9e829601185555890b9 diff --git a/modulefiles/gsi_orion.lua b/modulefiles/gsi_orion.lua index e75a01ef5e..e3d4a68b11 100644 --- a/modulefiles/gsi_orion.lua +++ b/modulefiles/gsi_orion.lua @@ -19,6 +19,16 @@ load(pathJoin("hpc-impi", hpc_impi_ver)) load(pathJoin("cmake", cmake_ver)) load("gsi_common") +setenv("crtm_ROOT","/work/noaa/da/eliu/JEDI-GDAS/crtm_v2.4.1-jedi.1-intel2022/build") +setenv("crtm_VERSION","2.4.1-jedi.1") +setenv("CRTM_INC","/work/noaa/da/eliu/JEDI-GDAS/crtm_v2.4.1-jedi.1-intel2022/build/module") +setenv("CRTM_LIB","/work/noaa/da/eliu/JEDI-GDAS/crtm_v2.4.1-jedi.1-intel2022/build/lib/libcrtm_static.a") +setenv("CRTM_FIX","/work/noaa/da/eliu/JEDI-GDAS/crtm_v2.4.1-jedi.1-fix_gdasapp/fix") +whatis("Name: crtm") +whatis("Version: 2.4.1-jedi.1") +whatis("Category: library") +whatis("Description: crtm library") + load(pathJoin("prod_util", prod_util_ver)) diff --git a/src/gsi/CMakeLists.txt b/src/gsi/CMakeLists.txt index af94224c05..325bfdb6a5 100644 --- a/src/gsi/CMakeLists.txt +++ b/src/gsi/CMakeLists.txt @@ -146,6 +146,7 @@ target_link_libraries(gsi_fortran_obj PUBLIC nemsio::nemsio) target_link_libraries(gsi_fortran_obj PUBLIC ncio::ncio) target_link_libraries(gsi_fortran_obj PUBLIC w3emc::w3emc_d) target_link_libraries(gsi_fortran_obj PUBLIC sp::sp_d) +add_library(crtm::crtm ALIAS crtm) target_link_libraries(gsi_fortran_obj PUBLIC bufr::bufr_d) target_link_libraries(gsi_fortran_obj PUBLIC crtm::crtm) if(GSI_MODE MATCHES "Regional") diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index 2305c84340..1b093d44a2 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -85,6 +85,7 @@ module crtm_interface public destroy_crtm ! Subroutine destroys initialization for crtm public sensorindex public surface +public atmosphere public isatid ! = 1 index of satellite id public itime ! = 2 index of analysis relative obs time public ilon ! = 3 index of grid relative obs location (x) @@ -125,6 +126,8 @@ module crtm_interface public idtw ! = 35/37 index of d(Tw) public idtc ! = 36/38 index of d(Tc) public itz_tr ! = 37/39 index of d(Tz)/d(Tr) +public n_clouds_fwd_wk +public n_absorbers ! For TMI and GMI public iedge_log ! = 32 ! index, if obs is to be obleted beause of locating near scan edges. @@ -202,6 +205,7 @@ module crtm_interface logical ,save :: mixed_use logical ,save :: use_gfdl_qsat integer(i_kind), parameter :: min_n_absorbers = 2 + integer(i_kind) :: n_absorbers integer(i_kind),save :: iedge_log integer(i_kind),save :: ilzen_ang2,ilazi_ang2,iscan_ang2,iszen_ang2,isazi_ang2 @@ -356,7 +360,6 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo ! ...all "additional absorber" variables integer(i_kind) :: j,icount integer(i_kind) :: ig - integer(i_kind) :: n_absorbers logical quiet logical print_verbose diff --git a/src/gsi/genstats_gps.f90 b/src/gsi/genstats_gps.f90 index acf5ca2756..9d41602228 100644 --- a/src/gsi/genstats_gps.f90 +++ b/src/gsi/genstats_gps.f90 @@ -64,6 +64,15 @@ module m_gpsStats integer(i_kind) :: idv,iob ! device id and obs index for sorting real (r_kind) :: elat, elon ! earth lat-lon for redistribution !real (r_kind) :: dlat, dlon ! earth lat-lon for redistribution + ! + real(r_kind),dimension(:),pointer :: tsenges => NULL() + real(r_kind),dimension(:),pointer :: tvirges => NULL() + real(r_kind),dimension(:),pointer :: sphmges => NULL() + real(r_kind),dimension(:),pointer :: hgtlges => NULL() + real(r_kind),dimension(:),pointer :: hgtiges => NULL() + real(r_kind),dimension(:),pointer :: prslges => NULL() + real(r_kind),dimension(:),pointer :: prsiges => NULL() + ! end type gps_all_ob_type type gps_all_ob_head @@ -427,7 +436,10 @@ subroutine genstats_gps(bwork,awork,toss_gps_sub,conv_diagsave,mype) end do END DO if(icnt > 0)then - nreal =22 +! +! nreal =22 + nreal =34 +! ioff =nreal if (lobsdiagsave) nreal=nreal+4*miter+1 if (save_jacobian) then @@ -439,8 +451,6 @@ subroutine genstats_gps(bwork,awork,toss_gps_sub,conv_diagsave,mype) end if endif - - ! Loop over data to apply final qc, superobs factors, accumulate ! statistics and (optionally) load diagnostic output arrays icnt=0 @@ -534,10 +544,12 @@ subroutine genstats_gps(bwork,awork,toss_gps_sub,conv_diagsave,mype) elev = gps_allptr%rdiag(7) dobs = gps_allptr%rdiag(17) if (toss_gps(kprof) > zero .and. (dobs == toss_gps(kprof) .or. elev < dobs_height(kprof))) then ! SR from obs + if(ratio_errors*data_ier > tiny_r_kind) then ! obs was good if (luse) then if(conv_diagsave) then - gps_allptr%rdiag(10) = four + !gps_allptr%rdiag(10) = four + gps_allptr%rdiag(10) = 8 gps_allptr%rdiag(12) = -one gps_allptr%rdiag(16) = zero if(lobsdiagsave) gps_allptr%rdiag(mreal+jiter) = -one @@ -704,7 +716,6 @@ subroutine genstats_gps(bwork,awork,toss_gps_sub,conv_diagsave,mype) endif endif - ! Destroy arrays holding gps data call destroy_genstats_gps contains @@ -760,19 +771,25 @@ subroutine contents_netcdf_diag_ ! Observation class character(7),parameter :: obsclass = ' gps' - call nc_diag_metadata("Station_ID", gps_allptr%cdiag ) - call nc_diag_metadata("Observation_Class", obsclass ) + call nc_diag_metadata("Station_ID", gps_allptr%cdiag ) + call nc_diag_metadata("Observation_Class", obsclass ) obstype = gps_allptr%rdiag(1) obssubtype = gps_allptr%rdiag(2) call nc_diag_metadata("Observation_Type", obstype ) - call nc_diag_metadata("Observation_Subtype", obssubtype ) +! rename it to be consistent with GMAO gsi_ncdiag +! call nc_diag_metadata("Observation_Subtype", obssubtype ) + call nc_diag_metadata("record_number", obssubtype ) call nc_diag_metadata_to_single("Latitude", gps_allptr%rdiag(3) ) call nc_diag_metadata_to_single("Longitude", gps_allptr%rdiag(4) ) call nc_diag_metadata_to_single("Incremental_Bending_Angle", gps_allptr%rdiag(5) ) - call nc_diag_metadata_to_single("Pressure", gps_allptr%rdiag(6) ) - call nc_diag_metadata_to_single("Height", gps_allptr%rdiag(7) ) +! call nc_diag_metadata_to_single("Pressure", gps_allptr%rdiag(6) ) !orig + call nc_diag_metadata_to_single("Pressure", gps_allptr%rdiag(6)*100.0 ) +! rename the variable as impact_height +! call nc_diag_metadata_to_single("Height", gps_allptr%rdiag(7) ) + call nc_diag_metadata_to_single("Impact_Height", gps_allptr%rdiag(7) ) call nc_diag_metadata_to_single("Time", gps_allptr%rdiag(8) ) - call nc_diag_metadata_to_single("Model_Elevation", gps_allptr%rdiag(9) ) +! call nc_diag_metadata_to_single("Model_Elevation", gps_allptr%rdiag(9) ) !orig + call nc_diag_metadata_to_single("surface_geopotential_height", gps_allptr%rdiag(9) ) call nc_diag_metadata_to_single("Setup_QC_Mark", gps_allptr%rdiag(10) ) call nc_diag_metadata_to_single("Prep_Use_Flag", gps_allptr%rdiag(11) ) call nc_diag_metadata_to_single("Analysis_Use_Flag", gps_allptr%rdiag(12) ) @@ -788,14 +805,33 @@ subroutine contents_netcdf_diag_ call nc_diag_metadata_to_single("Temperature_at_Obs_Location", gps_allptr%rdiag(18) ) call nc_diag_metadata_to_single("Specific_Humidity_at_Obs_Location",gps_allptr%rdiag(21) ) + call nc_diag_metadata("impact_parameter", sngl(gps_allptr%rdiag(23)) ) + call nc_diag_metadata("pccf", sngl(gps_allptr%rdiag(24)) ) + call nc_diag_metadata("reference_sat_id", int(gps_allptr%rdiag(25)) ) + call nc_diag_metadata("earth_radius_of_curvature", sngl(gps_allptr%rdiag(26)) ) + call nc_diag_metadata("geoid_height_above_reference_ellipsoid", sngl(gps_allptr%rdiag(27)) ) + call nc_diag_metadata("qfro", int(gps_allptr%rdiag(28)) ) + call nc_diag_metadata("ascending_flag", int(gps_allptr%rdiag(29)) ) + call nc_diag_metadata("sensor_azimuth_angle", sngl(gps_allptr%rdiag(30)) ) + call nc_diag_metadata("sat_constellation", int(gps_allptr%rdiag(31)) ) + call nc_diag_metadata("occulting_sat", int(gps_allptr%rdiag(32)) ) + call nc_diag_metadata("process_center", int(gps_allptr%rdiag(33)) ) + call nc_diag_metadata("atmospheric_refractivity", sngl(gps_allptr%rdiag(34)) ) + call nc_diag_metadata("surface_altitude", sngl(gps_allptr%rdiag(9)) ) + if (save_jacobian) then call readarray(dhx_dx, gps_allptr%rdiag(ioff+1:nreal)) call nc_diag_data2d("Observation_Operator_Jacobian_stind", dhx_dx%st_ind(1:dhx_dx%nind)) call nc_diag_data2d("Observation_Operator_Jacobian_endind", dhx_dx%end_ind(1:dhx_dx%nind)) call nc_diag_data2d("Observation_Operator_Jacobian_val", real(dhx_dx%val(1:dhx_dx%nnz),r_single)) endif - - + call nc_diag_data2d("atmosphere_pressure_coordinate", sngl(gps_allptr%prslges)) + call nc_diag_data2d("atmosphere_pressure_coordinate_interface", sngl(gps_allptr%prsiges)) + call nc_diag_data2d("air_temperature", sngl(gps_allptr%tsenges)) + call nc_diag_data2d("virtual_temperature", sngl(gps_allptr%tvirges)) + call nc_diag_data2d("specific_humidity", sngl(gps_allptr%sphmges)) + call nc_diag_data2d("geopotential_height", sngl(gps_allptr%hgtlges)) + call nc_diag_data2d("geopotential_height_levels", sngl(gps_allptr%hgtiges)) ! call nc_diag_data2d("T_Jacobian", gps_allptr%mmpoint%jac_t ) if (lobsdiagsave) then diff --git a/src/gsi/gesinfo.F90 b/src/gsi/gesinfo.F90 index 9d287de414..0b5fdd7a0a 100644 --- a/src/gsi/gesinfo.F90 +++ b/src/gsi/gesinfo.F90 @@ -588,6 +588,15 @@ subroutine gesinfo time_offset=time_offset+fha(3)/r60 #endif + write(6,*)'emily time: nhr_assimilation = ', nhr_assimilation + write(6,*)'emily time: idate nmin_an = ', idate , nmin_an + write(6,*)'emily time: ibdate = ', ibdate + write(6,*)'emily time: iadatebgn = ', iadateend + write(6,*)'emily time: iedate = ', iedate + write(6,*)'emily time: iadateend = ', iedate + write(6,*)'emily time: ianldate = ', ianldate + write(6,*)'emily time: iwinbgn = ', iwinbgn + ! Get information about date/time and number of guess files if (regional) then if(wrf_nmm_regional) then diff --git a/src/gsi/qcmod.f90 b/src/gsi/qcmod.f90 index 7146ceff3e..79c3e5894f 100644 --- a/src/gsi/qcmod.f90 +++ b/src/gsi/qcmod.f90 @@ -147,7 +147,7 @@ module qcmod ! !$$$ end documentation block - use kinds, only: i_kind,r_kind,r_double + use kinds, only: i_kind,r_kind,r_double, r_single use constants, only: zero,quarter,half,one,two,three,four,five,tiny_r_kind,rd,grav use constants, only: r0_01,r0_02,r0_03,r0_04,r0_05,r10,r60,r100,h300,r400,r1000,r2000,r2400,r4000 use constants, only: deg2rad,rad2deg,t0c,one_tenth,rearth_equator @@ -193,7 +193,7 @@ module qcmod public :: ifail_iland_det, ifail_isnow_det, ifail_iice_det, ifail_iwater_det,& ifail_imix_det, ifail_iomg_det, ifail_isst_det, ifail_itopo_det,& ifail_iwndspeed_det - public :: cao_check + public :: cao_check, ifail_clrfrac_geocsr_qc !emily public :: buddycheck_t,buddydiag_save public :: vadwnd_l2rw_qc public :: pvis,pcldch,scale_cv,estvisoe,estcldchoe,vis_thres,cldch_thres @@ -247,6 +247,7 @@ module qcmod integer(i_kind),parameter:: ifail_surface_qc=5 ! Reject due to gross check in specific qc routine integer(i_kind),parameter:: ifail_gross_routine_qc=6 + integer(i_kind),parameter:: ifail_gross_routine_nonsea_qc=16 ! Reject due to cloud > limit for channel in qc routine integer(i_kind),parameter:: ifail_cloud_qc=7 ! Reject due to inaccurate emissivity/surface temperature estimate in qc routine @@ -354,7 +355,8 @@ module qcmod ! QC_geocsr ! Reject because of standard deviation in subroutine qc_geocsr - integer(i_kind),parameter:: ifail_std_geocsr_qc=50 + integer(i_kind),parameter:: ifail_std_geocsr_qc=52 !emily change from 50 to 52 + integer(i_kind),parameter:: ifail_clrfrac_geocsr_qc=51 !emily ! QC_avhrr ! Reject because of too large surface temperature physical retrieval in qc routine: tz_retrieval (see tzr_qc) @@ -1085,6 +1087,15 @@ end subroutine tz_retrieval subroutine qc_ssmi(nchanl,nsig,ich,sfchgt,luse,sea,mixed, & temp,wmix,ts,pems,ierrret,kraintype,tpwc,clw,sgagl,tzbgr, & tbc,tbcnob,tb_ges,tnoise,ssmi,amsre_low,amsre_mid,amsre_hig,ssmis, & + varinv_after_grossroutinechk_over_ocean, & !emily + varinv_after_topo, & !emily + varinv_after_sfcchk, & !emily + varinv_after_ch2chk, & !emily + varinv_after_grossroutinechk,& !emily + varinv_after_scatteringchk, & !emily + varinv_after_nsstret, & !emily + varinv_after_jsfcchk, & !emily + pred9,pred10,pred11, & !emily varinv,errf,aivals,id_qc) ! varinv,errf,aivals,id_qc,radmod) ! all-sky @@ -1197,6 +1208,15 @@ subroutine qc_ssmi(nchanl,nsig,ich,sfchgt,luse,sea,mixed, & real(r_kind),dimension(nsig,nchanl),intent(in ) :: temp,wmix real(r_kind) ,dimension(nchanl),intent(inout) :: varinv,errf + real(r_kind) ,dimension(nchanl),intent( out) :: varinv_after_grossroutinechk_over_ocean !emily + real(r_kind) ,dimension(nchanl),intent( out) :: varinv_after_grossroutinechk !emily + real(r_kind) ,dimension(nchanl),intent( out) :: varinv_after_topo !emily + real(r_kind) ,dimension(nchanl),intent( out) :: varinv_after_sfcchk !emily + real(r_kind) ,dimension(nchanl),intent( out) :: varinv_after_ch2chk !emily + real(r_kind) ,dimension(nchanl),intent( out) :: varinv_after_scatteringchk !emily + real(r_kind) ,dimension(nchanl),intent( out) :: varinv_after_nsstret !emily + real(r_kind) ,dimension(nchanl),intent( out) :: varinv_after_jsfcchk !emily + real(r_kind) ,intent( out) :: pred9,pred10,pred11 !emily real(r_kind) ,dimension(40) ,intent(inout) :: aivals ! Declare local variables @@ -1204,7 +1224,7 @@ subroutine qc_ssmi(nchanl,nsig,ich,sfchgt,luse,sea,mixed, & integer(i_kind) :: l,i real(r_kind) :: efact,vfact,dtempf,dtbf,term real(r_kind),dimension(nchanl) :: demisf_mi,clwcutofx - real(r_kind) :: pred9,pred10,pred11 +! real(r_kind) :: pred9,pred10,pred11 !emily real(r_kind) :: dtz,ts_ave,xindx,tzchks !------------------------------------------------------------------ @@ -1248,6 +1268,15 @@ subroutine qc_ssmi(nchanl,nsig,ich,sfchgt,luse,sea,mixed, & efact =one vfact =one +!>>emily + do i = 1,24 + varinv_after_grossroutinechk_over_ocean(i) = varinv(i) !emily + varinv_after_topo(i) = varinv(i) !emily + varinv_after_sfcchk(i) = varinv(i) !emily + varinv_after_ch2chk(i) = varinv(i) !emily + varinv_after_grossroutinechk(i) = varinv(i) !emily + enddo +!<= 1.5_r_kind) then ! the data at cloud-affected channels are not used do i =1,2 varinv(i) = zero @@ -1368,13 +1401,16 @@ subroutine qc_ssmi(nchanl,nsig,ich,sfchgt,luse,sea,mixed, & varinv(i) = zero if(id_qc(i)== igood_qc ) id_qc(i)=ifail_ch2_qc end do + varinv_after_ch2chk(:) = varinv(:) !emily endif !General qc criteria for all channels do i = 1,24 if( abs(tbcnob(i)) >= 3.5_r_kind) then varinv(i) = zero - if(id_qc(i)== igood_qc ) id_qc(i)=ifail_gross_routine_qc + ! if(id_qc(i)== igood_qc ) id_qc(i)=ifail_gross_routine_qc + if(id_qc(i)== igood_qc ) id_qc(i)=ifail_gross_routine_nonsea_qc !emily end if + varinv_after_grossroutinechk(i) = varinv(i) !emily enddo end if @@ -1392,6 +1428,12 @@ subroutine qc_ssmi(nchanl,nsig,ich,sfchgt,luse,sea,mixed, & end if end if +!>>emily + do i = 1, nchanl + varinv_after_scatteringchk(i) = varinv(i) !emily + enddo +!< two) then varinv(9) =zero if(id_qc(9)== igood_qc) id_qc(9)=ifail_scatt_qc + varinv_after_scatteringchk(9) = varinv(9) !emily end if if(pred10 - tbcnob(10) - tb_ges(10) > two) then varinv(10)=zero if(id_qc(10)== igood_qc) id_qc(10)=ifail_scatt_qc + varinv_after_scatteringchk(10) = varinv(10) !emily end if if(pred11 - tbcnob(11) - tb_ges(11) > two) then varinv(11)=zero if(id_qc(11)== igood_qc) id_qc(11)=ifail_scatt_qc + varinv_after_scatteringchk(11) = varinv(11) !emily end if end if ! @@ -1436,6 +1481,11 @@ subroutine qc_ssmi(nchanl,nsig,ich,sfchgt,luse,sea,mixed, & enddo endif endif +!>>emily + do l=1,nchanl + varinv_after_nsstret(:) = varinv(:) !emily + enddo +!<>emily + do i=1,nchanl + varinv_after_wavenum(i) = varinv(i) + enddo +!< r1000 .or. tb_obs(i) <= zero) then varinv(i)=zero varinv_use(i)=zero + if(id_qc(i) == igood_qc)id_qc(i)=ifail_range_qc !emily_bugfix end if + varinv_after_rangechk(i) = varinv(i) !emily varinv(i) = varinv(i)*(one-(one-sfchgtfact)*ptau5(1,i)) varinv_use(i) = varinv_use(i)*(one-(one-sfchgtfact)*ptau5(1,i)) + varinv_after_topo(i) = varinv(i) !emily ! Modify error based on transmittance at top of model varinv(i)=varinv(i)*ptau5(nsig,i) varinv_use(i)=varinv_use(i)*ptau5(nsig,i) errf(i)=errf(i)*ptau5(nsig,i) + varinv_after_transmittop(i) = varinv(i) !emily ! QC based on presence/absence of cloud sum3=sum3+tbc(i)*tbc(i)*varinv_use(i) @@ -2752,6 +2817,12 @@ subroutine qc_avhrr(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse, & end do end if +!>>emily + do i=1,nchanl + varinv_after_clddet(i) = varinv(i) !emily + enddo +!<>emily + do i=1,nchanl + varinv_after_nsstret(i) = varinv(i) !emily + enddo +!< tiny_r_kind)varinv(i)=varinv(i)/(one+varinv(i)*term) end if end do - +!>>emily + do i=1,nchanl + varinv_after_jsfcchk(i) = varinv(i) !emily + end do +!< all-sky chk +!>>emily + qc4emiss_out = zero + if (qc4emiss) qc4emiss_out = one +!<>emily + do i = 1, nchanl + varinv_grosschk(i) = varinv(i) + enddo +!< r2000) then @@ -3423,6 +3525,12 @@ subroutine qc_amsua(nchanl,is,ndat,nsig,npred,sea,land,ice,snow,mixed,luse, & end if end do +!>>emily + do i=1,nchanl + varinv_after_jsfcchk(i) = varinv(i) !emily + end do +!<>emily + do i=1,nchanl + varinv_sdoei(i) = varinv(i) !emily + enddo +!<>emily + do i=1,nchanl + varinv_after_sdoei(i) = varinv(i) !emily + end do +!<>emily + id_qc,aivals,errf,varinv,varinv_use, & + varinv_after_sfcterrianchk, & + varinv_after_rangechk, & + varinv_after_topo, & + varinv_after_transmittop, & + varinv_after_clddet, & + varinv_after_stdchk, & + varinv_after_grossroutinechk, & + varinv_after_stdadj, & + varinv_after_nsstret, & + varinv_after_jsfcchk, & + cld,cldp,kmax,abi,ahi,seviri) +!< r1000 .or. tb_obs(i) <= zero) then varinv(i)=zero varinv_use(i)=zero + if(id_qc(i) == igood_qc)id_qc(i)=ifail_range_qc !emily_bugfix end if + varinv_after_rangechk(i) = varinv(i) !emily + tmp=one-(one-sfchgtfact)*ptau5(1,i) varinv(i) = varinv(i)*tmp varinv_use(i) = varinv_use(i)*tmp + varinv_after_topo(i) = varinv(i) !emily ! Modify error based on transmittance at top of model varinv(i)=varinv(i)*ptau5(nsig,i) varinv_use(i)=varinv_use(i)*ptau5(nsig,i) errf(i)=errf(i)*ptau5(nsig,i) + varinv_after_transmittop(i) = varinv(i) !emily ! QC based on presence/absence of cloud sum3=sum3+tbc(i)*tbc(i)*varinv_use(i) @@ -4482,7 +4640,6 @@ subroutine qc_geocsr(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse, & if(id_qc(i) == igood_qc)id_qc(i)=ifail_cloud_qc end if end do - ! If no clouds check surface temperature/emissivity else ! If no cloud was detected, do surface temp/emiss checks sum=zero @@ -4511,7 +4668,19 @@ subroutine qc_geocsr(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse, & end do end if endif +!>>emily + do i = 1, nchanl + varinv_after_clddet(i) = varinv(i) !emily + end do +!<>emily + do i = 1, nchanl + varinv_after_stdchk(i) = varinv(i) !emily + varinv_after_grossroutinechk(i) = varinv(i) !emily + varinv_after_stdadj(i) = varinv(i) !emily + enddo +!<= 0.5 for chn10.3 @@ -4532,6 +4701,7 @@ subroutine qc_geocsr(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse, & varinv(i)=zero end if end if + varinv_after_stdchk(i) = varinv(i) !emily ! QC_o-g: If abs(o-g) > 2.0 do not use if ( i/=2 .and. abs(tbc(i)) > two ) then varinv(i) = zero @@ -4539,6 +4709,7 @@ subroutine qc_geocsr(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse, & ! QC1 in statsrad if(luse)aivals(8,is)= aivals(8,is) + one !hliu check end if + varinv_after_grossroutinechk(i) = varinv(i) !emily end if ! adjust varinv according to the BT standard deviation @@ -4557,7 +4728,7 @@ subroutine qc_geocsr(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse, & if(seviri .or. ahi) then varinv(i)=varinv(i) end if - + varinv_after_stdadj(i) = varinv(i) !emily end do ! ! Apply Tz retrieval @@ -4586,6 +4757,12 @@ subroutine qc_geocsr(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse, & endif end if +!>>emily + do i=1,nchanl + varinv_after_nsstret(i) = varinv(i) !emily + enddo +!>>emily + ! Generate q.c. bounds and modified variances. do i=1,nchanl if(varinv(i) > tiny_r_kind)then @@ -4593,8 +4770,8 @@ subroutine qc_geocsr(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse, & term = dtbf*dtbf if(term > tiny_r_kind)varinv(i)=varinv(i)/(one+varinv(i)*term) end if + varinv_after_jsfcchk(i) = varinv(i) !emily end do - return diff --git a/src/gsi/read_bufrtovs.f90 b/src/gsi/read_bufrtovs.f90 index 2fc14b5cdf..d07933c133 100644 --- a/src/gsi/read_bufrtovs.f90 +++ b/src/gsi/read_bufrtovs.f90 @@ -1083,6 +1083,11 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& call count_obs(ndata,nele,ilat,ilon,data_all,nobs) write(lunout) obstype,sis,nreal,nchanl,ilat,ilon write(lunout) ((data_all(k,n),k=1,nele),n=1,ndata) +! write(6,*)'emily checking jsatid = ', jsatid +! write(6,*)'emily checking nread = ', nread +! write(6,*)'emily checking ndata = ', ndata +! write(6,*)'emily checking ndata*nchanl = ', ndata*nchanl +! write(6,*)'emily checking nodata = ', nodata end if ! Deallocate local arrays diff --git a/src/gsi/read_gps.f90 b/src/gsi/read_gps.f90 index c0cef658a6..65e7343e15 100644 --- a/src/gsi/read_gps.f90 +++ b/src/gsi/read_gps.f90 @@ -137,21 +137,30 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & bend_error,ref_error,bend_pccf,ref_pccf real(r_kind),allocatable,dimension(:,:):: cdata_all + - integer(i_kind),parameter:: n1ahdr=10 +! integer(i_kind),parameter:: n1ahdr=10 + integer(i_kind),parameter:: n1ahdr=13 + real(r_double),dimension(n1ahdr):: bfr1ahdr real(r_double),dimension(50,maxlevs):: data1b real(r_double),dimension(50,maxlevs):: data2a real(r_double),dimension(maxlevs):: nreps_this_ROSEQ2 + + real(r_kind):: azm_ang, sat_ascd, sat_constid, siid, ogce data lnbufr/10/ - data hdr1a / 'YEAR MNTH DAYS HOUR MINU PCCF ELRC SAID PTID GEODU' / +! data hdr1a / 'YEAR MNTH DAYS HOUR MINU PCCF ELRC SAID PTID GEODU' / + data hdr1a / 'YEAR MNTH DAYS HOUR MINU PCCF ELRC SAID PTID GEODU SCLF SIID OGCE' / data nemo /'QFRO'/ !*********************************************************************************** maxobs=2e6 - nreal=maxinfo +!> xuanli +! nreal=maxinfo + nreal=23 +!< xuanli nchanl=0 ilon=2 ilat=3 @@ -170,7 +179,6 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & write(6,*)'READ GPS: CONVINFO DOES NOT INCLUDE ANY ',trim(sis),' DATA' return end if - ! Open file for input, then read bufr data open(lnbufr,file=trim(infile),form='unformatted') call openbf(lnbufr,'IN',lnbufr) @@ -214,16 +222,19 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & call ufbint(lnbufr,qfro,1,1,iret,nemo) ! observation time in minutes - idate5(1) = bfr1ahdr(1) ! year - idate5(2) = bfr1ahdr(2) ! month - idate5(3) = bfr1ahdr(3) ! day - idate5(4) = bfr1ahdr(4) ! hour - idate5(5) = bfr1ahdr(5) ! minute - pcc=bfr1ahdr(6) ! profile per cent confidence - roc=bfr1ahdr(7) ! Earth local radius of curvature - said=bfr1ahdr(8) ! Satellite identifier - ptid=bfr1ahdr(9) ! Platform transmitter ID number - geoid=bfr1ahdr(10) ! Geoid undulation + idate5(1) = bfr1ahdr(1) ! year + idate5(2) = bfr1ahdr(2) ! month + idate5(3) = bfr1ahdr(3) ! day + idate5(4) = bfr1ahdr(4) ! hour + idate5(5) = bfr1ahdr(5) ! minute + pcc=bfr1ahdr(6) ! profile per cent confidence + roc=bfr1ahdr(7) ! Earth local radius of curvature + said=bfr1ahdr(8) ! Satellite identifier + ptid=bfr1ahdr(9) ! Platform transmitter ID number + geoid=bfr1ahdr(10) ! Geoid undulation + sat_constid=bfr1ahdr(11) ! Satellite classification + siid=bfr1ahdr(12) ! Satellite instrument + ogce = bfr1ahdr(13) ! Identification of originating/generating centre call w3fs21(idate5,minobs) ! Locate satellite id in convinfo file @@ -296,6 +307,15 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & endif endif +! ascending flag: when qfro bit3 is set, occultation is ascending +! bit3 is clear, occultation is descending + sat_ascd = 0.0 + call upftbv(lnbufr,nemo,qfro,mxib,ibit,nib) + if(nib > 0) then + do i=1,nib + if(ibit(i) .eq. 3) sat_ascd=1.0 + enddo + endif ! Read bending angle information ! Get the number of occurences of sequence ROSEQ2 in this subset @@ -345,6 +365,7 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & nread=nread+1 ! count observations rlat=data1b(1,k) ! earth relative latitude (degrees) rlon=data1b(2,k) ! earth relative longitude (degrees) + azm_ang=data1b(3,k) ! azimuth angle !xuanli height=data2a(1,k) ref=data2a(2,k) ref_error=data2a(4,k) @@ -368,7 +389,8 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & good=.true. if((abs(rlat)>90._r_kind).or.(abs(rlon)>r360).or.(height<=zero)) then good=.false. - else if (ref_obs) then + endif + if (ref_obs) then if ((ref>=1.e+9_r_kind).or.(ref<=zero).or.(height>=1.e+9_r_kind)) then good=.false. endif @@ -440,6 +462,14 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & cdata_all(15,ndata)= dlat_earth_deg ! earth relative latitude (degrees) cdata_all(16,ndata)= geoid ! geoid undulation (m) + cdata_all(17,ndata)= qfro ! qfro + cdata_all(18,ndata)= sat_ascd ! ascending flag + cdata_all(19,ndata)= azm_ang ! azimuth angle + cdata_all(20,ndata)= sat_constid ! satellite classification + cdata_all(21,ndata)= siid ! occulting satellite + cdata_all(22,ndata)= ogce ! Identification of originating/generating centre + cdata_all(23,ndata)= ref ! refractivity obs (units of N) + else notgood = notgood + 1 end if @@ -465,9 +495,8 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & write(6,*)'READ_GPS: # bad or missing data=', notgood do i=1,ngpsro_type if (nmrecs_id(i)>0) & - write(6,1021)'READ_GPS: LEO_id,nprof_gps = ',gpsro_itype(i),nmrecs_id(i) + write(6,1020)'READ_GPS: LEO_id,nprof_gps = ',gpsro_itype(i),nmrecs_id(i) end do -1021 format(a31,i6,i6) write(6,1020)'READ_GPS: ref_obs,nprof_gps= ',ref_obs,nprof_gps 1020 format(a31,L,i6) diff --git a/src/gsi/read_ozone.f90 b/src/gsi/read_ozone.f90 index 2ad95a858e..82725a9750 100644 --- a/src/gsi/read_ozone.f90 +++ b/src/gsi/read_ozone.f90 @@ -154,6 +154,7 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & integer(i_kind) JULIAN,IDAYYR,IDAYWK integer(i_kind) ikx integer(i_kind) decimal,binary(14),binary_mls(18) + integer(i_kind) :: iuseflag !emily 1: used; -1: not used integer(i_kind) itx,itt,ipoq7 @@ -220,7 +221,8 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & !************************************************************************** ! Set constants. Initialize variables rsat=999._r_kind - maxobs=1e6 + maxobs=1e6 !orig +! maxobs=1e7 !emily ilon=3 ilat=4 ipoq7=0 @@ -231,7 +233,8 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & if (obstype == 'sbuv2' .or. obstype == 'ompsnp') then - nreal=9 +! nreal=9 !orig + nreal=10 !emily open(lunin,file=trim(infile),form='unformatted') nmrecs=0 call openbf(lunin,'IN',lunin) @@ -272,6 +275,7 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & endif read_loop1: do + iuseflag = 1 !emily call readsb(lunin,iret) if (iret/=0) then call readmg(lunin,subset,jdate,iret) @@ -396,15 +400,20 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & if (version8) then call ufbint(lunin,toq,1,1,iret,'SBUVTOQ') call ufbint(lunin,poq,1,1,iret,'SBUVPOQ') - if (toq/=0 .and. toq/=2) cycle read_loop1 - if (poq/=0 .and. poq/=1 .and. poq/=ipoq7) cycle read_loop1 +!!emily if (toq/=0 .and. toq/=2) cycle read_loop1 +!!emily if (poq/=0 .and. poq/=1 .and. poq/=ipoq7) cycle read_loop1 +!>>emily + if (toq/=0 .and. toq/=2) iuseflag = -1 + if (poq/=0 .and. poq/=1 .and. poq/=ipoq7) iuseflag = -1 +!<badoz) cycle read_loop1 +!!emily if (poz(k)>badoz) cycle read_loop1 + if (poz(k)>badoz) iuseflag = -1 end do - + ! Write ozone record to output file ndata=min(ndata+1,maxobs) nodata=nodata+nloz+1 @@ -417,8 +426,10 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & ozout(7,ndata)=toq ! total ozone error flag ozout(8,ndata)=poq ! profile ozone error flag ozout(9,ndata)=solzen ! solar zenith angle + ozout(10,ndata)=iuseflag ! solar zenith angle !emily do k=1,nloz+1 - ozout(k+9,ndata)=poz(k) +! ozout(k+9,ndata)=poz(k) !orig + ozout(k+10,ndata)=poz(k) !emily end do ! Loop back to read next profile @@ -624,7 +635,10 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & call makegrids(rmesh,ithin,n_tbin=n_tbin) ! Set dependent variables and allocate arrays - nreal=14 +! nreal=14 !orig + nreal=15 !emily (add AFBO) + nreal=16 !emily (add iuseflag) + iuseflag=1 !emily nloz=0 nchanl=1 nozdat=nreal+nchanl @@ -641,6 +655,7 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & ihh=0 read_loop2: do + iuseflag = 1 call readsb(lunin,iret) if (iret/=0) then call readmg(lunin,subset,jdate,iret) @@ -702,7 +717,8 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & ! extract total ozone call ufbint(lunin,totoz,1,1,iret,'OZON') - if (totoz > badoz ) cycle read_loop2 + if (totoz > badoz ) iuseflag = -1 !emily +!!emily if (totoz > badoz ) cycle read_loop2 ! QC for omi_aura if (obstype == 'omi') then @@ -710,18 +726,22 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & ! Bit 10 in TOQF represents row anomaly. decimal=int(hdrozo2(6)) call dec2bin(decimal,binary,14) - if (binary(10) == 1 ) cycle read_loop2 + if (binary(10) == 1 ) iuseflag = -1 !emily +!!emily if (binary(10) == 1 ) cycle read_loop2 ! remove the bad scan position data: fovn beyond 25 - if (hdrozo2(7) >=25.0_r_double) cycle read_loop2 + if (hdrozo2(7) >=25.0_r_double) iuseflag = -1 !emily +!!emily if (hdrozo2(7) >=25.0_r_double) cycle read_loop2 end if ! only accept flag 0 1, flag 2 is high SZA data which is not used for now toq=hdrozo2(5) - if (toq/=0 .and. toq/=1) cycle read_loop2 + if (toq/=0 .and. toq/=1) iuseflag = -1 !emily +!!emily if (toq/=0 .and. toq/=1) cycle read_loop2 ! remove the data in which the C-pair algorithm ((331 and 360 nm) is used. - if (hdrozo2(8) == 3_r_double .or. hdrozo2(8) == 13_r_double) cycle read_loop2 + if (hdrozo2(8) == 3_r_double .or. hdrozo2(8) == 13_r_double) iuseflag = -1 !emily +!!emily if (hdrozo2(8) == 3_r_double .or. hdrozo2(8) == 13_r_double) cycle read_loop2 ! thin OMI/OMPS-NM(or TC8) data @@ -746,14 +766,18 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & ozout(7,itx)=hdrozo2(5) ! total ozone quality code ozout(8,itx)=hdrozo(10) ! solar zenith angle if (obstype == 'omi') then - ozout(9,itx)=binary(10) ! row anomaly flag +! ozout(9,itx)=binary(10) ! row anomaly flag !orig + ozout(9,itx)=hdrozo2(6) ! total ozone quality flag !emily end if ozout(10,itx)=hdrozo2(1) ! cloud amount ozout(11,itx)=hdrozo2(4) ! vzan ozout(12,itx)=hdrozo2(2) ! aerosol index ozout(13,itx)=hdrozo2(3) ! ascending/descending ozout(14,itx)=hdrozo2(7) ! scan position - ozout(15,itx)=totoz + ozout(15,itx)=hdrozo2(8) ! AFBO !emily + ozout(16,itx)=iuseflag !emily + ozout(17,itx)=totoz !emily +!orig ozout(15,itx)=totoz ! End of loop over observations end do read_loop2 @@ -1252,7 +1276,6 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & ndata=kk nodata=ndata endif - ! Write header record and data to output file for further processing call count_obs(ndata,nozdat,ilat,ilon,ozout,nobs) write(lunout) obstype,sis,nreal,nchanl,ilat,ilon diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90 index eaece05451..1b4358838d 100644 --- a/src/gsi/read_prepbufr.f90 +++ b/src/gsi/read_prepbufr.f90 @@ -478,11 +478,13 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& if(tob)then nreal=25 else if(uvob) then - nreal=26 + nreal=34 +! nreal=26 !orig else if(spdob) then nreal=24 else if(psob) then - nreal=20 +! nreal=20 !orig + nreal=21 !emily else if(qob) then nreal=26 else if(pwob) then @@ -635,6 +637,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! Time offset if(nmsg == 0) call time_4dvar(idate,toff) +! write(6,*)'emily checking toff = ', toff nmsg=nmsg+1 if (nmsg>nmsgmax) then write(6,*)'READ_PREPBUFR: messages exceed maximum ',nmsgmax @@ -1049,6 +1052,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& !------------------------------------------------------------------------ +! write(6,*)'emily checking offtime_data = ', offtime_data if(offtime_data) then ! in time correction for observations to account for analysis @@ -1079,6 +1083,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& if (.not. (aircraft_t_bc .and. acft_profl_file)) then timeobs=real(real(hdr(4),r_single),r_double) t4dv=timeobs + toff + zeps=1.0e-8_r_kind if (t4dv -zeps) t4dv=zero if (t4dv>winlen.and.t4dv>emily +!>>orig +! if (isflg /= 0) cycle loop_readsb +! if (tsavg <= 273.0_r_kind) cycle loop_readsb +!<179.and.kx<190).or.(kx>=280.and.kx<=290).or. & @@ -1627,10 +1636,11 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! If temperature ob, extract information regarding virtual ! versus sensible temperature - if(tob) then +! if(tob) then !orig + if(tob .or. psob ) then !emily ! use tvirtual if tsensible flag not set, and not in either 2Dregional or global_2m DA mode if ( (.not. tsensible) .and. .not. (twodvar_regional .or. global_2m_land) ) then - + !write(6,*)'emily checking vtcd = ', vtcd do k=1,levs tvflg(k)=one ! initialize as sensible do j=1,20 @@ -1803,6 +1813,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& if ( obsdat(1,k)< r500) qm=100 zqm=idnint(qcmark(4,k)) if (zqm>=lim_zqm .and. zqm/=15 .and. zqm/=9) qm=9 + ! if (zqm>=lim_zqm .and. zqm/=15 .and. zqm/=9) pqm(k)=9 !emily_test endif ! if(convobs .and. pqm(k) >=lim_qm .and. qm/=15 .and. qm/=9 )cycle loop_k_levs @@ -2320,6 +2331,9 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! Surface pressure else if(psob) then + qtflg=tvflg(k) !emily + write(6,*)'emily checkint psob t4dv = ', t4dv + poe=obserr(1,k)*one_tenth ! convert from mb to cb if (inflate_error) poe=poe*r1_2 cdata_all(1,iout)=poe ! surface pressure error (cb) @@ -2346,7 +2360,9 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& cdata_all(18,iout)=r_prvstg(1,1) ! provider name cdata_all(19,iout)=r_sprvstg(1,1) ! subprovider name cdata_all(20,iout)=var_jb(1,k) ! non linear qc b parameter - if(perturb_obs)cdata_all(21,iout)=ran01dom()*perturb_fact ! ps perturbation + cdata_all(21,iout)=qtflg ! qtflg=0 (virtual) qtflg=1 (sensible) !emily +! if(perturb_obs)cdata_all(21,iout)=ran01dom()*perturb_fact ! ps perturbation !orig + if(perturb_obs)cdata_all(22,iout)=ran01dom()*perturb_fact ! ps perturbation if (twodvar_regional) & call adjust_error(cdata_all(14,iout),cdata_all(15,iout),cdata_all(11,iout),cdata_all(1,iout)) diff --git a/src/gsi/read_satwnd.f90 b/src/gsi/read_satwnd.f90 index 943cf4d47b..0077b4a867 100644 --- a/src/gsi/read_satwnd.f90 +++ b/src/gsi/read_satwnd.f90 @@ -272,7 +272,8 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis ! Set lower limits for observation errors werrmin=one nsattype=0 - nreal=26 + nreal=34 +! nreal=26 !orig if(perturb_obs ) nreal=nreal+2 ntread=1 ntmatch=0 @@ -1665,13 +1666,26 @@ subroutine read_satwnd(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,sis cdata_all(23,iout)=r_sprvstg(1,1) ! subprovider name cdata_all(25,iout)=var_jb ! non linear qc parameter cdata_all(26,iout)=one ! hilbert curve weight + ! extra variables for satwind qc for brett + cdata_all(27,iout)=obsdat(5) ! AMVQ for GOES-17 mitig.AMVs + cdata_all(28,iout)=hdrdat(9) ! wind computation method + cdata_all(29,iout)=hdrdat(10) ! satellite zenith angle + cdata_all(30,iout)=hdrdat(1) ! satellite identifier + cdata_all(31,iout)=qifn ! QI without forecast + cdata_all(32,iout)=qify ! QI with forecast + cdata_all(33,iout)=ee ! expected error + cdata_all(34,iout)=pct1 if(perturb_obs)then - cdata_all(27,iout)=ran01dom()*perturb_fact ! u perturbation - cdata_all(28,iout)=ran01dom()*perturb_fact ! v perturbation + cdata_all(35,iout)=ran01dom()*perturb_fact ! u perturbation + cdata_all(36,iout)=ran01dom()*perturb_fact ! v perturbation +!>>orig +! if(perturb_obs)then +! cdata_all(27,iout)=ran01dom()*perturb_fact ! u perturbation +! cdata_all(28,iout)=ran01dom()*perturb_fact ! v perturbation +!< obsLL(:) @@ -328,6 +343,21 @@ subroutine setupbend(obsLL,odiagLL, & ilone=14 ! index of earth relative longitude (degrees) ilate=15 ! index of earth relative latitude (degrees) igeoid=16 ! index of geoid undulation (a value per profile, m) +! + iqfro=17 ! index of qfro (integer) + iascd=18 ! index of ascending flag (integer) + iazm=19 ! index of azimuth angle + iconstid=20 ! index of classification ID (integer) + isiid=21 ! index of occulting sat (integer) + iogce=22 ! index of identification of originating (integer) + iref=23 ! index of refractivity + qc_superr(:) = 0 + qc_stat(:) = 0 + qc_largeBA(:) = 0 + qc_metop(:) = 0 + qc_layer(:) = 0 + qc_ddnj(:) = 0 + kprof = 0 ! initialize kprof otherwise it is assigned as 787086608 ! Intialize variables nsig_up=nsig+nsig_ext ! extend nsig_ext levels above interface level nsig @@ -345,7 +375,8 @@ subroutine setupbend(obsLL,odiagLL, & allocate(ddnj(grids_dim),grid_s(grids_dim),ref_rad_s(grids_dim)) ! Allocate arrays for output to diagnostic file - mreal=22 +! mreal=22 ! + mreal=34 ! more arrays nreal=mreal if (lobsdiagsave) nreal=nreal+4*miter+1 if (save_jacobian) then @@ -451,7 +482,20 @@ subroutine setupbend(obsLL,odiagLL, & call tintrp2a11(ges_z,zsges,dlat,dlon,dtime,hrdifsig,& mype,nfldsig) - prsltmp_o(1:nsig,i)=prsltmp(1:nsig) ! needed in minimization + prsltmp_o(1:nsig,i)=prsltmp(1:nsig) ! needed in minimization + + call tintrp2a1(ges_tsen, Tsen(1:nsig,i), dlat,dlon,dtime,hrdifsig, & + nsig, mype,nfldsig) + call tintrp2a1(geop_hgtl, hgtl(1:nsig,i), dlat,dlon,dtime,hrdifsig, & + nsig, mype,nfldsig) + call tintrp2a1(ges_lnprsl,prslnl(1:nsig,i),dlat,dlon,dtime,hrdifsig, & + nsig, mype,nfldsig) + + Tvir(1:nsig,i) = tges(1:nsig) ! virtual temperature + sphm(1:nsig,i) = qges(1:nsig) ! specific humidity + hgtl(1:nsig,i) = hgtl(1:nsig,i) + zsges ! mid level geopotential height + hgti(1:nsig+1,i) = hges(1:nsig+1) + zsges ! interface level geopotential height + prslni(1:nsig+1,i) = prsltmp(1:nsig+1) ! interface level log(pressure) ! Compute refractivity index-radius product at interface ! @@ -522,6 +566,7 @@ subroutine setupbend(obsLL,odiagLL, & do k=nsigstart,1,-1 ! check for model SR layer at obs location grad_mod=1000.0_r_kind*(nrefges(k+1,i)-nrefges(k,i))/(rges(k+1,i)-rges(k,i)) + if (abs(grad_mod)>= half*crit_grad) then ! SR - likely, to be used in obs SR qc qc_layer_SR=.true. !SR-likely layer detected endif @@ -573,7 +618,10 @@ subroutine setupbend(obsLL,odiagLL, & rdiagbuf(2,i) = data(iprof,i) ! profile identifier rdiagbuf(3,i) = data(ilate,i) ! lat in degrees rdiagbuf(4,i) = data(ilone,i) ! lon in degrees +!> modified imph in the diag file. In jedi: imph=impp-roc-geoid rdiagbuf(7,i) = tpdpres(i)-rocprof ! impact height in meters +! rdiagbuf(7,i) = tpdpres(i)-rocprof-unprof ! impact height in meters -- defination in JEDI + ! rdiagbuf(7,i) = tpdpres(i) ! impact parameter in meters rdiagbuf(8,i) = dtime-time_offset ! obs time (hours relative to analysis time) ! rdiagbuf(9,i) = data(ipctc,i) ! input bufr qc - index of per cent confidence @@ -582,6 +630,18 @@ subroutine setupbend(obsLL,odiagLL, & rdiagbuf(17,i) = data(igps,i) ! bending angle observation (radians) rdiagbuf(19,i) = hob ! model vertical grid (interface) if monotone grid rdiagbuf(22,i) = 1.e+10_r_kind ! spread (filled in by EnKF) + rdiagbuf(23,i) = tpdpres(i) ! impact parameter in meters + rdiagbuf(24,i) = data(ipctc,i) ! input bufr qc - index of per cent confidence + rdiagbuf(25,i) = data(iptid,i) ! transmitter occ id + rdiagbuf(26,i) = rocprof ! local radius of curvature (m) + rdiagbuf(27,i) = unprof ! geoid undulation (m) + rdiagbuf(28,i) = data(iqfro,i) ! qfro + rdiagbuf(29,i) = data(iascd,i) ! ascending flag + rdiagbuf(30,i) = data(iazm,i) ! azimuth angle + rdiagbuf(31,i) = data(iconstid,i) ! satellite classification + rdiagbuf(32,i) = data(isiid,i) ! occulting satellite + rdiagbuf(33,i) = data(iogce,i) ! Identification of processing center + rdiagbuf(34,i) = data(iref,i) ! refractivity if(ratio_errors(i) > tiny_r_kind) then ! obs inside model grid @@ -590,6 +650,7 @@ subroutine setupbend(obsLL,odiagLL, & if ((tpdpres(i)ref_rad(top_layer_SR+5)) then ! obs above SR/close-to-SR layer qcfail(i)=.false. if(hob < top_layer_SR+1) then !location might be aliased to the lower section of the non-monotonicity @@ -601,6 +662,7 @@ subroutine setupbend(obsLL,odiagLL, & endif else ! obs inside model SR/shadow or close-to-SR layer qcfail(i)=.true. + qc_superr(i)=1 endif endif @@ -608,6 +670,7 @@ subroutine setupbend(obsLL,odiagLL, & if ( data(igps,i) >= 0.03_r_kind .and. qc_layer_SR) then kprof = data(iprof,i) toss_gps_sub(kprof) = max (toss_gps_sub(kprof),data(igps,i)) + qc_layer(i)=1 endif endif @@ -624,8 +687,10 @@ subroutine setupbend(obsLL,odiagLL, & qrefges=qges_o(k1)*(one-delz)+qges_o(k2)*delz !Lidia rdiagbuf( 6,i) = ten*exp(dpressure) ! pressure at obs location (hPa) if monotone grid + ! atmosphere_pressure_coordinate rdiagbuf(18,i) = trefges ! temperature at obs location (Kelvin) if monotone grid rdiagbuf(21,i) = qrefges ! specific humidity at obs location (kg/kg) if monotone grid + commdat=.false. if (data(isatid,i)>=265 .and. data(isatid,i)<=269) commdat=.true. if (.not. qcfail(i)) then ! not SR @@ -716,7 +781,6 @@ subroutine setupbend(obsLL,odiagLL, & enddo muse(i)=nint(data(iuse,i)) <= jiter - ! Get refractivity index-radius and [d(ln(n))/dx] in new grid. intloop: do j=1,grids_dim ref_rad_s(j)=sqrt(grid_s(j)*grid_s(j)+tpdpres(i)*tpdpres(i)) !x_j @@ -752,8 +816,15 @@ subroutine setupbend(obsLL,odiagLL, & ihob=ihob-1 endif ddnj(j)=dot_product(dw4,nrefges(ihob-1:ihob+2,i))!derivative (dN/dx)_j + if(ddnj(j)>zero) then +!-- set dbend, rdiagbuf(5) as JEDI missing + dbend=-3.368795e+38 + rdiagbuf( 5,i) = -3.368795e+38 +! rdiagbuf( 17,i) = -3.368795e+38 +! qcfail(i)=.true. + qc_ddnj(i)=1 data(ier,i) = zero ratio_errors(i) = zero muse(i)=.false. @@ -778,6 +849,11 @@ subroutine setupbend(obsLL,odiagLL, & end do intloop if (obs_check) then ! reject observation +! -- set dbend, rdiagbuf(5) as JEDI missing + rdiagbuf( 5,i) = -3.368795e+38 +! rdiagbuf( 17,i) = -3.368795e+38 + dbend = -3.368795e+38 +! qcfail(i)=.true. data(ier,i) = zero ratio_errors(i) = zero @@ -800,9 +876,15 @@ subroutine setupbend(obsLL,odiagLL, & ! Remove very large simulated values if(dbend > 0.05_r_kind) then +! -- set dbend, rdiagbuf(5) as JEDI missing + rdiagbuf( 5,i) = -3.368795e+38 +! rdiagbuf( 17,i) = -3.368795e+38 + dbend = -3.368795e+38 +! data(ier,i) = zero ratio_errors(i) = zero qcfail(i)=.true. + qc_largeBA(i)=1 muse(i)=.false. endif @@ -871,6 +953,7 @@ subroutine setupbend(obsLL,odiagLL, & if(abs(rdiagbuf(5,i)) > cutoff) then qcfail(i)=.true. + qc_stat(i)=1 data(ier,i) = zero ratio_errors(i) = zero muse(i) = .false. @@ -890,6 +973,7 @@ subroutine setupbend(obsLL,odiagLL, & if( (alt <= eight) .and. & ((data(isatid,i)==4).or.(data(isatid,i)==3).or.(data(isatid,i)==5))) then qcfail(i)=.true. + qc_metop(i)=1 data(ier,i) = zero ratio_errors(i) = zero muse(i)=.false. @@ -951,6 +1035,12 @@ subroutine setupbend(obsLL,odiagLL, & if(qcfail(i)) rdiagbuf(10,i) = four !modified in genstats due to toss_gps_sub if(qcfail_loc(i) == one) rdiagbuf(10,i) = one if(qcfail_high(i) == one) rdiagbuf(10,i) = two + if(qc_superr(i) == 1) rdiagbuf(10,i) = 7 ! print out SR + !if(qc_layer(i) == 1) rdiagbuf(10,i) = 8 ! print out SR2, this is marked in genstats_gps + if(qc_stat(i) == 1) rdiagbuf(10,i) = 10 ! print out cutoff + if(qc_ddnj(i) == 1) rdiagbuf(10,i) = 9 ! print out ddnj>0 + if(qc_largeBA(i) == 1) rdiagbuf(10,i) = 11 ! print out large BA + if(qc_metop(i) == 1) rdiagbuf(10,i) = 12 ! print out metop 8 km if(muse(i)) then ! modified in genstats_gps due to toss_gps_sub rdiagbuf(12,i) = one ! minimization usage flag (1=use, -1=not used) @@ -1027,6 +1117,28 @@ subroutine setupbend(obsLL,odiagLL, & gps_alltail(ibin)%head%elat= data(ilate,i) gps_alltail(ibin)%head%elon= data(ilone,i) +! 2 dimensional geovals for JEDI + allocate(gps_alltail(ibin)%head%tvirges(nsig),stat=istatus) + allocate(gps_alltail(ibin)%head%tsenges(nsig),stat=istatus) + allocate(gps_alltail(ibin)%head%sphmges(nsig),stat=istatus) + allocate(gps_alltail(ibin)%head%hgtlges(nsig),stat=istatus) + allocate(gps_alltail(ibin)%head%hgtiges(nsig+1),stat=istatus) + allocate(gps_alltail(ibin)%head%prsiges(nsig+1),stat=istatus) + allocate(gps_alltail(ibin)%head%prslges(nsig),stat=istatus) + + do j= 1, nsig + gps_alltail(ibin)%head%tvirges(j) = Tvir(j,i) + gps_alltail(ibin)%head%tsenges(j) = Tsen(j,i) + gps_alltail(ibin)%head%sphmges(j) = sphm(j,i) + gps_alltail(ibin)%head%hgtlges(j) = hgtl(j,i) + gps_alltail(ibin)%head%prslges(j) = 1000.0*exp(prslnl(j,i)) + end do + + do j= 1, nsig + 1 + gps_alltail(ibin)%head%hgtiges(j) = hgti(j,i) + gps_alltail(ibin)%head%prsiges(j) = 1000.0*exp(prslni(j,i)) + end do + allocate(gps_alltail(ibin)%head%rdiag(nreal),stat=istatus) if (istatus/=0) write(6,*)'SETUPBEND: allocate error for gps_alldiag, istatus=',istatus @@ -1268,6 +1380,7 @@ subroutine setupbend(obsLL,odiagLL, & gps_alltail(ibin)%head%dataerr = data(ier,i)*data(igps,i) gps_alltail(ibin)%head%muse = muse(i) ! logical endif ! (last_pass) + end do ! i=1,nobs deallocate(ddnj,grid_s,ref_rad_s) ! Release memory of local guess arrays diff --git a/src/gsi/setupoz.f90 b/src/gsi/setupoz.f90 index b16a33b414..04b76b3629 100644 --- a/src/gsi/setupoz.f90 +++ b/src/gsi/setupoz.f90 @@ -205,7 +205,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& real(r_kind),dimension(nlevs):: pobs,gross,tnoise real(r_kind),dimension(nreal+nlevs,nobs):: data real(r_kind),dimension(nsig+1)::prsitmp - real(r_kind),dimension(nsig)::ozgestmp + real(r_kind),dimension(nsig)::ozgestmp ! GeoVaLs for JEDI/UFO real(r_single),dimension(nlevs):: pob4,grs4,err4 real(r_single),dimension(ireal,nobs):: diagbuf real(r_single),allocatable,dimension(:,:,:)::rdiagbuf @@ -223,7 +223,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& integer(i_kind) k1,k2,k,j,nz,jc,idia,irdim1,istatus,ioff0,ioff1 integer(i_kind) ioff,itoss,ikeep,ierror_toq,ierror_poq integer(i_kind) isolz,ifovn,itoqf - integer(i_kind) mm1,itime,ilat,ilon,ilate,ilone,itoq,ipoq + integer(i_kind) mm1,itime,ilat,ilon,ilate,ilone,itoq,ipoq,iafbo,iuseflag !emily integer(i_kind),dimension(iint,nobs):: idiagbuf integer(i_kind),dimension(nlevs):: ipos,iouse,ikeepk @@ -258,6 +258,9 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& call init_vars_ mm1=mype+1 + + write(6,*)'emily checking: you are here ...', myname, obstype + ! !********************************************************************************* ! Initialize arrays @@ -303,6 +306,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& end do nlev=jc + write(6,*)'emily nlev nlevs = ', obstype, nlev, nlevs ! Handle error conditions if (nlevs>nlev) write(6,*)'SETUPOZLAY: level number reduced for ',obstype,' ', & nlevs,' --> ',nlev @@ -352,6 +356,9 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& isolz=8 ! index of solar zenith angle (gome and omi only) itoqf=9 ! index of row anomaly (omi only) ifovn=14 ! index of scan position (gome and omi only) + iafbo=15 ! index of algorithm flag for best ozone (for omi, ompsnm, and ompstc8) !emily + iuseflag=16 ! index of useflag !emily + if (obstype == 'ompsnp') iuseflag=10 !emily ! If requested, save data for diagnostic ouput @@ -369,7 +376,12 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& dlat=data(ilat,i) dlon=data(ilon,i) dtime=data(itime,i) - + +!>>emily_test + ierror_toq = nint(data(itoq,i)) + ierror_poq = nint(data(ipoq,i)) +!<0) then ! write(6,*)'setupozlay: nobskeep',nobskeep @@ -416,7 +428,7 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& ozp_omi(nloz_omi) = prsitmp(1) call grdcrd1(ozp_omi(nloz_omi),prsitmp,nsig+1,-1) end if - + ! GeoVaLs for JEDI/UFO call tintrp2a1(ges_oz,ozgestmp,dlat,dlon,dtime,hrdifsig,& nsig,mype,nfldsig) @@ -431,16 +443,18 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& nlevs,mype,doz_dz) endif +!>>orig ! Check scan position errors in ompstc8 - if(obstype == "ompstc8") then - if(data(ifovn,i) == 1 .or. data(ifovn,i) == 2 .or. & - data(ifovn,i) == 3 .or. data(ifovn,i) == 4 .or. & - data(ifovn,i) == 35) then - if(abs(data(ilate,i)) > 50._r_kind)then - luse(i) = .false. - endif - endif - endif +! if(obstype == "ompstc8") then +! if(data(ifovn,i) == 1 .or. data(ifovn,i) == 2 .or. & +! data(ifovn,i) == 3 .or. data(ifovn,i) == 4 .or. & +! data(ifovn,i) == 35) then +! if(abs(data(ilate,i)) > 50._r_kind)then +! luse(i) = .false. +! endif +! endif +! endif +!<>emily + if(data(iuseflag,i) < 0) then + varinv3(k)=zero + ratio_errors(k)=zero + rat_err2=zero + if(luse(i))stats_oz(2,j) = stats_oz(2,j) + one ! number of obs tossed + endif +! Check scan position errors in ompstc8 + if(obstype == "ompstc8") then + if(data(ifovn,i) == 1 .or. data(ifovn,i) == 2 .or. & + data(ifovn,i) == 3 .or. data(ifovn,i) == 4 .or. & + data(ifovn,i) == 35) then + if(abs(data(ilate,i)) > 50._r_kind)then + varinv3(k)=zero + ratio_errors(k)=zero + rat_err2=zero + if(luse(i))stats_oz(2,j) = stats_oz(2,j) + one ! number of obs tossed + endif + endif + endif +!<tiny_r_kind .or. & @@ -550,6 +585,11 @@ subroutine setupozlay(obsLL,odiagLL,lunin,mype,stats_oz,nlevs,nreal,nobs,& if (rat_err2*varinv3(k)>tiny_r_kind .and. luse(i)) & stats_oz(7,j) = stats_oz(7,j) + one +!>>emily_added + varinv4diag(k)=varinv3(k) + rat_err4diag=rat_err2 +!<>emily + if (obstype == 'omi' .or. obstype == 'ompstc8' .or. obstype == 'ompsnm') then + call nc_diag_metadata("Algorithm_Flag_For_Best_Ozone", sngl(data(iafbo,i))) + endif +!< 0)then do i=1,nobs @@ -299,9 +316,31 @@ subroutine setupps(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsa ! Check for duplicate observations at same location dup=one do k=1,nobs + rstation_id = data(id,k) + ikx = nint(data(ikxx,k)) + itype = ictype(ikx) + write(10000+mype, *) ' ' + write(10000+mype, *) 'new obs ... ' + write(10000+mype, *) 'k lat/lon D ', k, data(ilate,k), data(ilone,k) + write(10000+mype, *) 'k lat/lon S ', k, slat(k), slon(k) + write(10000+mype, *) 'k station_id ', k, station_id + write(10000+mype, *) 'k itype ', k, itype + write(10000+mype, *) 'k time ', k, data(itime,k) + write(10000+mype, *) 'k obserr ', k, data(ier,k) + write(10000+mype, *) 'k muse ', k, muse(k), data(iuse,k) + write(10000+mype, *) 'k luse ', k, luse(k) do l=k+1,nobs - if(data(ilat,k) == data(ilat,l) .and. & - data(ilon,k) == data(ilon,l) .and. & + rstation_id = data(id,l) + ikx = nint(data(ikxx,l)) + itype = ictype(ikx) +!>>orig +! if(data(ilat,k) == data(ilat,l) .and. & +! data(ilon,k) == data(ilon,l) .and. & +! data(ier,k) < r1000 .and. data(ier,l) < r1000 .and. & +! muse(k) .and. muse(l))then +!< hPa --> Pa ikx=nint(data(ikxx,i)) itype=ictype(ikx) isubtype=icsubtype(ikx) @@ -380,7 +435,8 @@ subroutine setupps(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsa ! Link obs to diagnostics structure if (luse_obsdiag) then my_diag => obsdiagLList_nextNode(my_diagLL ,& - create = .not.lobsdiag_allocated ,& + + create = .not.lobsdiag_allocated ,& idv = is ,& iob = ioid(i) ,& ich = 1 ,& @@ -399,6 +455,11 @@ subroutine setupps(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsa obserror = max(cermin(ikx)*one_tenth,& min(cermax(ikx)*one_tenth,data(ier,i))) +!>>emily + obserr_step1 = data(ier,i)*r1000 +!<>emily + if (pob > 1.0e8 .or. dhgt > 1.0e8) & + write(6,*) 'emily checking pob dhgt ', pob, dhgt + write(6,*) 'emily checking missing value ', bmiss, bmiss*one_tenth + if (pob >= bmiss*one_tenth .or. dhgt >= bmiss) then + obserr_factor_prschk = 1.0_r_kind !emily_test + endif +!<>emily + rstation_id = data(id,i) + ikx = nint(data(ikxx,i)) + itype = ictype(ikx) + write(30000+mype,*) 'emily ................................................... ' + write(30000+mype,*) 'emily lat/lon D ', data(ilate,i), data(ilone,i) + write(30000+mype,*) 'emily lat/lon S ', slat(i), slon(i) + write(30000+mype,*) 'emily luse ', luse(i) + write(30000+mype,*) 'emily station_id ', station_id + write(30000+mype,*) 'emily obs typ ', itype + write(30000+mype,*) 'emily obs dtime ', data(itime,i) + write(30000+mype,*) 'emily obs ps ', data(ipres,i)*1000. + write(30000+mype,*) 'emily model elevation ', zsges + write(30000+mype,*) 'emily obs hgt ', dhgt + write(30000+mype,*) 'emily obs selev ', data(istnelv,i) + write(30000+mype,*) 'emily obs tv ', data(itemp,i) + write(30000+mype,*) 'emily ges tv ', tvges + write(30000+mype,*) 'emily zsges ', zsges + write(30000+mype,*) 'emily avg_tv0 ', tges0 + write(30000+mype,*) 'emily avg_tv ', tges + write(30000+mype,*) 'emily model_temp_sfc ', tvges + write(30000+mype,*) 'emily ps ', pgesorig*1000. + write(30000+mype,*) 'emily ps_cor ', pges*1000. + write(30000+mype,*) 'emily ....... ' + write(30000+mype,*) 'emily obs ps ', data(ipres,i)*1000. + write(30000+mype,*) 'emily tvflg ', data(iqt,i), iqtflg + write(30000+mype,*) 'emily model elevation ', zsges + write(30000+mype,*) 'emily ob elevation (ob hgt)', dhgt + write(30000+mype,*) 'emily obs_temp_sfc(obs tv) ', data(itemp,i) + write(30000+mype,*) 'emily model_temp_sfc ', tvges + write(30000+mype,*) 'emily tges2 ', tges2 + write(30000+mype,*) 'emily tges ', tges + write(30000+mype,*) 'emily orig oberr ', data(ier2,i)*r100 + write(30000+mype,*) 'emily current oberr ', data(ier,i)*r100 + write(30000+mype,*) 'emily rdelz ', rdelz + write(30000+mype,*) 'emily drbx ', drbx + write(30000+mype,*) 'emily drdp ', drdp*r1000 + write(30000+mype,*) 'emily obserr_factor_prschk ', obserr_factor_prschk +!< qcgross .or. ratio_errors < tiny_r_kind) then if (luse(i)) awork(6) = awork(6)+one error = zero ratio_errors = zero + gross_check_flag = 1 !emily else ratio_errors = ratio_errors/sqrt(dup(i)) end if @@ -650,6 +791,7 @@ subroutine setupps(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsa else err_final = huge_single endif + if (.not. muse(i)) err_final = huge_single !emily errinv_input = huge_single errinv_adjst = huge_single @@ -765,6 +907,61 @@ subroutine init_vars_ write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) endif +! get q ... + varname='q' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_q))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_q(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + ges_q(:,:,:,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) + ges_q(:,:,:,ifld)=rank3 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif + +! get u ... + varname='u' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_u))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_u(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + ges_u(:,:,:,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) + ges_u(:,:,:,ifld)=rank3 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif +! get v ... + varname='v' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_v))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_v(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + ges_v(:,:,:,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) + ges_v(:,:,:,ifld)=rank3 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif else write(6,*) trim(myname), ': inconsistent vector sizes (nfldsig,size(metguess_bundle) ',& nfldsig,size(gsi_metguess_bundle) @@ -882,6 +1079,7 @@ subroutine contents_binary_diag_(odiag) end subroutine contents_binary_diag_ subroutine contents_netcdf_diag_(odiag) + use obsmod, only: bmiss !emily type(obs_diag),pointer,intent(in):: odiag ! Observation class character(7),parameter :: obsclass = ' ps' @@ -891,12 +1089,20 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_metadata("Observation_Type", ictype(ikx) ) call nc_diag_metadata("Observation_Subtype", icsubtype(ikx) ) !Replace direct calls to nc_diag_metadata with the screening subroutine - call nc_diag_metadata_to_single("Latitude", data(ilate,i) ) - call nc_diag_metadata_to_single("Longitude", data(ilone,i) ) +! call nc_diag_metadata_to_single("Latitude", data(ilate,i) ) !orig +! call nc_diag_metadata_to_single("Longitude", data(ilone,i) ) !orig + call nc_diag_metadata("Latitude", slat(i) ) !emily + call nc_diag_metadata("Longitude", slon(i) ) !emily call nc_diag_metadata_to_single("Station_Elevation",data(istnelv,i) ) call nc_diag_metadata_to_single("Pressure", data(ipres,i),r10,'*' ) call nc_diag_metadata_to_single("Height", dhgt ) + call nc_diag_metadata_to_single("ObserrFactorDup", sqrt(dup(i)) ) !emily + call nc_diag_metadata_to_single("ObserrFactorPrsChk", obserr_factor_prschk) !emily + call nc_diag_metadata_to_single("GrossChkBound", grosschkbound ) !emily call nc_diag_metadata_to_single("Time", dtime,time_offset,'-' ) + call nc_diag_metadata_to_single("deltObsTime", data(itime,i) ) !emily_test + call nc_diag_metadata("GrossCheckFlag", sngl(gross_check_flag))!emily_test + call nc_diag_metadata_to_single("Setup_QC_Mark", data(iqt,i) ) !emily_test call nc_diag_metadata_to_single("Prep_QC_Mark", data(iqc,i) ) call nc_diag_metadata_to_single("Prep_Use_Flag", data(iuse,i) ) call nc_diag_metadata_to_single("Nonlinear_QC_Var_Jb",var_jb ) @@ -907,6 +1113,8 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_metadata("Analysis_Use_Flag", -1.0_r_single ) endif + call nc_diag_metadata_to_single("Error_Intermediate1", obserr_step1 ) !emily: Pa + call nc_diag_metadata_to_single("Error_Input", obserr0 ) !emily: Pa call nc_diag_metadata_to_single("Errinv_Input", errinv_input ) call nc_diag_metadata_to_single("Errinv_Adjust", errinv_adjst ) call nc_diag_metadata_to_single("Errinv_Final", errinv_final ) @@ -915,6 +1123,24 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_metadata_to_single("Obs_Minus_Forecast_adjusted",pob,pges,'-') call nc_diag_metadata_to_single("Obs_Minus_Forecast_unadjusted",pob,pgesorig,'-') +!>>emily + if (iqtflg) then ! this is virtual temperature + call nc_diag_metadata("Observation_Virtual_Temperature", sngl(dtemp) ) + call nc_diag_metadata("Observation_Sensible_Temperature", sngl(bmiss)) + else ! this is sensible temperature + call nc_diag_metadata("Observation_Virtual_Temperature", sngl(bmiss) ) + call nc_diag_metadata("Observation_Sensible_Temperature", sngl(dtemp)) + endif +!<>orig +! call nc_diag_metadata("Obs_Minus_Forecast_adjusted", sngl(pob-pges) ) +! call nc_diag_metadata("Obs_Minus_Forecast_unadjusted", sngl(pob-pgesorig)) +!<>emily + call nc_diag_metadata("Forecast_adjusted", sngl(pges) ) + call nc_diag_metadata("Forecast_unadjusted", sngl(pgesorig)) +!<>emily + do k = 1, nsig + kk = nsig-k+1 + utmp_reverse(kk) = utmp(k) + vtmp_reverse(kk) = vtmp(k) + ttmp_reverse(kk) = ttmp(k) + qtmp_reverse(kk) = qtmp(k) + hsges_reverse(kk) = hsges(k) + prsltmp2_reverse(kk) = prsltmp2(k) + enddo + do k = 1, nsig+1 + kk = (nsig+1)-k+1 + prsitmp_reverse(kk) = prsitmp(k) + enddo + call nc_diag_data2d("atmosphere_pressure_coordinate", sngl(prsltmp2_reverse*r1000)) + call nc_diag_data2d("atmosphere_pressure_coordinate_interface", sngl(prsitmp_reverse*r1000)) + call nc_diag_data2d("air_temperature", sngl(ttmp_reverse)) + call nc_diag_data2d("specific_humidity", sngl(qtmp_reverse)) + call nc_diag_data2d("eastward_wind", sngl(utmp_reverse)) + call nc_diag_data2d("northward_wind", sngl(vtmp_reverse)) + call nc_diag_data2d("geopotential_height", sngl(hsges_reverse) ) + !<>orig +! call nc_diag_data2d("atmosphere_pressure_coordinate", sngl(prsltmp2*r1000)) +! call nc_diag_data2d("atmosphere_pressure_coordinate_interface", sngl(prsitmp*r1000)) +! call nc_diag_data2d("air_temperature", sngl(ttmp)) +! call nc_diag_data2d("specific_humidity", sngl(qtmp)) +! call nc_diag_data2d("eastward_wind", sngl(utmp)) +! call nc_diag_data2d("northward_wind", sngl(vtmp)) +! call nc_diag_data2d("geopotential_height", sngl(hsges) ) +!<>emily + isli = data(idomsfc,i) + water_frac = 0.0 + land_frac = 0.0 + ice_frac = 0.0 + snow_frac = 0.0 + if (isli == 0) water_frac = 1.0 + if (isli == 1) land_frac = 1.0 + if (isli == 2) ice_frac = 1.0 + if (isli == 3) snow_frac = 1.0 + call nc_diag_metadata("Water_Fraction", sngl(water_frac)) + call nc_diag_metadata("Land_Fraction", sngl(land_frac)) + call nc_diag_metadata("Ice_Fraction", sngl(ice_frac)) + call nc_diag_metadata("Snow_Fraction", sngl(snow_frac)) +!<>emily + real(r_kind),dimension(nchanl):: varinv_grosschk, varinv_sdoei + real(r_kind),dimension(nchanl):: varinv_after_jsfcchk,varinv_after_sdoei + real(r_kind),dimension(nchanl):: varinv_after_grosschk + real(r_kind),dimension(nchanl):: varinv_after_wavenum,varinv_after_rangechk,varinv_after_topo,varinv_after_transmittop + real(r_kind),dimension(nchanl):: varinv_after_clddet,varinv_after_jsfcchk_land,varinv_after_nsstret + real(r_kind),dimension(nchanl):: varinv_after_grossroutinechk_over_ocean + real(r_kind),dimension(nchanl):: varinv_after_grossroutinechk + real(r_kind),dimension(nchanl):: varinv_after_sfcchk + real(r_kind),dimension(nchanl):: varinv_after_ch2chk + real(r_kind),dimension(nchanl):: varinv_after_scatteringchk + real(r_kind),dimension(nchanl):: varinv_after_sfcterrianchk + real(r_kind),dimension(nchanl):: varinv_after_clrfracchk + real(r_kind),dimension(nchanl):: varinv_after_stdchk + real(r_kind),dimension(nchanl):: varinv_after_stdadj + real(r_kind) :: pred9,pred10,pred11 !emily +!<>emily + do k = 1, nsig + kk = nsig-k+1 + utmp_reverse(kk) = utmp(k) + vtmp_reverse(kk) = vtmp(k) + ttmp_reverse(kk) = ttmp(k) + qtmp_reverse(kk) = qtmp(k) + hsges_reverse(kk) = hsges(k) + prsltmp2_reverse(kk) = prsltmp2(k) + enddo + do k = 1, nsig+1 + kk = (nsig+1)-k+1 + prsitmp_reverse(kk) = prsitmp(k) + enddo + call nc_diag_data2d("atmosphere_pressure_coordinate", sngl(prsltmp2_reverse*r1000)) + call nc_diag_data2d("atmosphere_pressure_coordinate_interface", sngl(prsitmp_reverse*r1000)) + call nc_diag_data2d("air_temperature", sngl(ttmp_reverse)) + call nc_diag_data2d("specific_humidity", sngl(qtmp_reverse)) + call nc_diag_data2d("eastward_wind", sngl(utmp_reverse)) + call nc_diag_data2d("northward_wind", sngl(vtmp_reverse)) + call nc_diag_data2d("geopotential_height", sngl(hsges_reverse) ) + !<>orig +! call nc_diag_data2d("atmosphere_pressure_coordinate", sngl(prsltmp2*r1000)) +! call nc_diag_data2d("atmosphere_pressure_coordinate_interface", sngl(prsitmp*r1000)) +! call nc_diag_data2d("air_temperature", sngl(ttmp)) +! call nc_diag_data2d("specific_humidity", sngl(qtmp)) +! call nc_diag_data2d("eastward_wind", sngl(utmp)) +! call nc_diag_data2d("northward_wind", sngl(vtmp)) +! call nc_diag_data2d("geopotential_height", sngl(hsges) ) +!<>emily + isli = data(idomsfc,i) + water_frac = 0.0 + land_frac = 0.0 + ice_frac = 0.0 + snow_frac = 0.0 + if (isli == 0) water_frac = 1.0 + if (isli == 1) land_frac = 1.0 + if (isli == 2) ice_frac = 1.0 + if (isli == 3) snow_frac = 1.0 + call nc_diag_metadata("Water_Fraction", sngl(water_frac)) + call nc_diag_metadata("Land_Fraction", sngl(land_frac)) + call nc_diag_metadata("Ice_Fraction", sngl(ice_frac)) + call nc_diag_metadata("Snow_Fraction", sngl(snow_frac)) +!<> JEDI + real(r_kind),dimension(nsig):: prsltmp2 + real(r_kind),dimension(nsig+1):: prsitmp + real(r_kind),dimension(nsig):: ttmp,qtmp,utmp,vtmp,hsges + ! GSI profiles are stored with bottom up index; output the profiles + ! with top down index + real(r_kind),dimension(nsig):: ttmp_reverse,tvtmp_reverse,qtmp_reverse,utmp_reverse,vtmp_reverse + real(r_kind),dimension(nsig):: hsges_reverse, zges_reverse,prsltmp2_reverse + real(r_kind),dimension(nsig):: zges_read_reverse, zges_geometric_reverse + real(r_kind),dimension(nsig+1):: prsitmp_reverse + real(r_kind) psges2 + !<< JEDI integer(i_kind) i,nchar,nreal,k,j,l,ii,itype,ijb ! Variables needed for new polar winds QC based on Log Normalized Vector Departure (LNVD) @@ -292,7 +309,8 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav integer(i_kind) ihgt,ier2,iuse,ilate,ilone integer(i_kind) izz,iprvd,isprvd integer(i_kind) idomsfc,isfcr,iskint,iff10 - integer(i_kind) ibb,ikk,ihil,idddd + integer(i_kind) kk,ibb,ikk,ihil,idddd,iamvq + integer(i_kind) water_frac,land_frac,ice_frac, snow_frac !emily integer(i_kind) num_bad_ikx,iprev_station @@ -333,6 +351,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav real(r_kind),allocatable,dimension(:,:,: ) :: ges_ps real(r_kind),allocatable,dimension(:,:,: ) :: ges_z + real(r_kind),allocatable,dimension(:,:,:,:) :: ges_q real(r_kind),allocatable,dimension(:,:,:,:) :: ges_u real(r_kind),allocatable,dimension(:,:,:,:) :: ges_v real(r_kind),allocatable,dimension(:,:,:,:) :: ges_tv @@ -501,6 +520,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav dlon=data(ilon,i) rstation_id = data(id,i) error=data(ier2,i) + obserr0=data(ier2,i) !emily var_jb=data(ijb,i) if(ikx < 1 .or. ikx > nconvtype) then num_bad_ikx=num_bad_ikx+1 @@ -577,6 +597,53 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav call tintrp2a1(ges_lnprsl,prsltmp,dlat,dlon,dtime,hrdifsig,& nsig,mype,nfldsig) +!>>JEDI +! GEOVALS for UFO eval + psges2 = psges ! keep in cb + prsltmp2 = exp(prsltmp) + call tintrp2a1(ges_prsi,prsitmp,dlat,dlon,dtime,hrdifsig,& + nsig+1,mype,nfldsig) + call tintrp2a1(ges_tsen,ttmp,dlat,dlon,dtime,hrdifsig,& + nsig,mype,nfldsig) + call tintrp2a1(ges_q,qtmp,dlat,dlon,dtime,hrdifsig,& + nsig,mype,nfldsig) + call tintrp2a1(ges_u,utmp,dlat,dlon,dtime,hrdifsig,& + nsig,mype,nfldsig) + call tintrp2a1(ges_v,vtmp,dlat,dlon,dtime,hrdifsig,& + nsig,mype,nfldsig) + call tintrp2a1(geop_hgtl,hsges,dlat,dlon,dtime,hrdifsig,& !orig + nsig,mype,nfldsig) + ! geopotential height at obs location (lat/lon) and time + call tintrp2a1(geop_hgtl,zges_read,dlat,dlon,dtime,hrdifsig, & + nsig,mype,nfldsig) + ! model virtual temperature (ges_tv) at obs location (lat/lon) + ! obs time, and model surface (lower model level) + call tintrp31(ges_tv,sfctges,dlat,dlon,log(psges),dtime, & + hrdifsig,mype,nfldsig) + +! Convert geopotential height at layer midpoints to geometric +! height using equations (17, 20, 23) in MJ Mahoney's note +! "A discussion of various measures of altitude" (2001). +! Available on the web at +! http://mtp.jpl.nasa.gov/notes/altitude/altitude.html +! +! termg = equation 17 +! termr = equation 21 +! termrg = first term in the denominator of equation 23 +! zges = equation 23 + + slat = data(ilate,i)*deg2rad + sin2 = sin(slat)*sin(slat) + termg = grav_equator * & + ((one+somigliana*sin2)/sqrt(one-eccentricity*eccentricity*sin2)) + termr = semi_major_axis /(one + flattening + grav_ratio - & + two*flattening*sin2) + termrg = (termg/grav)*termr + do k=1,nsig + zges_geometric(k) = (termr*zges_read(k)) / (termrg-zges_read(k)) ! eq (23) + end do +! END GEOVALS +!<>orig +!! Get guess surface elevation and geopotential height profile +!! at observation location. +! call tintrp2a1(geop_hgtl,zges,dlat,dlon,dtime,hrdifsig,& +! nsig,mype,nfldsig) +!<=223 .and. itype<=228) .or. itype == 218 .or. sfc_data) .and. .not.twodvar_regional) then -! Convert geopotential height at layer midpoints to geometric -! height using equations (17, 20, 23) in MJ Mahoney's note -! "A discussion of various measures of altitude" (2001). -! Available on the web at -! http://mtp.jpl.nasa.gov/notes/altitude/altitude.html +!>>orig +!! Convert geopotential height at layer midpoints to geometric +!! height using equations (17, 20, 23) in MJ Mahoney's note +!! "A discussion of various measures of altitude" (2001). +!! Available on the web at +!! http://mtp.jpl.nasa.gov/notes/altitude/altitude.html +!! +!! termg = equation 17 +!! termr = equation 21 +!! termrg = first term in the denominator of equation 23 +!! zges = equation 23 ! -! termg = equation 17 -! termr = equation 21 -! termrg = first term in the denominator of equation 23 -! zges = equation 23 - - slat = data(ilate,i)*deg2rad - sin2 = sin(slat)*sin(slat) - termg = grav_equator * & - ((one+somigliana*sin2)/sqrt(one-eccentricity*eccentricity*sin2)) - termr = semi_major_axis /(one + flattening + grav_ratio - & - two*flattening*sin2) - termrg = (termg/grav)*termr - do k=1,nsig - zges(k) = (termr*zges(k)) / (termrg-zges(k)) ! eq (23) - end do +! slat = data(ilate,i)*deg2rad +! sin2 = sin(slat)*sin(slat) +! termg = grav_equator * & +! ((one+somigliana*sin2)/sqrt(one-eccentricity*eccentricity*sin2)) +! termr = semi_major_axis /(one + flattening + grav_ratio - & +! two*flattening*sin2) +! termrg = (termg/grav)*termr +! do k=1,nsig +! zges(k) = (termr*zges(k)) / (termrg-zges(k)) ! eq (23) +! zges_geometric(k) = (termr*zges_read(k)) / (termrg-zges_read(k)) ! eq (23) !emily_gmao +! end do +!<>JEDI + ! GEOVALS + !>>emily + do k = 1, nsig + kk = nsig-k+1 + utmp_reverse(kk) = utmp(k) + vtmp_reverse(kk) = vtmp(k) + ttmp_reverse(kk) = ttmp(k) + tvtmp_reverse(kk) = tges(k) !emily + qtmp_reverse(kk) = qtmp(k) + hsges_reverse(kk) = hsges(k) + zges_read_reverse(kk)= zges_read(k) + zges_geometric_reverse(kk)= zges_geometric(k) + zges_reverse(kk) = zges(k) + prsltmp2_reverse(kk) = prsltmp2(k) + enddo + do k = 1, nsig+1 + kk = (nsig+1)-k+1 + prsitmp_reverse(kk) = prsitmp(k) + enddo + call nc_diag_data2d("atmosphere_pressure_coordinate", sngl(prsltmp2_reverse*r1000)) + call nc_diag_data2d("atmosphere_pressure_coordinate_interface", sngl(prsitmp_reverse*r1000)) + call nc_diag_data2d("air_temperature", sngl(ttmp_reverse)) + call nc_diag_data2d("virtual_temperature", sngl(tvtmp_reverse)) !emily + call nc_diag_data2d("specific_humidity", sngl(qtmp_reverse)) + call nc_diag_data2d("eastward_wind", sngl(utmp_reverse)) + call nc_diag_data2d("northward_wind", sngl(vtmp_reverse)) +! call nc_diag_data2d("geopotential_height", sngl(hsges_reverse) ) !orig + call nc_diag_data2d("geopotential_height", sngl(zges_read_reverse) ) !emily + call nc_diag_data2d("geometric_height", sngl(zges_geometric_reverse) ) !emily + !<>orig +! call nc_diag_data2d("atmosphere_pressure_coordinate", sngl(prsltmp2*r1000)) +! call nc_diag_data2d("atmosphere_pressure_coordinate_interface", sngl(prsitmp*r1000)) +! call nc_diag_data2d("air_temperature", sngl(ttmp)) +! call nc_diag_data2d("specific_humidity", sngl(qtmp)) +! call nc_diag_data2d("eastward_wind", sngl(utmp)) +! call nc_diag_data2d("northward_wind", sngl(vtmp)) +! call nc_diag_data2d("geopotential_height", sngl(hsges) ) +! call nc_diag_data2d("model_height", sngl(zges) ) +!<>emily + water_frac = 0.0 + land_frac = 0.0 + ice_frac = 0.0 + snow_frac = 0.0 + if (isli == 0) water_frac = 1.0 + if (isli == 1) land_frac = 1.0 + if (isli == 2) ice_frac = 1.0 + if (isli == 3) snow_frac = 1.0 + call nc_diag_metadata("Water_Fraction", sngl(water_frac)) + call nc_diag_metadata("Land_Fraction", sngl(land_frac)) + call nc_diag_metadata("Ice_Fraction", sngl(ice_frac)) + call nc_diag_metadata("Snow_Fraction", sngl(snow_frac)) + !<= r40 ) then itrop_k = itrp_pv elseif (slatd(i) >= r20) then @@ -197,6 +198,7 @@ subroutine tpause(mype,method) trop_pv(i,j) = prs(itrp_pv)*r0_01 !hPa trop_oz(i,j) = prs(itrp_oz)*r0_01 !hPa trop_pvoz(i,j) = prs(itrop_k)*r0_01 !hPa + end do end do