diff --git a/src/gsi/genstats_gps.f90 b/src/gsi/genstats_gps.f90 index acf5ca275..2c83c1686 100644 --- a/src/gsi/genstats_gps.f90 +++ b/src/gsi/genstats_gps.f90 @@ -228,6 +228,7 @@ subroutine genstats_gps(bwork,awork,toss_gps_sub,conv_diagsave,mype) ! new module interface name. gpsStats_destroy(). ! 2016-11-29 shlyaeva - increase the size of nreal for saving linearized Hx for EnKF ! 2016-12-09 mccarty - add ncdiag writing support +! 2024-12-17 li - add new hybrid obs error model by Chris Riedel ! ! input argument list: ! toss_gps_sub - array of qc'd profile heights @@ -247,7 +248,7 @@ subroutine genstats_gps(bwork,awork,toss_gps_sub,conv_diagsave,mype) use m_obsdiagNode, only: obs_diag use m_obsdiagNode, only: obsdiagNode_set use obsmod, only: nprof_gps,lobsdiag_forenkf - use obsmod, only: lobsdiagsave,luse_obsdiag + use obsmod, only: lobsdiagsave,luse_obsdiag,nobs_gps use obsmod, only: binary_diag,netcdf_diag,dirname,ianldate use nc_diag_write_mod, only: nc_diag_init, nc_diag_header, nc_diag_metadata, & nc_diag_write, nc_diag_data2d, nc_diag_metadata_to_single @@ -255,10 +256,11 @@ subroutine genstats_gps(bwork,awork,toss_gps_sub,conv_diagsave,mype) use gridmod, only: nsig,regional use constants, only: tiny_r_kind,half,wgtlim,one,two,zero,five,four use qcmod, only: npres_print,ptop,pbot - use mpimod, only: ierror,mpi_comm_world,mpi_rtype,mpi_sum,mpi_max + use mpimod, only: ierror,mpi_comm_world,mpi_rtype,mpi_sum,mpi_max,mpi_integer,npe use jfunc, only: jiter,miter,jiterstart use gsi_4dvar, only: nobs_bins use convinfo, only: nconvtype + use m_gpsrhs, only: rdiagbuf implicit none ! Declare passed variables @@ -287,12 +289,33 @@ subroutine genstats_gps(bwork,awork,toss_gps_sub,conv_diagsave,mype) real(r_kind),dimension(max(1,nprof_gps)):: high_gps,high_gps_sub real(r_kind),dimension(max(1,nprof_gps)):: dobs_height,dobs_height_sub + real(r_kind),dimension(nobs_gps) :: complete_hght,complete_gps,complete_qc, & + complete_prof + real(r_single),allocatable,dimension(:,:)::sdiag character(8),allocatable,dimension(:):: cdiag + integer(i_kind) :: iprof,counts,total_obs + real(r_kind) :: gpsOb,gpsHght + !integer(i_kind), dimension(max(1,nprof_gps)) :: hold_profs + !integer(i_kind),allocatable, dimension(:) :: unique_IDs + !real(r_kind),dimension(nobs_gps) :: array_hght,array_gps,array_qc, & + ! array_prof + real(r_kind),allocatable,dimension(:) :: array_hght,array_gps + real(r_kind),allocatable,dimension(:) :: array_qc,array_prof + real(r_kind),allocatable,dimension(:) :: collect_hght,collect_gps,collect_qc, & + collect_prof + real(r_kind),dimension(max(1,nobs_gps)) :: holder_hght,holder_gps,holder_qc, & + holder_prof,STD4060 + real(r_kind),allocatable,dimension(:) :: profile_benda,profile_impact, & + sorted_profile_benda,sorted_profile_impact + real(r_kind) :: old_err,new_err,input_impact,input_std4060,input_FracLsw, & + relative_error + character(len=20) :: hold_place + !!!!!!!!!!!!!!!!!!!!! type(obs_diag), pointer :: obsptr => NULL() - - integer(i_kind) :: nnz, nind + integer(i_kind),allocatable,dimension(:) :: revcounts,displs,indices + integer(i_kind) :: nnz, nind,nobs,total_size,ind,prof_size type(gps_ob_type), pointer:: gpsptr type(gps_all_ob_type), pointer:: gps_allptr @@ -304,7 +327,139 @@ subroutine genstats_gps(bwork,awork,toss_gps_sub,conv_diagsave,mype) if (mype==0) write(6,*)'GENSTATS_GPS: no profiles to process (nprof_gfs=',nprof_gps,'), EXIT routine' return endif - +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + call mpi_comm_size(mpi_comm_world,total_size,ierror) + nobs = 0 + ind = 1 + holder_hght(:) = 0.0_r_kind;holder_gps(:) = 0.0_r_kind + holder_qc(:) = 0.0_r_kind;holder_prof(:) = 0.0_r_kind + DO ii=1,nobs_bins + gps_allptr => gps_allhead(ii)%head + do while (associated(gps_allptr)) + luse = gps_allptr%luse + if (luse) then + holder_hght(ind) = gps_allptr%rdiag(7) + holder_gps(ind) = gps_allptr%rdiag(17) + holder_qc(ind) = gps_allptr%rdiag(10) + holder_prof(ind) = gps_allptr%rdiag(2) + nobs = nobs + 1 + ind = ind + 1 + endif + gps_allptr => gps_allptr%llpoint + end do + END DO + allocate(collect_hght(nobs), source=holder_hght(1:ind)) + allocate(collect_gps(nobs), source=holder_gps(1:ind)) + allocate(collect_qc(nobs), source=holder_qc(1:ind)) + allocate(collect_prof(nobs), source=holder_prof(1:ind)) + + if (mype == 0) allocate(revcounts(total_size),array_hght(nobs_gps),array_gps(nobs_gps), & + array_qc(nobs_gps),array_prof(nobs_gps)) + call mpi_gather(nobs,1,mpi_integer, & + revcounts,1,mpi_integer,0,mpi_comm_world,ierror) + if (mype == 0) then + allocate(displs(total_size)) + displs(1) = 0 + do ii=2,total_size + displs(ii) = displs(ii-1) + revcounts(ii-1) + enddo + !write(6,*) "CHECK DISPLS: ",displs + array_hght(:) = 0.0_r_kind + array_gps(:) = 0.0_r_kind + array_qc(:) = 0.0_r_kind + array_prof(:) = 0.0_r_kind + endif + call mpi_gatherv(collect_hght,nobs,mpi_rtype,& + array_hght,revcounts,displs,mpi_rtype, & + 0,mpi_comm_world,ierror) + call mpi_gatherv(collect_gps,nobs,mpi_rtype,& + array_gps,revcounts,displs,mpi_rtype, & + 0,mpi_comm_world,ierror) + call mpi_gatherv(collect_qc,nobs,mpi_rtype,& + array_qc,revcounts,displs,mpi_rtype, & + 0,mpi_comm_world,ierror) + call mpi_gatherv(collect_prof,nobs,mpi_rtype,& + array_prof,revcounts,displs,mpi_rtype, & + 0,mpi_comm_world,ierror) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!! Compute STD4060 for hybrid error model + if (mype == 0) then + STD4060(:) = 0.0_r_kind + do ii=1,nprof_gps + prof_size = count(((array_prof(:) == ii).and.(array_hght(:)>=40.E3).and.(array_hght(:)<=60.E3).and.(array_qc/=6.0_r_kind))) + if (prof_size <= 15) then + STD4060(ii) = 40.0_r_kind + cycle + endif + allocate(profile_benda(prof_size),profile_impact(prof_size),indices(prof_size),sorted_profile_benda(prof_size),& + sorted_profile_impact(prof_size)) + profile_benda(:) = pack(array_gps,((array_prof == ii).and.(array_hght(:)>=40.E3).and.(array_hght(:)<=60.E3).and.(array_qc/=6.0_r_kind))) + profile_impact(:) = pack(array_hght,((array_prof == ii).and.(array_hght>=40.E3).and.(array_hght<=60.E3).and.(array_qc/=6.0_r_kind))) + call sort(prof_size,profile_impact,indices) + do k=1,prof_size + sorted_profile_benda(k) = profile_benda(indices(k)) + sorted_profile_impact(k) = profile_impact(indices(k)) + enddo + call cal_std4060(prof_size,sorted_profile_benda,sorted_profile_impact,STD4060(ii)) + if (STD4060(ii) > 40.0_r_kind) then + STD4060(ii) = 40.0_r_kind + endif + deallocate(profile_benda,profile_impact,indices,sorted_profile_benda, & + sorted_profile_impact) + enddo + endif + call mpi_bcast(STD4060,nprof_gps,mpi_rtype,0,mpi_comm_world,ierror) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + deallocate(collect_hght,collect_gps,collect_qc,collect_prof) + if (mype == 0) deallocate(revcounts,displs,array_hght,array_gps, & + array_qc,array_prof) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! compute error + DO ii=1,nobs_bins + gps_allptr => gps_allhead(ii)%head + do while (associated(gps_allptr)) + muse = gps_allptr%muse + if (muse) then + if (gps_allptr%rdiag(25) == 2.0) then + input_impact = gps_allptr%rdiag(7) + call compute_threeCH_uncert(input_impact*r1em3,relative_error) + new_err = (relative_error/100._r_kind) * gps_allptr%rdiag(17) + if (new_err < 1.0D-6) new_err = 1.0D-6 + input_std4060 = STD4060(int(gps_allptr%rdiag(2))) + gps_allptr%rdiag(24) = input_std4060 + old_err = (1.0_r_kind/gps_allptr%rdiag(14)) + gps_allptr%ratio_err = old_err/new_err + gpsptr => gps_allptr%mmpoint + if (associated(gpsptr)) then + gpsptr%raterr2 = gps_allptr%ratio_err **2 + endif + gps_allptr%rdiag(15) = gps_allptr%ratio_err * gps_allptr%rdiag(14) + gps_allptr%rdiag(16) = one/new_err + else + input_impact = gps_allptr%rdiag(7) + input_FracLsw = gps_allptr%rdiag(23) + if (input_FracLsw > 40.0_r_kind) then + input_FracLsw = 40.0_r_kind + endif + input_std4060 = STD4060(int(gps_allptr%rdiag(2))) + call error_model(input_impact,input_FracLsw,input_std4060,relative_error) + new_err = (relative_error/100._r_kind) * gps_allptr%rdiag(17) + if (new_err < 1.0D-6) new_err = 1.0D-6 + gps_allptr%rdiag(24) = input_std4060 + old_err = (1.0_r_kind/gps_allptr%rdiag(14)) + gps_allptr%ratio_err = old_err/new_err + gpsptr => gps_allptr%mmpoint + if (associated(gpsptr)) then + gpsptr%raterr2 = gps_allptr%ratio_err **2 + endif + gps_allptr%rdiag(15) = gps_allptr%ratio_err * gps_allptr%rdiag(14) + gps_allptr%rdiag(16) = one/new_err + endif + endif + gps_allptr => gps_allptr%llpoint + end do + END DO +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Reduce sub-domain specific QC'd profile height cutoff values to ! maximum global value for each profile toss_gps=zero @@ -427,7 +582,7 @@ subroutine genstats_gps(bwork,awork,toss_gps_sub,conv_diagsave,mype) end do END DO if(icnt > 0)then - nreal =22 + nreal =25 ioff =nreal if (lobsdiagsave) nreal=nreal+4*miter+1 if (save_jacobian) then @@ -708,6 +863,287 @@ subroutine genstats_gps(bwork,awork,toss_gps_sub,conv_diagsave,mype) ! Destroy arrays holding gps data call destroy_genstats_gps contains +subroutine sort(num,x,indx) +integer(i_kind), intent(in) :: num +real(r_kind), intent(in) :: x(num) +integer(i_kind), intent(inout) :: indx(num) + +integer(i_kind) :: ind, i, j, value_indx, line +real(r_kind) :: values + +! Initialize the index array to input order +do i = 1, num + indx(i) = i +end do + +! Only one element, just send it back +if(num <= 1) return + +line = num / 2 + 1 +ind = num +do + if(line > 1) then + line = line - 1 + values = x(indx(line)) + value_indx = indx(line) + else + values = x(indx(ind)) + value_indx = indx(ind) + + indx(ind) = indx(1) + ind = ind - 1 + if(ind == 1) then + indx(1) = value_indx + return + endif + endif + + i = line + j = 2 * line + + do while(j <= ind) + if(j < ind) then + if(x(indx(j)) < x(indx(j + 1))) j = j + 1 + endif + if(values < x(indx(j))) then + indx(i) = indx(j) + i = j + j = 2 * j + else + j = ind + 1 + endif + end do + + indx(i) = value_indx + +end do +end subroutine sort +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +subroutine lsqs(n,x,y,m,c) + ! Subroutine computes a least squares fit + ! Returns the slope and intercept + integer(i_kind) ,intent(in) :: n + real(r_kind) ,intent(in) :: x(n) + real(r_kind) ,intent(in) :: y(n) + real(r_kind) ,intent(out) :: m,c + !!!! + real(r_kind) :: top,bot + !real(r_kind) :: bot + integer(i_kind) :: ii + !!!! + real(r_kind) x_avg,y_avg + !!!! + x_avg = sum(x)/float(n) + y_avg = sum(y)/float(n) + !Sum of terms + top = 0.0 + bot = 0.0 + do ii=1,n + top = top + ((x(ii) - x_avg)*y(ii)) + bot = bot + ((x(ii) - x_avg)**2) + enddo + ! Calculate Slope + m = top/bot + ! Calculate intercept + c = y_avg - (m*x_avg) + return +end subroutine lsqs +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +subroutine cal_std4060(num_pts,benda,impacts,stddev) + integer(i_kind) ,intent(in) :: num_pts + real(r_kind) ,intent(in) :: benda(num_pts) + real(r_kind) ,intent(in) :: impacts(num_pts) + real(r_kind) ,intent(out) :: stddev + ! + real(r_kind) least_squares_fit(num_pts),norm_fit(num_pts) + real(r_kind) slope,intercept,mean + ! + call lsqs(num_pts, impacts, log(benda),slope,intercept) + ! + least_squares_fit = exp(slope*impacts + intercept) + norm_fit = (least_squares_fit - benda)/least_squares_fit * 100.0_r_kind + ! + mean = sum(norm_fit)/float(num_pts) + stddev = sqrt(sum((norm_fit - mean)**2)/(float(num_pts))) + return +end subroutine cal_std4060 +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +subroutine std4060_errmodel(input_hght,input_std4060,error) + real(r_kind) ,intent(in) :: input_hght,input_std4060 + real(r_kind) ,intent(out) :: error + ! + integer(i_kind) :: n_coefs_x = 4 + integer(i_kind) :: n_coefs_y = 5 + integer(i_kind) :: aa,bb + real(r_kind), dimension(4,5) :: coefs + ! + coefs(:,:) = 0.0_r_kind + coefs(1,:) = (/ 1.25000000000000000000000000D+00,0.00000000000000000000000000D+00,& + 0.00000000000000000000000000D+00,0.00000000000000000000000000D+00,& + 0.00000000000000000000000000D+00 /) + coefs(2,:) = (/ 1.99128045049602960634847933D-06, 6.75863476049698542422747641D-06,& + 8.25883030063171995856295508D-07,-2.74722471254849586614206177D-08,& + 1.82378331607358961583961066D-10 /) + coefs(3,:) = (/ -6.55621669612269976911399984D-09,3.15300096566169468056300701D-09,& + -3.77925688851882344702164472D-10,1.06712353466978524015944812D-11,& + -8.65642894028274826191107202D-14 /) + coefs(4,:) = (/ 2.66873781200705724938244662D-13,-7.28569101066463132754885337D-14,& + 1.68201368786068197554617494D-14,-5.49152377448293150781835625D-16,& + 5.43287924985643301111631898D-18 /) + !!!! + error = 0.0_r_kind + do aa=1, n_coefs_x + do bb=1,n_coefs_y + error = error + (coefs(aa,bb)*(input_hght**(aa-1))*(input_std4060**(bb-1))) + enddo + enddo + return +end subroutine std4060_errmodel +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +subroutine lsw_errmodel(input_hght,input_lsw,error) + real(r_kind) ,intent(in) :: input_hght,input_lsw + real(r_kind) ,intent(out) :: error + ! + integer(i_kind) :: n_coefs_x = 4 + integer(i_kind) :: n_coefs_y = 4 + integer(i_kind) :: aa,bb,a,b + real(r_kind), dimension(4,4) :: coefs + ! + coefs(:,:) = 0.0_r_kind + coefs(1,:) = (/ 1.25000000000000000000000000D+00,0.00000000000000000000000000D+00,& + 0.00000000000000000000000000D+00,0.00000000000000000000000000D+00 /) + coefs(2,:) = (/ 8.79644062400577173878313264D-04,5.00524082314834873816411509D-04,& + -7.06272962581374830816913560D-06,2.41695480315650456302628568D-08 /) + coefs(3,:) = (/-1.32890117750533099169737636D-07,-1.00068031233479483283742046D-07,& + 8.15492853573235497611204411D-10,1.37351779555849664534897326D-12 /) + coefs(4,:) = (/ 6.84487441321994421597774228D-12,6.99590380726247359733035140D-12,& + -3.31568154481349481116764491D-14,-2.52493347933840834626797529D-16 /) + + !!!! + error = 0.0_r_kind + do aa=1, n_coefs_x + do bb=1,n_coefs_y + error = error + (coefs(aa,bb)*(input_hght**(aa-1))*(input_lsw**(bb-1))) + enddo + enddo + return +end subroutine lsw_errmodel +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +subroutine error_model(impact_hght,lsw,std4060,errstd) + real(r_kind) ,intent(in) :: impact_hght,lsw,std4060 + real(r_kind) ,intent(inout) :: errstd + ! + real(r_kind),parameter :: min_error = 1.25_r_kind + real(r_kind),parameter :: lsw_boundary = 10.E3 + real(r_kind),parameter :: std4060_boundary = 30.E3 + !!!! + if (impact_hght <= lsw_boundary) then + call lsw_errmodel(lsw_boundary - impact_hght,lsw,errstd) + else if (impact_hght >= std4060_boundary) then + call std4060_errmodel(impact_hght-std4060_boundary,std4060,errstd) + else + errstd = min_error + endif + return +end subroutine error_model +!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! Subroutine to compute statistical error (from 3CHMethod) for a given +!! height using a spline fit +subroutine compute_threeCH_uncert(input_impact,error) + use kinds, only: r_kind,i_kind,r_single + real(r_kind), intent(in) :: input_impact + real(r_kind), intent(out) :: error + !!Define coefs arrays + real(r_kind),dimension(66) :: a,b,c,d + real(r_single),dimension(67) :: interp_impact + integer(i_kind) :: l,level_index + real(r_kind) :: interp_err,min_error + !!!!!!!!!!!!! + a(:) = 0.0_r_kind + a(:) = (/ -1.8072472782194637_r_kind,-4.553293905977256_r_kind,-1.042632730866977_r_kind,9.100931165806573_r_kind,1.6912359147188525_r_kind,-1.2982719737865533_r_kind,0.9843419625507028_r_kind,-0.7317115785860473_r_kind,& + -0.4110048129739834_r_kind,1.183432121346243_r_kind,0.07889575704898277_r_kind,0.1037405051596565_r_kind,0.12554091733179007_r_kind,-0.13836375493037956_r_kind,-0.31793666720456226_r_kind,1.2790616977422875_r_kind,& + 0.009323454981171042_r_kind,-0.03068971657715175_r_kind,-0.012434437302992873_r_kind,0.06721442044491856_r_kind,3.7505350493542786D-06,-0.0334244677836707_r_kind,-0.02819257166121869_r_kind,-0.2493119206629355_r_kind,& + -0.016830473861980624_r_kind,0.4977689706039675_r_kind,0.17375087190815963_r_kind,-0.01710203102076152_r_kind,-0.00047596910806305126_r_kind,0.005636654598357445_r_kind,0.0004980015635230243_r_kind,-0.006028788035461538_r_kind,& + -0.007278154670571774_r_kind,0.005612112069916075_r_kind,-0.003520194640361976_r_kind,-0.010140903000390833_r_kind,0.012442340314416822_r_kind,-0.013974982052712498_r_kind,-0.006212583721363296_r_kind,0.009663352106410006_r_kind,& + 0.0006611375699540623_r_kind,-0.021742887710101355_r_kind,0.01441853441306784_r_kind,-0.007862465965627119_r_kind,0.017098553464325428_r_kind,0.011225188873242864_r_kind,-0.052329516422997324_r_kind,0.02090873344068611_r_kind,& + -0.026410224720464726_r_kind,0.01045825824505231_r_kind,0.03388155526222891_r_kind,0.007517729045328103_r_kind,-0.05995059920209567_r_kind,-0.0418887916237658_r_kind,0.0779203213611881_r_kind,-0.06077763068895248_r_kind,& + 0.07575414794899071_r_kind,-0.014759941559792988_r_kind,-0.01213449120653376_r_kind,0.020914920476412746_r_kind,-0.04953774742002803_r_kind,0.25641102955203987_r_kind,-0.23665477425886117_r_kind,-0.03603425208201827_r_kind,& + 0.14615239845937822_r_kind,-0.06624648006639511_r_kind /) + + b(:) = 0.0_r_kind + b(:) = (/ -2.0561518447819296_r_kind,0.22049510820670193_r_kind,-0.4190389657895395_r_kind,-9.198242153793373_r_kind,-0.9470909055722503_r_kind,1.1951837100964111_r_kind,-0.6014140841997015_r_kind,0.6781444896434916_r_kind,& + 0.15003497758938922_r_kind,-0.8993086265952344_r_kind,0.20466466226031566_r_kind,0.18582750259752912_r_kind,0.170732464302648_r_kind,0.29451752503827366_r_kind,0.28822391802870695_r_kind,-0.6705595296494151_r_kind,& + 0.13353892042182558_r_kind,0.07289920169604404_r_kind,0.0208726883903759_r_kind,-0.07184111814190572_r_kind,0.0003379043049377586_r_kind,0.09345868741198415_r_kind,0.09717121911947513_r_kind,0.31906967975033607_r_kind,& + -0.018060435327087093_r_kind,-0.9433878688650862_r_kind,-0.13252807790515786_r_kind,0.06912762861774696_r_kind,0.015100757244949484_r_kind,0.0004800603807383688_r_kind,0.0010927527090346537_r_kind,0.015791494516346934_r_kind,& + 0.01563726374305644_r_kind,-0.00441866911778142_r_kind,0.01567255621904334_r_kind,0.023098270970024787_r_kind,-0.00936941457649551_r_kind,0.038235607387624285_r_kind,0.018612906103360458_r_kind,-0.0028718956893171443_r_kind,& + 0.01997358604096397_r_kind,0.04864454552748257_r_kind,-0.009024393794918734_r_kind,0.031210993807183396_r_kind,0.000167081939952074_r_kind,0.031472883074656566_r_kind,0.12331185037652542_r_kind,-0.0008671938265990553_r_kind,& + 0.07384170712382765_r_kind,0.011975287250526767_r_kind,0.0023017716173047154_r_kind,0.07679146481147481_r_kind,0.1722967313253665_r_kind,0.09943744879226468_r_kind,-0.06185583905726744_r_kind,0.16991926963328696_r_kind,& + -0.024082417185988136_r_kind,0.16481175816820093_r_kind,0.16889848541950192_r_kind,0.14476259759272248_r_kind,0.26699662506041877_r_kind,-0.06743664171335473_r_kind,0.8258835018193778_r_kind,0.4525933320273996_r_kind,& + 0.2969649654467448_r_kind,0.7880092025646874_r_kind /) + + c(:) = 0.0_r_kind + c(:) = (/ 7.80707613866254_r_kind,4.395488835216012_r_kind,1.2010135139397722_r_kind,0.0_r_kind,-2.3725437794384443_r_kind,-2.051207748971555_r_kind,-1.8297280192150587_r_kind,-1.692885631501733_r_kind,-1.5635248257977767_r_kind,& + -1.721743457938875_r_kind,-1.733477993524427_r_kind,-1.4696415134773744_r_kind,-1.206008632010103_r_kind,-0.9411204797086128_r_kind,-0.7503757708681239_r_kind,-0.7006043532428387_r_kind,-0.41186760958553814_r_kind,& + -0.11681940379837381_r_kind,-0.06309015013774097_r_kind,-0.0586480852659678_r_kind,-0.0006870602150235801_r_kind,0.0_r_kind,0.08664397147295622_r_kind,0.19640869472825037_r_kind,0.08661229224011606_r_kind,0.0_r_kind,& + -0.3934688259182698_r_kind,-0.13727236600410664_r_kind,-0.050323201830897266_r_kind,-0.021549594665187452_r_kind,-0.0036795101086383804_r_kind,0.0_r_kind,0.013496624926309255_r_kind,0.022936688400706822_r_kind,0.030935686374892208_r_kind,& + 0.05172021489189296_r_kind,0.06749404783077004_r_kind,0.08608223962102947_r_kind,0.12062850823814054_r_kind,0.13921656928077156_r_kind,0.16246283422136729_r_kind,0.2043934190131574_r_kind,0.23645384693781848_r_kind,& + 0.26166066258718457_r_kind,0.30049525230467_r_kind,0.35212507657755043_r_kind,0.4487464093465922_r_kind,0.5383815608306511_r_kind,0.5993733734995114_r_kind,0.6678261135857724_r_kind,0.723151462821983_r_kind,0.8293996718432791_r_kind,& + 1.0055357886022132_r_kind,1.1702774536466594_r_kind,1.243485976359891_r_kind,1.3535352623289207_r_kind,1.511040909528637_r_kind,1.690138519003633_r_kind,1.975482210660656_r_kind,2.2768757078800586_r_kind,2.6291456644947417_r_kind,& + 3.014525672355495_r_kind,3.6488854775849053_r_kind,4.590688158447078_r_kind,5.387772066255823_r_kind,6.420159192527448_r_kind /) + d(:) = 0.0_r_kind + d(:) = (/ 6.6174092389358306_r_kind,9.781003437294185_r_kind,11.46470989370671_r_kind,11.830127817870839_r_kind,10.668183675148317_r_kind,9.456543548375889_r_kind,8.567451604690895_r_kind,7.625276819352278_r_kind,6.856906178689028_r_kind,& + 6.061276908565739_r_kind,5.123507038115774_r_kind,4.317796176549762_r_kind,3.642399858605414_r_kind,3.0977712733424982_r_kind,2.683544945381463_r_kind,2.3406709560540073_r_kind,1.9826116092380202_r_kind,1.7136063750554786_r_kind,& + 1.6389964563759971_r_kind,1.5843445573256392_r_kind,1.5210697743626842_r_kind,1.5207243689876477_r_kind,1.5807585886159612_r_kind,1.7363812075471738_r_kind,2.002547661362825_r_kind,2.054269044413873_r_kind,1.6086501461527545_r_kind,& + 1.2564041142374864_r_kind,1.1711573458303652_r_kind,1.1354589321363544_r_kind,1.1200260524502628_r_kind,1.117937296614182_r_kind,1.1277000030950675_r_kind,1.1495557370938614_r_kind,1.1736858684467029_r_kind,1.2167739164002764_r_kind,& + 1.2814514992618034_r_kind,1.3520184728304947_r_kind,1.462361337786436_r_kind,1.5953901684065737_r_kind,1.741398194104438_r_kind,1.9244957519367234_r_kind,2.155790828767262_r_kind,2.3976388163232296_r_kind,2.6826480067519705_r_kind,& + 3.000408894460918_r_kind,3.395232042986368_r_kind,3.914960786286488_r_kind,4.473383886731226_r_kind,5.120188742634101_r_kind,5.810448401715452_r_kind,6.569783191416969_r_kind,7.483492057117051_r_kind,8.601373977842535_r_kind,& + 9.829200088657693_r_kind,11.088750547321505_r_kind,12.55142744859476_r_kind,14.1141400888864_r_kind,15.95433042449844_r_kind,18.086576629372065_r_kind,20.52912985532126_r_kind,23.37573439745639_r_kind,26.57923445765057_r_kind,& + 30.817348662795993_r_kind,35.82459590118845_r_kind,41.6554853313504_r_kind /) + + + interp_impact(:) = (/ 2._r_single, 2.5_r_single, 3._r_single,3.5_r_single, 4._r_single, 4.5_r_single, 5._r_single, 5.5_r_single,6._r_single, & + 6.5_r_single, 7._r_single, 7.5_r_single,8._r_single, 8.5_r_single, 9._r_single, 9.5_r_single, 10._r_single,11._r_single, & + 12._r_single, 13._r_single, 14._r_single,15._r_single, 16._r_single, 17._r_single, 18._r_single, 19._r_single,20._r_single, & + 21._r_single, 22._r_single, 23._r_single,24._r_single, 25._r_single, 26._r_single, 27._r_single, 28._r_single,29._r_single, & + 30._r_single, 31._r_single, 32._r_single,33._r_single, 34._r_single, 35._r_single, 36._r_single, 37._r_single,38._r_single, 39._r_single,& + 40._r_single, 41._r_single, 42._r_single, 43._r_single, 44._r_single, 45._r_single, 46._r_single, 47._r_single,48._r_single,49._r_single,& + 50._r_single, 51._r_single, 52._r_single, 53._r_single,54._r_single, 55._r_single, 56._r_single, 57._r_single, 58._r_single,59._r_single, 60._r_single /) + !!!!!!!!!!!!!!! + min_error = 1.25_r_kind + + if (input_impact <= interp_impact(1)) then + error = a(1)*(interp_impact(1) - interp_impact(1))**3 + b(1)*(interp_impact(1)-interp_impact(1))**2 + & + c(1)*(interp_impact(1)-interp_impact(1)) + d(1) + else if (input_impact > interp_impact(67)) then + error = 48.79740724637614_r_kind + else + level_index = 1 + do l=1,size(interp_impact) + if (interp_impact(l) >= input_impact) then + level_index = l + exit + else + continue + endif + enddo + interp_err = a(level_index-1)*(input_impact - interp_impact(level_index-1))**3 + & + b(level_index-1)*(input_impact-interp_impact(level_index-1))**2 + & + c(level_index-1)*(input_impact-interp_impact(level_index-1)) + d(level_index-1) + if (input_impact > 9.0_r_kind .and. input_impact < 11.0_r_kind) then + error = ((11.0_r_kind - input_impact)/2.0_r_kind)*interp_err + ((input_impact - 9.0_r_kind)/2.0_r_kind)*min_error + else if (input_impact > 29.0_r_kind .and. input_impact < 31.0_r_kind) then + error = ((31.0_r_kind - input_impact)/2.0_r_kind)*min_error + ((input_impact - 29.0_r_kind)/2.0_r_kind)*interp_err + else + if (input_impact <= 9.0_r_kind .or. input_impact >= 31.0_r_kind) then + error = interp_err + else + error = min_error + endif + endif + endif + return +end subroutine compute_threeCH_uncert subroutine init_netcdf_diag_ integer(i_kind) ncd_fileid, ncd_nobs @@ -788,6 +1224,10 @@ 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_to_single("Frac_LSW", gps_allptr%rdiag(23) ) + call nc_diag_metadata_to_single("STD4060", gps_allptr%rdiag(24) ) + call nc_diag_metadata_to_single("LSWFlag", gps_allptr%rdiag(25) ) + 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)) diff --git a/src/gsi/m_gpsrhs.F90 b/src/gsi/m_gpsrhs.F90 index 80f0b5b83..2049f6efe 100644 --- a/src/gsi/m_gpsrhs.F90 +++ b/src/gsi/m_gpsrhs.F90 @@ -68,11 +68,16 @@ module m_gpsrhs public:: cdiagbuf public:: qcfail - public:: qcfail_loc - public:: qcfail_high - public:: qcfail_gross public:: qcfail_jac + public:: qcfail_one + public:: qcfail_two + public:: qcfail_three + public:: qcfail_five + public:: qcfail_six + public:: qcfail_seven + public:: qcfail_eight + public:: data_ier public:: data_igps public:: data_ihgt @@ -117,11 +122,16 @@ module m_gpsrhs character(len=8), pointer, dimension( :):: cdiagbuf => null() logical , pointer, dimension( :):: qcfail => null() - real(r_single ), pointer, dimension( :):: qcfail_loc => null() - real(r_single ), pointer, dimension( :):: qcfail_high => null() - real(r_single ), pointer, dimension( :):: qcfail_gross=> null() real(r_single ), pointer, dimension( :):: qcfail_jac=> null() + real(r_single ), pointer, dimension( :):: qcfail_one => null() + real(r_single ), pointer, dimension( :):: qcfail_two => null() + real(r_single ), pointer, dimension( :):: qcfail_three=> null() + real(r_single ), pointer, dimension( :):: qcfail_five=> null() + real(r_single ), pointer, dimension( :):: qcfail_six=> null() + real(r_single ), pointer, dimension( :):: qcfail_seven=> null() + real(r_single ), pointer, dimension( :):: qcfail_eight=> null() + real(r_kind ), pointer, dimension( :):: data_ier => null() real(r_kind ), pointer, dimension( :):: data_igps => null() real(r_kind ), pointer, dimension( :):: data_ihgt => null() @@ -152,7 +162,10 @@ module m_gpsrhs character(len=8), pointer, dimension( :), save:: cdiagbuf logical , pointer, dimension( :), save:: qcfail - real(r_single ), pointer, dimension( :), save:: qcfail_loc,qcfail_high,qcfail_gross,qcfail_jac + real(r_single ), pointer, dimension( :), save:: qcfail_jac + + real(r_single ), pointer, dimension( :), save:: qcfail_one,qcfail_two,qcfail_three,qcfail_five + real(r_single ), pointer, dimension( :), save:: qcfail_six,qcfail_seven,qcfail_eight real(r_kind ), pointer, dimension( :), save:: data_ier real(r_kind ), pointer, dimension( :), save:: data_igps @@ -264,17 +277,27 @@ subroutine gpsrhs_alloc(is,class,nobs,nsig,nreal,grids_dim,nsig_ext) b%cdiagbuf( :)="" allocate(b%qcfail (nobs)) - allocate(b%qcfail_loc (nobs)) - allocate(b%qcfail_high (nobs)) - allocate(b%qcfail_gross (nobs)) allocate(b%qcfail_jac (nobs)) + allocate(b%qcfail_one (nobs)) + allocate(b%qcfail_two (nobs)) + allocate(b%qcfail_three (nobs)) + allocate(b%qcfail_five (nobs)) + allocate(b%qcfail_six (nobs)) + allocate(b%qcfail_seven (nobs)) + allocate(b%qcfail_eight (nobs)) + b%qcfail=.false. - b%qcfail_loc =zero - b%qcfail_high =zero - b%qcfail_gross =zero b%qcfail_jac =zero + b%qcfail_one =zero + b%qcfail_two =zero + b%qcfail_three =zero + b%qcfail_five =zero + b%qcfail_six =zero + b%qcfail_seven =zero + b%qcfail_eight =zero + allocate(b%data_ier (nobs)) allocate(b%data_igps(nobs)) allocate(b%data_ihgt(nobs)) @@ -347,11 +370,16 @@ subroutine gpsrhs_dealloc(is) deallocate(b%cdiagbuf) deallocate(b%qcfail ) - deallocate(b%qcfail_loc ) - deallocate(b%qcfail_high ) - deallocate(b%qcfail_gross ) deallocate(b%qcfail_jac ) + deallocate(b%qcfail_one ) + deallocate(b%qcfail_two ) + deallocate(b%qcfail_three ) + deallocate(b%qcfail_five ) + deallocate(b%qcfail_six ) + deallocate(b%qcfail_seven ) + deallocate(b%qcfail_eight ) + deallocate(b%data_ier ) deallocate(b%data_igps) deallocate(b%data_ihgt) @@ -405,11 +433,16 @@ subroutine gpsrhs_aliases(is) cdiagbuf => b%cdiagbuf qcfail => b%qcfail - qcfail_loc => b%qcfail_loc - qcfail_high => b%qcfail_high - qcfail_gross => b%qcfail_gross qcfail_jac => b%qcfail_jac + qcfail_one => b%qcfail_one + qcfail_two => b%qcfail_two + qcfail_three => b%qcfail_three + qcfail_five => b%qcfail_five + qcfail_six => b%qcfail_six + qcfail_seven => b%qcfail_seven + qcfail_eight => b%qcfail_eight + data_ier => b%data_ier data_igps => b%data_igps data_ihgt => b%data_ihgt @@ -442,8 +475,10 @@ subroutine gpsrhs_unaliases(is) nullify(rges,gp2gm,prsltmp_o,tges_o) nullify(error,error_adjst,ratio_errors) nullify(rdiagbuf,cdiagbuf) - nullify(qcfail,qcfail_loc,qcfail_gross,qcfail_jac) - nullify(qcfail_high) + nullify(qcfail,qcfail_jac) + + nullify(qcfail_one,qcfail_two,qcfail_three,qcfail_five) + nullify(qcfail_six,qcfail_seven,qcfail_eight) nullify(data_ier,data_igps,data_ihgt) _EXIT_(myname_) end subroutine gpsrhs_unaliases diff --git a/src/gsi/obsmod.F90 b/src/gsi/obsmod.F90 index dade37f1a..b898190db 100644 --- a/src/gsi/obsmod.F90 +++ b/src/gsi/obsmod.F90 @@ -448,6 +448,7 @@ module obsmod public :: lsaveobsens public :: iout_cldch, mype_cldch public :: nprof_gps,time_offset,ianldate,tcp_box + public :: nobs_gps public :: iout_oz,iout_co,dsis,ref_obs,obsfile_all,lobserver,tcp_posmatch,perturb_obs,ditype,dsfcalc,dplat public :: time_window,dval,dtype,dfile,dirname,obs_setup,oberror_tune,offtime_data public :: lobsdiagsave,lobsdiag_forenkf,blacklst,hilbert_curve,lobskeep,time_window_max,sfcmodel,ext_sonde @@ -593,7 +594,7 @@ module obsmod real(r_kind),dimension(50):: dmesh real(r_kind) r_hgt_fed integer(i_kind) nchan_total,ianldate - integer(i_kind) ndat,ndat_types,ndat_times,nprof_gps + integer(i_kind) ndat,ndat_types,ndat_times,nprof_gps,nobs_gps integer(i_kind) lunobs_obs,nloz_v6,nloz_v8,nobskeep,nloz_omi integer(i_kind) nlco,use_limit integer(i_kind) iout_rad,iout_pcp,iout_t,iout_q,iout_uv, & @@ -944,6 +945,7 @@ subroutine init_obsmod_dflts ! precipitation rate observations nprof_gps = 0 + nobs_gps = 0 hilbert_curve=.false. neutral_stability_windfact_2dvar=.false. diff --git a/src/gsi/read_gps.f90 b/src/gsi/read_gps.f90 index c0cef658a..a21f6eb3f 100644 --- a/src/gsi/read_gps.f90 +++ b/src/gsi/read_gps.f90 @@ -61,6 +61,7 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & ! 2017-11-16 dutta - addition of profile quality flags for KOMPSAT5 GPSRO. ! 2019-08-21 Shao - add qc flags input for METOP-C, COSMIC-2 and PAZ ! 2020-05-21 Shao - add qc flags input for commercial GNSSRO data +! 2025-01-21 Li - add LSW variables ! ! input argument list: ! infile - unit from which to read BUFR data @@ -102,7 +103,7 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & ! Declare local parameters integer(i_kind),parameter:: maxlevs=500 - integer(i_kind),parameter:: maxinfo=16 + integer(i_kind),parameter:: maxinfo=18 real(r_kind),parameter:: r10000=10000.0_r_kind real(r_kind),parameter:: r360=360.0_r_kind @@ -129,12 +130,20 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & integer(i_kind),allocatable,dimension(:):: gpsro_itype,gpsro_ikx,nmrecs_id - + !!!! LSW Check Variables + real(r_kind) real_lsw,real_flsw,fraclsw,lsw_flag_hold + integer(i_kind) keep_level,keep_level_sorted,count_check + real(r_kind),allocatable,dimension(:) :: bend4060,lsw_flag,& + sorted_lsw,sorted_fraclsw,sorted_impact, & + test_qc,array_lsw,array_fraclsw,array_impact, & + sorted_lswflag + integer(i_kind),allocatable,dimension(:) :: indices + !!!! real(r_kind) timeo,t4dv real(r_kind) pcc,qfro,usage,dlat,dlat_earth,dlon,dlon_earth,freq_chk,freq real(r_kind) dlat_earth_deg,dlon_earth_deg real(r_kind) height,rlat,rlon,ref,bend,impact,roc,geoid,& - bend_error,ref_error,bend_pccf,ref_pccf + bend_error,ref_error,bend_pccf,ref_pccf,bend_lsw real(r_kind),allocatable,dimension(:,:):: cdata_all @@ -258,10 +267,11 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & endif ! Check profile quality flags +! GRACE-FO was generated by GFZ (the same as Terra-SAR and Tandem SAID = 42 and 43) if ( ((said > 739).and.(said < 746)).or.(said == 820).or.(said == 786).or. & ((said > 749).and.(said < 756)).or.(said == 825).or.(said == 44) .or. & (said == 265).or.(said == 266).or.(said == 267).or.(said == 268).or. & - (said == 269)) then !CDAAC processing + (said == 269).or.(said == 803).or.(said == 804).or.(said == 66)) then !CDAAC processing if(pcc==zero) then ! write(6,*)'READ_GPS: bad profile said=',said,'ptid=',ptid,& ! ' SKIP this report' @@ -340,6 +350,41 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & if(mod(nmrecs,ncnumgrp(ikx))== ncgroup(ikx)-1)usage=ncmiter(ikx) end if +! New LSW Check + allocate(array_lsw(levs),array_fraclsw(levs),array_impact(levs), & + sorted_lsw(levs),sorted_fraclsw(levs),sorted_impact(levs), & + indices(levs),lsw_flag(levs),sorted_lswflag(levs)) + lsw_flag(:) = 0.0_r_kind + sorted_lswflag(:) = 0.0_r_kind + do k=1, levs + do i=1,nreps_this_ROSEQ2(k) + m=(6*i)-2 + freq_chk=data1b(m,k) ! frequency (hertz) + if(nint(freq_chk).ne.0) cycle ! do not want non-zero freq., go on to next replication of ROSEQ2 + freq=data1b(m,k) + array_impact(k)=(data1b(m+1,k) - roc)/1000._r_kind ! impactparameter (m) + array_lsw(k) = data1b(m+4,k) ! RMSE in bending angle (rad) + array_fraclsw(k) = (data1b(m+4,k)/data1b(m+2,k))*100._r_kind + enddo + enddo + ! + call sort(levs,array_impact,indices) + do k=1,levs + sorted_impact(k) = array_impact(indices(k)) + sorted_lsw(k) = array_lsw(indices(k)) + sorted_fraclsw(k) = array_fraclsw(indices(k)) + if (array_lsw(indices(k)) > 1.0e8 .and. array_impact(indices(k)) <= 10._r_kind) then + lsw_flag(indices(k)) = 2.0_r_kind + sorted_lswflag(k) = 2.0_r_kind + else if (array_fraclsw(indices(k)) > 35._r_kind .and. array_lsw(indices(k)) < 1.0e8 & + .and. array_impact(indices(k)) <= 10._r_kind) then + lsw_flag(indices(k)) = 1.0_r_kind + sorted_lswflag(k) = 1.0_r_kind + endif + enddo + deallocate(array_lsw,array_fraclsw,array_impact, & + sorted_lsw,sorted_fraclsw,sorted_impact,indices,sorted_lswflag) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Loop over levs in profile do k=1, levs nread=nread+1 ! count observations @@ -359,6 +404,8 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & impact=data1b(m+1,k) ! impact parameter (m) bend=data1b(m+2,k) ! bending angle (rad) bend_error=data1b(m+4,k) ! RMSE in bending angle (rad) + bend_lsw = data1b(m+4,k) ! Bend Lsw + lsw_flag_hold = lsw_flag(k) enddo bend_pccf=data1b((6*nreps_this_ROSEQ2(k))+4,k) ! percent confidence for this ROSEQ1 replication @@ -368,7 +415,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 @@ -420,11 +468,15 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & cdata_all(4,ndata) = height ! geometric height above geoid (m) cdata_all(5,ndata) = ref ! refractivity obs (units of N) ! cdata_all(9,ndata) = ref_pccf ! per cent confidence (%) + cdata_all(17,ndata) = 0.0_r_kind + cdata_all(18,ndata) = -99.0_r_kind else cdata_all(1,ndata) = bend_error ! gps bending error (radians) cdata_all(4,ndata) = impact ! impact parameter (m) cdata_all(5,ndata) = bend ! bending angle obs (radians) ! cdata_all(9,ndata) = bend_pccf ! per cent confidence (%) + cdata_all(17,ndata) = bend_lsw + cdata_all(18,ndata) = lsw_flag_hold endif cdata_all(9,ndata) = pcc ! profile per cent confidence (0 or 100) cdata_all(2,ndata) = dlon ! grid relative longitude @@ -447,7 +499,7 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & ! End of k loop over levs end do - + deallocate(lsw_flag) enddo read_loop ! subsets enddo ! messages @@ -475,6 +527,65 @@ subroutine read_gps(nread,ndata,nodata,infile,lunout,obstype,twind, & deallocate(gpsro_ctype,gpsro_itype,gpsro_ikx,nmrecs_id) return +contains +subroutine sort(num,x,indx) +integer(i_kind), intent(in) :: num +real(r_kind), intent(in) :: x(num) +integer(i_kind), intent(inout) :: indx(num) + +integer(i_kind) :: ind, i, j, value_indx, line +real(r_kind) :: values + +! Initialize the index array to input order +do i = 1, num + indx(i) = i +end do + +! Only one element, just send it back +if(num <= 1) return + +line = num / 2 + 1 +ind = num + +! Keep looping until finished +do + ! Keep going down lines until bottom + if(line > 1) then + line = line - 1 + values = x(indx(line)) + value_indx = indx(line) + else + values = x(indx(ind)) + value_indx = indx(ind) + + indx(ind) = indx(1) + ind = ind - 1 + if(ind == 1) then + indx(1) = value_indx + return + endif + endif + + i = line + j = 2 * line + + do while(j <= ind) + if(j < ind) then + if(x(indx(j)) < x(indx(j + 1))) j = j + 1 + endif + if(values < x(indx(j))) then + indx(i) = indx(j) + i = j + j = 2 * j + else + j = ind + 1 + endif + end do + + indx(i) = value_indx +end do +end subroutine sort + end subroutine read_gps diff --git a/src/gsi/read_obs.F90 b/src/gsi/read_obs.F90 index 05c15b0ef..984e72aa8 100644 --- a/src/gsi/read_obs.F90 +++ b/src/gsi/read_obs.F90 @@ -405,8 +405,9 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) (said == 44) .or. (said == 5) .or. (said == 41) .or. & (said == 42) .or. (said == 43) .or. (said == 722) .or. & (said == 723).or. (said == 265).or. (said == 266) .or. & - (said == 267).or. (said == 268).or. (said == 269)) then - lexist=.true. + (said == 267).or. (said == 268).or. (said == 269) .or. & + (said == 803).or. (said == 804).or. (said == 66)) then + lexist=.true. exit gpsloop end if nread = nread + 1 @@ -733,7 +734,8 @@ subroutine read_obs(ndata,mype) use obsmod, only: iadate,ndat,time_window,dplat,dsfcalc,dfile,dthin, & dtype,dval,dmesh,obsfile_all,ref_obs,nprof_gps,dsis,ditype,& perturb_obs,lobserver,lread_obs_save,obs_input_common, & - reduce_diag,nobs_sub,dval_use,hurricane_radar,l2rwthin + reduce_diag,nobs_sub,dval_use,hurricane_radar,l2rwthin, & + nobs_gps use gsi_nstcouplermod, only: nst_gsi ! use gsi_nstcouplermod, only: gsi_nstcoupler_set use hdraobmod, only: read_hdraob,nhdt,nhdq,nhduv,nhdps,hdtlist,hdqlist,hduvlist,hdpslist,nodet,nodeq,nodeuv,nodeps @@ -799,6 +801,7 @@ subroutine read_obs(ndata,mype) character(20):: sis integer(i_kind) i,j,k,ii,nmind,lunout,isfcalc,ithinx,ithin,nread,npuse,nouse integer(i_kind) nprof_gps1,npem1,krsize,len4file,npemax,ilarge,nlarge,npestart + integer(i_kind) nobs_gps1 integer(i_llong) :: lenbytes integer(i_kind):: npetot,npeextra,mmdat,nodata integer(i_kind):: iworld,iworld_group,next_mype,mm1,iix @@ -845,6 +848,7 @@ subroutine read_obs(ndata,mype) end do npem1=npe-1 nprof_gps1=0 + nobs_gps1 = 0 if(njqc) then call converr_ps_read(mype) @@ -1930,6 +1934,7 @@ subroutine read_obs(ndata,mype) else if (ditype(i) == 'gps')then call read_gps(nread,npuse,nouse,infile,lunout,obstype,twind, & nprof_gps1,sis,nobs_sub1(1,i)) + nobs_gps1 = nouse string='READ_GPS' ! Process aerosol data @@ -2042,6 +2047,7 @@ subroutine read_obs(ndata,mype) ! Collect number of gps profiles (needed later for qc) call mpi_allreduce(nprof_gps1,nprof_gps,1,mpi_integer,mpi_sum,mpi_comm_world,ierror) + call mpi_allreduce(nobs_gps1,nobs_gps,1,mpi_integer,mpi_sum,mpi_comm_world,ierror) call mpi_allreduce(nobs_sub1,nobs_sub,npe*ndat,mpi_integer,mpi_sum,mpi_comm_world,& ierror) diff --git a/src/gsi/setupbend.f90 b/src/gsi/setupbend.f90 index 607311f34..a29200608 100644 --- a/src/gsi/setupbend.f90 +++ b/src/gsi/setupbend.f90 @@ -108,6 +108,9 @@ subroutine setupbend(obsLL,odiagLL, & ! 2021-11-05 cucurull - update QCs and optimize/improve forward operator; bug fixes ! 2022-01-28 cucurull - add Sentinel-6, PAZ ! 2022-04-06 collard - reintroduce Jacbian QC as an option (default off) +! 2024-12-04 Li - turn on MetOp data below 8 km +! 2024-12-04 Li - add GRACE-FO (803&804) data below 8 km +! 2024-12-04 Li - add new obs error model by Chris Riedel ! ! input argument list: ! lunin - unit from which to read observations @@ -168,7 +171,9 @@ subroutine setupbend(obsLL,odiagLL, & use m_gpsrhs, only: ratio_errors use m_gpsrhs, only: rdiagbuf,cdiagbuf use m_gpsrhs, only: qcfail - use m_gpsrhs, only: qcfail_loc,qcfail_high,qcfail_gross,qcfail_jac + use m_gpsrhs, only: qcfail_jac + use m_gpsrhs, only: qcfail_one,qcfail_two,qcfail_three,qcfail_five + use m_gpsrhs, only: qcfail_six,qcfail_seven,qcfail_eight use m_gpsrhs, only: data_ier,data_igps,data_ihgt use m_gpsrhs, only: gpsrhs_alloc use m_gpsrhs, only: gpsrhs_dealloc @@ -196,6 +201,7 @@ subroutine setupbend(obsLL,odiagLL, & ! Declare local parameters real(r_kind),parameter:: r240 = 240.0_r_kind real(r_kind),parameter:: six = 6.0_r_kind + real(r_kind),parameter:: seven = 7.0_r_kind real(r_kind),parameter:: ten = 10.0_r_kind real(r_kind),parameter:: eight = 8.0_r_kind real(r_kind),parameter:: nine = 9.0_r_kind @@ -228,7 +234,8 @@ subroutine setupbend(obsLL,odiagLL, & real(r_kind),dimension(4) :: w4,dw4,dw4_TL integer(i_kind) ier,ilon,ilat,ihgt,igps,itime,ikx,iuse, & - iprof,ipctc,iroc,isatid,iptid,ilate,ilone,ioff,igeoid + iprof,ipctc,iroc,isatid,iptid,ilate,ilone,ioff,igeoid, & + ilsw,ilswflag integer(i_kind) i,j,k,kk,mreal,nreal,jj,ikxx,ibin integer(i_kind) mm1,nsig_up,ihob,istatus,nsigstart integer(i_kind) kprof,istat,k1,k2,nobs_out,top_layer_SR,bot_layer_SR,count_SR @@ -301,6 +308,8 @@ subroutine setupbend(obsLL,odiagLL, & !268 => PlanetiQ GNOMES-B !269 => Spire Lemur 3U CubeSat !66 => Sentinel-6 +!803=> GRACE C +!804=> GRACE D ! Check to see if required guess fields are available call check_vars_(proceed) @@ -328,6 +337,8 @@ 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) + ilsw=17 ! index of bending angle LSW + ilswflag=18 ! index of LSW flag ! Intialize variables nsig_up=nsig+nsig_ext ! extend nsig_ext levels above interface level nsig @@ -345,7 +356,7 @@ 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=25 nreal=mreal if (lobsdiagsave) nreal=nreal+4*miter+1 if (save_jacobian) then @@ -397,9 +408,12 @@ subroutine setupbend(obsLL,odiagLL, & muse(:)=.false. qcfail=.false. - qcfail_loc=zero;qcfail_gross=zero - qcfail_high=zero qcfail_jac=zero + + qcfail_one=zero;qcfail_two=zero + qcfail_three=zero;qcfail_five=zero + qcfail_six=zero;qcfail_seven=zero + qcfail_eight=zero toss_gps_sub=zero dbend_loc=zero @@ -550,7 +564,7 @@ subroutine setupbend(obsLL,odiagLL, & data(ier,i) = zero ratio_errors(i) = zero muse(i)=.false. - qcfail_loc(i)=one + qcfail_one(i) = one endif ! Increment obs counter along with low and high obs counters @@ -582,6 +596,9 @@ 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) = (data(ilsw,i)/data(igps,i))*100._r_kind ! Fractional LSW + rdiagbuf(24,i) = 0.0_r_kind ! holder for STD4060 in hybrid error model + rdiagbuf(25,i) = data(ilswflag,i) ! LSW flag if(ratio_errors(i) > tiny_r_kind) then ! obs inside model grid @@ -639,7 +656,8 @@ subroutine setupbend(obsLL,odiagLL, & (data(isatid,i)==42) .or.(data(isatid,i)==3) .or. & (data(isatid,i)==821).or.(data(isatid,i)==421).or. & (data(isatid,i)==440).or.(data(isatid,i)==43) .or. & - (data(isatid,i)==5).or.(data(isatid,i)==66)) then + (data(isatid,i)==5).or.(data(isatid,i)==66) .or. & + (data(isatid,i)==803).or.(data(isatid,i)==804)) then if((data(ilate,i)> r40).or.(data(ilate,i)< -r40)) then if(alt>r12) then @@ -754,6 +772,7 @@ subroutine setupbend(obsLL,odiagLL, & ddnj(j)=dot_product(dw4,nrefges(ihob-1:ihob+2,i))!derivative (dN/dx)_j if(ddnj(j)>zero) then qcfail(i)=.true. + qcfail_eight(i) = one data(ier,i) = zero ratio_errors(i) = zero muse(i)=.false. @@ -802,6 +821,7 @@ subroutine setupbend(obsLL,odiagLL, & if(dbend > 0.05_r_kind) then data(ier,i) = zero ratio_errors(i) = zero + qcfail_five(i) = one qcfail(i)=.true. muse(i)=.false. endif @@ -821,7 +841,7 @@ subroutine setupbend(obsLL,odiagLL, & if (luse(i)) then awork(4) = awork(4)+one endif - qcfail_gross(i)=one + qcfail_three(i) = one data(ier,i) = zero ratio_errors(i) = zero muse(i)=.false. @@ -871,6 +891,7 @@ subroutine setupbend(obsLL,odiagLL, & if(abs(rdiagbuf(5,i)) > cutoff) then qcfail(i)=.true. + qcfail_six(i) = one data(ier,i) = zero ratio_errors(i) = zero muse(i) = .false. @@ -882,18 +903,19 @@ subroutine setupbend(obsLL,odiagLL, & if((alt > gpstop) .or. (commdat .and. (alt > commgpstop))) then data(ier,i) = zero ratio_errors(i) = zero - qcfail_high(i)=one + qcfail_two(i)=one muse(i)=.false. endif -! Remove MetOP/GRAS data below 8 km - if( (alt <= eight) .and. & - ((data(isatid,i)==4).or.(data(isatid,i)==3).or.(data(isatid,i)==5))) then - qcfail(i)=.true. - data(ier,i) = zero - ratio_errors(i) = zero - muse(i)=.false. - endif +! Turn on MetOP/GRAS data below 8 km +! if( (alt <= eight) .and. & +! ((data(isatid,i)==4).or.(data(isatid,i)==3).or.(data(isatid,i)==5))) then +! qcfail(i)=.true. +! qcfail_seven(i) = one +! data(ier,i) = zero +! ratio_errors(i) = zero +! muse(i)=.false. +! endif end if ! obs above super-refraction and shadow layers end if ! obs inside the vertical grid @@ -911,7 +933,9 @@ subroutine setupbend(obsLL,odiagLL, & do i=1,nobs - if (qcfail(i)) then + if (qcfail(i) .or. qcfail_five(i) > zero .or. & + qcfail_six(i) > zero .or. qcfail_seven(i) > zero .or. & + qcfail_eight(i) > zero) then data(ier,i) = zero ratio_errors(i) = zero muse(i) = .false. @@ -947,10 +971,14 @@ subroutine setupbend(obsLL,odiagLL, & ! flags for observations that failed qc checks ! zero = observation is good - if(qcfail_gross(i) == one) rdiagbuf(10,i) = three - 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(qcfail(i)) rdiagbuf(10,i) = four !modified in genstats due to toss_gps_sub + if(qcfail_one(i) == one) rdiagbuf(10,i) = one + if(qcfail_two(i) == one) rdiagbuf(10,i) = two + if(qcfail_three(i) == one) rdiagbuf(10,i) = three + if(qcfail_five(i) == one) rdiagbuf(10,i) = five + if(qcfail_six(i) == one) rdiagbuf(10,i) = six + if(qcfail_seven(i) == one) rdiagbuf(10,i) = seven + if(qcfail_eight(i) == one) rdiagbuf(10,i) = eight 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) diff --git a/src/gsi/setupref.f90 b/src/gsi/setupref.f90 index 4f240d756..cf5e267d9 100644 --- a/src/gsi/setupref.f90 +++ b/src/gsi/setupref.f90 @@ -164,7 +164,7 @@ subroutine setupref(obsLL,odiagLL,lunin,mype,awork,nele,nobs,toss_gps_sub,is,ini use m_gpsrhs, only: ratio_errors,dpresl use m_gpsrhs, only: rdiagbuf,cdiagbuf use m_gpsrhs, only: qcfail - use m_gpsrhs, only: qcfail_loc,qcfail_high,qcfail_gross + use m_gpsrhs, only: qcfail_one,qcfail_two,qcfail_three use m_gpsrhs, only: data_ier,data_igps,data_ihgt use m_gpsrhs, only: gpsrhs_alloc use m_gpsrhs, only: gpsrhs_dealloc @@ -341,9 +341,9 @@ subroutine setupref(obsLL,odiagLL,lunin,mype,awork,nele,nobs,toss_gps_sub,is,ini ! these statements should moved into gpsrhs_alloc(), but ! left unchanged for verification purposes. qcfail=.false. - qcfail_loc=zero - qcfail_gross=zero - qcfail_high=zero + qcfail_one=zero + qcfail_three=zero + qcfail_two=zero muse(:)=.false. @@ -494,7 +494,7 @@ subroutine setupref(obsLL,odiagLL,lunin,mype,awork,nele,nobs,toss_gps_sub,is,ini data(ier,i) = zero ratio_errors(i) = zero muse(i)=.false. - qcfail_loc(i)=one + qcfail_one(i)=one endif ! Increment obs counter along with low and high obs counters @@ -565,7 +565,7 @@ subroutine setupref(obsLL,odiagLL,lunin,mype,awork,nele,nobs,toss_gps_sub,is,ini if (luse(i)) then awork(4) = awork(4)+one endif - qcfail_gross(i)=one + qcfail_three(i)=one data(ier,i) = zero ratio_errors(i) = zero muse(i)=.false. @@ -613,7 +613,7 @@ subroutine setupref(obsLL,odiagLL,lunin,mype,awork,nele,nobs,toss_gps_sub,is,ini if(alt > gpstop) then data(ier,i) = zero ratio_errors(i) = zero - qcfail_high(i)=one + qcfail_two(i)=one muse(i)=.false. endif @@ -690,13 +690,13 @@ subroutine setupref(obsLL,odiagLL,lunin,mype,awork,nele,nobs,toss_gps_sub,is,ini if(r1em3*rdiagbuf(7,i) <= 3.0_r_kind) then ! check for SR-likely conditions below 3 km kprof = data(iprof,i) jprof = data(iprof,i+1) - if( kprof == jprof .and. .not. qcfail(i) .and. qcfail_loc(i) == zero .and. qcfail_loc(i+1) == zero .and. qcfail_gross(i) == zero)then + if( kprof == jprof .and. .not. qcfail(i) .and. qcfail_one(i) == zero .and. qcfail_one(i+1) == zero .and. qcfail_three(i) == zero)then grad_mod=1000.0_r_kind*& ((rdiagbuf(17,i+1)*(one-rdiagbuf(5,i+1)))-(rdiagbuf(17,i)*(one-rdiagbuf(5,i))))/(rdiagbuf(7,i+1)-rdiagbuf(7,i)) grad_obs=1000.0_r_kind*(rdiagbuf(17,i+1)-rdiagbuf(17,i))/(rdiagbuf(7,i+1)-rdiagbuf(7,i)) if ((abs(grad_mod)>= half*crit_grad) .or. (abs(grad_obs)>=half*crit_grad)) then qcfail(i) = .true. - if( qcfail_gross(i+1) == zero .and. .not. qcfail(i+1)) then + if( qcfail_three(i+1) == zero .and. .not. qcfail(i+1)) then qcfail(i+1)= .true. end if end if @@ -707,7 +707,7 @@ subroutine setupref(obsLL,odiagLL,lunin,mype,awork,nele,nobs,toss_gps_sub,is,ini kprof = data(iprof,i) do j=1,nobs jprof = data(iprof,j) - if( kprof == jprof .and. .not. qcfail(j) .and. qcfail_loc(j) == zero .and. qcfail_gross(j) == zero)then + if( kprof == jprof .and. .not. qcfail(j) .and. qcfail_one(j) == zero .and. qcfail_three(j) == zero)then ! Remove data below if(r1em3*rdiagbuf(7,j) < r1em3*rdiagbuf(7,i))then @@ -780,10 +780,10 @@ subroutine setupref(obsLL,odiagLL,lunin,mype,awork,nele,nobs,toss_gps_sub,is,ini ! flags for observations that failed qc checks ! zero = observation is good - if(qcfail_gross(i) == one) rdiagbuf(10,i) = three + if(qcfail_three(i) == one) rdiagbuf(10,i) = three 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(qcfail_one(i) == one) rdiagbuf(10,i) = one + if(qcfail_two(i) == one) rdiagbuf(10,i) = two 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)