Index: arch/preamble =================================================================== --- arch/preamble (revision 374) +++ arch/preamble (working copy) @@ -28,7 +28,7 @@ NETCDF_INC = -I${NETCDF}/include -NCARG_LIBS = -L${NCARG_ROOT}/lib -lncarg -lncarg_gks -lncarg_c -lX11 -lm +NCARG_LIBS = -L${NCARG_ROOT}/lib -lncarg -lncarg_gks -lncarg_c -lX11 -lm -lcairo #### Architecture specific settings #### Index: namelist.oa =================================================================== --- namelist.oa (revision 374) +++ namelist.oa (working copy) @@ -1,21 +1,22 @@ &record1 - start_year = 2000 - start_month = 01 - start_day = 24 + start_year = 2006 + start_month = 12 + start_day = 19 start_hour = 12 - end_year = 2000 - end_month = 01 - end_day = 25 + end_year = 2006 + end_month = 12 + end_day = 20 end_hour = 12 - interval = 21600 + interval = 10800 / &record2 grid_id = 1 - obs_filename = '../JAN00/obs' + obs_filename = './obs' remove_data_above_qc_flag = 32768 remove_unverified_data = .TRUE. / + obs_filename = '../JAN00/obs' trim_domain = .TRUE. trim_value = 5 @@ -43,6 +44,11 @@ max_p_extend_t = 1300 max_p_extend_w = 1300 / + max_error_dewpoint = 20 + max_buddy_dewpoint = 20 + qc_psfc = .TRUE. + use_p_tolerance_one_lev = .TRUE. + max_p_tolerance_one_lev_qc = 700 &record5 print_obs_files = .TRUE. @@ -86,6 +92,10 @@ oa_3D_option = 1 oa_3D_type = 'Cressman' radius_influence = 5,4,3,2, + scale_cressman_rh_decreases = .TRUE. + radius_influence_sfc_mult = 0.4 + oa_psfc = .FALSE. + max_p_tolerance_one_lev_oa = 700 &plot_sounding Index: src/Makefile =================================================================== --- src/Makefile (revision 374) +++ src/Makefile (working copy) @@ -9,7 +9,7 @@ SRC = driver.F error_handler.F linpak.F main.F \ proc_final_analysis.F proc_first_guess.F proc_get_info_header.F \ proc_header.F proc_namelist.F proc_oa.F \ - proc_ob_access.F proc_obs_sort.F proc_qc.F sample.F yx2xy.F + proc_ob_access.F proc_obs_sort.F proc_qc.F sample.F yx2xy.F OBJS = $(SRC:.F=.o) @@ -19,7 +19,7 @@ namelist_module.o oa_module.o \ ob_access_module.o obs_sort_module.o printout_module.o \ qc0_module.o qc1_module.o qc2_module.o qc3_module.o \ - small_header_data_module.o tinterp_module.o + small_header_data_module.o tinterp_module.o get_fg_at_point_module.o INCS = big_header.inc constants.inc error.inc first_guess.inc first_guess_size.inc \ map.inc missing.inc namelist.inc proc_get_info_header.inc @@ -81,6 +81,7 @@ first_guess_module.o: error.inc error.int first_guess_size.inc first_guess.inc big_header.inc \ input_data_module.o date_pack_module.o small_header_data_module.o \ oa_module.o +get_fg_at_point_module.o: missing.inc constants.inc map_utils_module.o map_utils_helper_module.o input_data_module.o: small_header_data_module.o linpak.o: main.o: error.inc big_header.inc proc_header.int error.int \ @@ -91,11 +92,13 @@ make_sfc_fdda.o: diags_module.o namelist_module.o: error.inc error.int namelist.inc namelist.common oa_module.o: error.int error.inc first_guess_size.inc first_guess.inc -ob_access_module.o: error.int error.inc missing.inc \ - obs_sort_module.o date_pack_module.o +ob_access_module.o: error.int error.inc missing.inc constants.inc \ + obs_sort_module.o date_pack_module.o map_utils_module.o map_utils_helper_module.o \ + get_fg_at_point_module.o obs_sort_module.o: missing.inc error.inc error.int proc_get_info_header.inc map.inc \ constants.inc \ - date_pack_module.o map_utils_module.o map_utils_helper_module.o + date_pack_module.o map_utils_module.o map_utils_helper_module.o \ + get_fg_at_point_module.o printout_module.o: proc_final_analysis.o: first_guess_size.inc first_guess.inc proc_get_info_header.inc \ proc_get_info_header.int sample.int header.common \ @@ -117,7 +120,8 @@ qc0_module.o: obs_sort_module.o qc1_module.o: qc0_module.o qc2_module.o: qc0_module.o -qc3_module.o: obs_sort_module.o +qc3_module.o: constants.inc \ + obs_sort_module.o qc0_module.o sample.o: printout_module.o small_header_data_module.o: tinterp_module.o: small_header_data_module.o input_data_module.o Index: src/constants.inc =================================================================== --- src/constants.inc (revision 374) +++ src/constants.inc (working copy) @@ -2,3 +2,10 @@ cp = 1004. , & ! specific heat at constant pressure gasr = 287.95 , & ! gas constant for dry air rcp = gasr / cp +!BPR BEGIN + REAL , PARAMETER :: L_over_Rv = 5418.12 + REAL , PARAMETER :: Rv_over_L = 1. / L_over_Rv + REAL , PARAMETER :: t_at_sea_level_std_atm = 288.15 + REAL , PARAMETER :: dTdz_std_atm_lt_11000m = -0.0065 + REAL , PARAMETER :: very_small_rh = 0.0001 +!BPR END Index: src/driver.F90 =================================================================== --- src/driver.F90 (revision 374) +++ src/driver.F90 (working copy) @@ -43,6 +43,16 @@ TYPE (report) :: obs_sort_tmp CHARACTER ( LEN = 132 ) :: dummy_filename, tmp_filename + ! BPR BEGIN + TYPE (report) , ALLOCATABLE , DIMENSION ( : ) :: obs_before_sort + INTEGER , PARAMETER :: INT_14DIGITS = selected_int_kind ( 14 ) + INTEGER (kind=INT_14DIGITS) , ALLOCATABLE , DIMENSION ( : ) :: date_int + INTEGER , ALLOCATABLE , DIMENSION ( : ) :: obs_order_index + INTEGER :: obs_order_sort_tmp + INTEGER :: date_int_sort_tmp + INTEGER :: sort_type + ! BPR END + ! The NAMELIST variables, with a small amount of error checks ! from the main program. @@ -85,6 +95,9 @@ INTEGER :: obs_file_count INTEGER :: mqd_count , mqd_abs_min + !BPR BEGIN + LOGICAL :: use_grid_relative_winds + !BPR END ! If we are doing sfc FDDA, the first time in we do the traditional analysis AND the surface ! FDDA fields for that time (fdda_loop_max=2). All of the subsequent times, we process the @@ -225,6 +238,7 @@ iew_alloc , jns_alloc , kbu_alloc , pressure , & nml%record_5%print_analysis , & current_date_8 , current_time_6 , date_char , icount ) + #ifdef PRESSURE pressure(:kbu_alloc) = pressure(:kbu_alloc) * 100. #endif @@ -251,6 +265,25 @@ dxd = dxd * 0.001 END IF + ! BPR BEGIN + ! If the user chose to implement pressure tolerances for single-level + ! above-surface obs AND the tolerance is really a tolerance (i.e. > 1 + ! Pa) then print out the effective pressure coverage given the + ! tolerances and the pressures on which the first guess field exists + ! Only do this once per Obsgrid run. This assumes that the vertical + ! levels in the first-guess field are the same for each time in a given + ! Obsgrid run. + IF ( icount .EQ. 1 ) THEN + IF( nml%record_4%use_p_tolerance_one_lev .AND. & + ( ( nml%record_4%max_p_tolerance_one_lev_qc .GT. 1 ) .OR. & + ( nml%record_9%max_p_tolerance_one_lev_oa .GT. 1 ) ) ) THEN + CALL one_lev_coverage_check ( kbu_alloc , pressure , & + nml%record_4%max_p_tolerance_one_lev_qc , nml%record_9%max_p_tolerance_one_lev_oa ) + END IF + END IF + + ! BPR END + ! All of the following routines assume this pressure to be in hPa. pressure(:kbu_alloc) = NINT ( pressure(:kbu_alloc) * 0.01 ) @@ -347,6 +380,14 @@ ALLOCATE ( obs ( nml%record_3%max_number_of_obs ) ) ALLOCATE ( index ( nml%record_3%max_number_of_obs ) ) + + !BPR BEGIN + !Mark all time/date stamps with an obviously impossible time + !If this is done, array entries that are not filled will have + !nonsense that is not necessarily an integer which can cause problems + !when the array is read into an integer array + obs(:)%valid_time%date_char = '-9999999999999' + !BPR END ! The observation data can now be ingested, sorted, duplicates ! removed. On return, we have all of the sorted observations @@ -361,7 +402,10 @@ CALL proc_obs_sort ( dummy_filename , unit+10 , & obs , number_of_obs , nml%record_3%max_number_of_obs , nml%record_3%fatal_if_exceed_max_obs , total_dups , & index , nml%record_5%print_found_obs , nml%record_5%print_obs_files , & - kbu_alloc , pressure , & + !BPR BEGIN + !kbu_alloc , pressure , & + kbu_alloc , pressure , slp_x , t , nml%record_4%use_p_tolerance_one_lev , & + !BPR END nml%record_4%max_p_extend_t , nml%record_4%max_p_extend_w , & h , iew_alloc , jns_alloc , map_projection , current_date_8 , current_time_6 , fdda_loop ) ELSE @@ -369,7 +413,10 @@ CALL proc_obs_sort ( dummy_filename , unit+10 , & obs , number_of_obs , nml%record_3%max_number_of_obs , nml%record_3%fatal_if_exceed_max_obs , total_dups , & index , nml%record_5%print_found_obs , nml%record_5%print_obs_files , & - kbu_alloc , pressure , & + !BPR BEGIN + !kbu_alloc , pressure , & + kbu_alloc , pressure , slp_x , t , nml%record_4%use_p_tolerance_one_lev , & + !BPR END nml%record_4%max_p_extend_t , nml%record_4%max_p_extend_w , & h , iew_alloc , jns_alloc , map_projection , fdda_date_8 , fdda_time_6 , fdda_loop ) END IF @@ -382,15 +429,32 @@ total_dups , map_projection , & nml%record_4%qc_test_error_max , nml%record_4%qc_test_buddy , & nml%record_4%qc_test_vert_consistency , nml%record_4%qc_test_convective_adj , & + !BPR BEGIN + nml%record_4%qc_psfc , & + !BPR END nml%record_4%max_error_t , nml%record_4%max_error_uv , & nml%record_4%max_error_z , nml%record_4%max_error_rh , & + !BPR BEGIN + nml%record_4%max_error_dewpoint , & + !BPR END nml%record_4%max_error_p , nml%record_5%print_error_max , & nml%record_4%max_buddy_t , nml%record_4%max_buddy_uv , & nml%record_4%max_buddy_z , nml%record_4%max_buddy_rh , & - nml%record_4%max_buddy_p , nml%record_5%print_buddy , & + !BPR BEGIN + nml%record_4%max_buddy_dewpoint , & + !BPR END + !BPR BEGIN + !nml%record_4%max_buddy_p , nml%record_5%print_buddy , & + nml%record_4%max_buddy_p , & + nml%record_4%use_p_tolerance_one_lev , nml%record_4%max_p_tolerance_one_lev_qc , & + nml%record_5%print_buddy , & + !BPR END nml%record_5%print_found_obs , & nml%record_5%print_qc_vert , nml%record_5%print_qc_dry , & - pressure , current_date_8 , current_time_6 , dxd , 1. , & + !BPR BEGIN + !pressure , current_date_8 , current_time_6 , dxd , 1. , & + pressure , pres , current_date_8 , current_time_6 , dxd , 1. , & + !BPR END obs , index , nml%record_3%max_number_of_obs , & t , u , v , h , rh , slp_x , sst , tobbox , odis ) ELSE @@ -398,15 +462,32 @@ total_dups , map_projection , & nml%record_4%qc_test_error_max , nml%record_4%qc_test_buddy , & nml%record_4%qc_test_vert_consistency , nml%record_4%qc_test_convective_adj , & + !BPR BEGIN + nml%record_4%qc_psfc , & + !BPR END nml%record_4%max_error_t , nml%record_4%max_error_uv , & nml%record_4%max_error_z , nml%record_4%max_error_rh , & + !BPR BEGIN + nml%record_4%max_error_dewpoint , & + !BPR END nml%record_4%max_error_p , nml%record_5%print_error_max , & nml%record_4%max_buddy_t , nml%record_4%max_buddy_uv , & nml%record_4%max_buddy_z , nml%record_4%max_buddy_rh , & - nml%record_4%max_buddy_p , nml%record_5%print_buddy , & + !BPR BEGIN + nml%record_4%max_buddy_dewpoint , & + !BPR END + !BPR BEGIN + !nml%record_4%max_buddy_p , nml%record_5%print_buddy , & + nml%record_4%max_buddy_p , & + nml%record_4%use_p_tolerance_one_lev , nml%record_4%max_p_tolerance_one_lev_qc , & + nml%record_5%print_buddy , & + !BPR END nml%record_5%print_found_obs , & nml%record_5%print_qc_vert , nml%record_5%print_qc_dry , & - pressure , fdda_date_8 , fdda_time_6 , dxd , 1. , & + !BPR BEGIN + !pressure , fdda_date_8 , fdda_time_6 , dxd , 1. , & + pressure , pres , fdda_date_8 , fdda_time_6 , dxd , 1. , & + !BPR END obs , index , nml%record_3%max_number_of_obs , & t , u , v , h , rh , slp_x , sst , tobbox , odis ) END IF @@ -427,22 +508,34 @@ ! April 2009 - name changed from qc_out to qc_obs_raw ! We also add domain info here WRITE (tmp_filename,'("qc_obs_raw.d",i2.2,".")') nml%record_2%grid_id - CALL output_obs ( obs , 2 , trim(tmp_filename)//dt_char , number_of_obs , & - 1 , .TRUE., .TRUE., 200000, .FALSE., kbu_alloc, pressure, & - filename ) + !BPR BEGIN + use_grid_relative_winds = .TRUE. + CALL output_obs ( obs , 2 , trim(tmp_filename)//dt_char , number_of_obs , & + 1 , .TRUE., .TRUE., 200000, .FALSE., kbu_alloc, pressure , & + filename, use_grid_relative_winds ) + ! filename ) + !BPR END END IF ! With the observations properly QC'ed and stored, and the first guess ! background field already ingested, we can proceed with the objective ! analysis. - IF ( ( .NOT. nml%record_7%f4d ) .OR. & ( ( nml%record_7%f4d ) .AND. ( fdda_loop .EQ. 1 ) ) ) THEN - CALL proc_oa ( t , u , v , rh , slp_x , pressure , & + !BPR BEGIN + !Note that pressure is 2D pressure and just has 1001 for surface + !pressure whereas pres is 3D pressure and so has the full surface + !pressure field + !CALL proc_oa ( t , u , v , rh , slp_x , pressure , & + CALL proc_oa ( t , u , v , rh , slp_x , pressure , pres , & + !BPR END iew_alloc , jns_alloc , kbu_alloc , & current_date_8 , current_time_6 , fdda_loop , mqd_count , mqd_abs_min , & nml%record_3%max_number_of_obs , number_of_obs , total_dups , & map_projection , obs , dxd , lat_center , & + !BPR BEGIN + nml%record_4%use_p_tolerance_one_lev , & + !BPR END nml%record_5%print_oa , nml%record_5%print_found_obs , & nml%record_5%print_obs_files , & nml%record_7%use_first_guess , & @@ -454,7 +547,12 @@ nml%record_9%mqd_minimum_num_obs , nml%record_9%mqd_maximum_num_obs , & nml%record_9%oa_max_switch , nml%record_9%radius_influence , & nml%record_9%oa_min_switch , nml%record_9%oa_3D_option , & - nml%record_2%grid_id ) + !BPR BEGIN + !nml%record_2%grid_id ) + nml%record_2%grid_id , terrain, h, nml%record_9%scale_cressman_rh_decreases , & + nml%record_9%radius_influence_sfc_mult, nml%record_9%oa_psfc , & + nml%record_9%max_p_tolerance_one_lev_oa ) + !BPR END ! Store the final analysis back into the all_3d and all_2d arrays if we are doing ! SFC FDDA. Why? So that when we do the LAGTEM or temporal interpolation, we are @@ -464,11 +562,20 @@ CALL store_fa ( t , u , v , rh , slp_x , iew_alloc , jns_alloc , kbu_alloc , num3d , num2d , icount ) END IF ELSE - CALL proc_oa ( t , u , v , rh , slp_x , pressure , & + !BPR BEGIN + !Note that pressure is 2D pressure and just has 1001 for surface + !pressure whereas pres is 3D pressure and so has the full surface + !pressure field + !CALL proc_oa ( t , u , v , rh , slp_x , pressure , & + CALL proc_oa ( t , u , v , rh , slp_x , pressure , pres , & + !BPR END iew_alloc , jns_alloc , kbu_alloc , & fdda_date_8 , fdda_time_6 , fdda_loop , mqd_count , mqd_abs_min , & nml%record_3%max_number_of_obs , number_of_obs , total_dups , & map_projection , obs , dxd , lat_center , & + !BPR BEGIN + nml%record_4%use_p_tolerance_one_lev , & + !BPR END nml%record_5%print_oa , nml%record_5%print_found_obs , & nml%record_5%print_obs_files , & nml%record_7%use_first_guess , & @@ -480,7 +587,12 @@ nml%record_9%mqd_minimum_num_obs , nml%record_9%mqd_maximum_num_obs , & nml%record_9%oa_max_switch , nml%record_9%radius_influence , & nml%record_9%oa_min_switch , nml%record_9%oa_3D_option , & - nml%record_2%grid_id ) + !BPR BEGIN + !nml%record_2%grid_id ) + nml%record_2%grid_id , terrain, h, nml%record_9%scale_cressman_rh_decreases , & + nml%record_9%radius_influence_sfc_mult, nml%record_9%oa_psfc , & + nml%record_9%max_p_tolerance_one_lev_oa ) + !BPR END END IF END IF @@ -497,15 +609,34 @@ CALL make_date ( fdda_date_8 , fdda_time_6 , dt_char ) END IF IF ( nml%record_5%print_obs_files ) THEN + !BPR BEGIN + + !Since grid-relative speed/direction was just calculated it is now + !trustworthy, unlike in the standard version of OBSGRID + !BPR END + WRITE (tmp_filename,'("qc_obs_used.d",i2.2,".")') nml%record_2%grid_id + !BPR BEGIN + use_grid_relative_winds = .TRUE. CALL output_obs ( obs , 2 , trim(tmp_filename)//dt_char , number_of_obs , 1 , & .TRUE., .TRUE., nml%record_2%remove_data_above_qc_flag, & - nml%record_2%remove_unverified_data, kbu_alloc, pressure, & - filename ) + nml%record_2%remove_unverified_data, kbu_alloc, pressure , & + filename, use_grid_relative_winds ) + ! filename ) + !Write LITTLE_R file that is identical to qc_obs_used but with winds + !earth-relative. This is needed tor ingestion by the MET software. + WRITE (tmp_filename,'("qc_obs_used_earth_relative.d",i2.2,".")') nml%record_2%grid_id + use_grid_relative_winds = .FALSE. + CALL output_obs ( obs , 2 , trim(tmp_filename)//dt_char , number_of_obs , 1 , & + .TRUE., .TRUE., nml%record_2%remove_data_above_qc_flag, & + nml%record_2%remove_unverified_data, kbu_alloc, pressure , & + filename, use_grid_relative_winds ) + !BPR END END IF 101 CONTINUE ! come here if end of fdda loop, so we don't repeat the analysis time + IF ( ( .NOT. nml%record_7%f4d ) .OR. & ( ( nml%record_7%f4d ) .AND. ( fdda_loop .EQ. 1 ) ) ) THEN CALL proc_final_analysis ( filename , filename_out , & @@ -522,7 +653,11 @@ nml%record_4%buddy_weight , date_char , & nml%record_2%fg_filename , nml%record_9%oa_3D_option , & nml%record_7%intf4d , nml%record_7%lagtem , & - nml%record_9%oa_type , nml%record_9%oa_3D_type , nml%record_9%radius_influence ) + !BPR BEGIN + !nml%record_9%oa_type , nml%record_9%oa_3D_type , nml%record_9%radius_influence ) + nml%record_9%oa_type , nml%record_9%oa_3D_type , nml%record_9%radius_influence , & + nml%record_9%oa_psfc ) + !BPR END tobbox_ana = tobbox odis_ana = odis ELSE @@ -544,59 +679,138 @@ nml%record_4%buddy_weight , date_char , & nml%record_2%fg_filename , nml%record_9%oa_3D_option , & nml%record_7%intf4d , nml%record_7%lagtem , & - nml%record_9%oa_type , nml%record_9%oa_3D_type , nml%record_9%radius_influence ) + !BPR BEGIN + !nml%record_9%oa_type , nml%record_9%oa_3D_type , nml%record_9%radius_influence ) + nml%record_9%oa_type , nml%record_9%oa_3D_type , nml%record_9%radius_influence , & + nml%record_9%oa_psfc ) + !BPR END !nml%record_4%buddy_weight , nml%record_1%start_date ) END IF !!! Let's make sure the output to OBS_DOMAINnxx is sorted in time. IF ( dummy_filename(1:4) .NE. 'null' .AND. fdda_loop .NE. fdda_loop_max ) THEN - i_found_time = 0 - loop_sort_time : DO - DO i_num = 1, number_of_obs-1 - - IF ( obs(i_num)%valid_time%date_char > obs(i_num+1)%valid_time%date_char ) THEN - !print*,"FOUND one - ", i_num, " of ", number_of_obs, " (", i_found_time,")" - !print* ,"DATES:",obs(i_num)%valid_time%date_char ," - ", obs(i_num+1)%valid_time%date_char - obs_sort_tmp = obs(i_num) - obs(i_num) = obs(i_num+1) - obs(i_num+1) = obs_sort_tmp - i_found_time = i_found_time + 1 - END IF - - END DO - if (i_found_time == 0 ) exit loop_sort_time - if (i_found_time > 0 ) i_found_time = 0 - END DO loop_sort_time + !BPR BEGIN + !Modifications made to speed up this sorting of obs by time since this + !checking contributes notably to total run times (at least in cases with + !many observations) + !The altered code allows the time comparison to be done on integers + !rather than characters and prevents the obs from being reshuffled + !everytime two are found from being out of order. Instead, it first + !figures out the final order and then fills the derived type array based + !on the final order. + ALLOCATE ( obs_before_sort ( nml%record_3%max_number_of_obs ) ) + ALLOCATE ( date_int ( nml%record_3%max_number_of_obs ) ) + ALLOCATE ( obs_order_index ( nml%record_3%max_number_of_obs ) ) + !Save a copy of the obs before sorting them + obs_before_sort = obs + !Put all of the date strings of the obs into an integer array + !It is faster to do integer comparisons than string comparisons (at + !least for a string this length) + READ( obs(:)%valid_time%date_char,'(I14)') date_int + !Set the time in the obs array to a clearly wrong number so that if we + !fail to fill any entry in obs we will know + obs(:)%valid_time%date_char = '-9999999999999' + obs_order_index(:) = 0 + + !Choose sort type + !After each method is noted how the computation time for the method + !scales for n observations in 1) the best case, 2) the average case, 3) + !the worse case + !1 = bubble sort (original method): n, n**2, n**2 + !2 = comb sort: n, n*log(n), n**2 + sort_type = 2 + !Bubble sort + IF(sort_type.EQ.1) THEN + obs_order_index(:) = (/(i_num,i_num=1,nml%record_3%max_number_of_obs)/) + !BPR END + i_found_time = 0 + loop_sort_time : DO + DO i_num = 1, number_of_obs-1 + + !BPR BEGIN + !IF ( obs(i_num)%valid_time%date_char > obs(i_num+1)%valid_time%date_char ) THEN + IF( date_int(i_num) > date_int(i_num+1) ) THEN + !BPR END + !print*,"FOUND one - ", i_num, " of ", number_of_obs, " (", i_found_time,")" + !print* ,"DATES:",obs(i_num)%valid_time%date_char ," - ", obs(i_num+1)%valid_time%date_char + !BPR BEGIN + date_int_sort_tmp = date_int(i_num) + obs_order_sort_tmp = obs_order_index(i_num) + date_int(i_num) = date_int(i_num+1) + obs_order_index(i_num) = obs_order_index(i_num+1) + date_int(i_num+1) = date_int_sort_tmp + obs_order_index(i_num+1) = obs_order_sort_tmp + !obs_sort_tmp = obs(i_num) + !obs(i_num) = obs(i_num+1) + !obs(i_num+1) = obs_sort_tmp + !BPR END + i_found_time = i_found_time + 1 + END IF + + END DO + + if (i_found_time == 0 ) exit loop_sort_time + if (i_found_time > 0 ) i_found_time = 0 + END DO loop_sort_time + !BPR BEGIN + !Comb sort + ELSEIF( sort_type .eq. 2 ) THEN + CALL COMBSORT_INT(date_int(1:number_of_obs),obs_order_index(1:number_of_obs)) + ELSE + PRINT *,'Invalid sort_type in driver.F90 = ',sort_type + STOP + ENDIF + !Put the obs stored in obs_before_sort into obs in the order that the + !sorting code found + DO i_num = 1, number_of_obs + obs(i_num) = obs_before_sort(obs_order_index(i_num)) + ENDDO + + DEALLOCATE ( obs_order_index ) + DEALLOCATE ( date_int ) + DEALLOCATE ( obs_before_sort ) + !BPR END + ENDIF - IF ( fdda_loop.EQ.1) THEN obs_file_count = (icount-1)*2 + 1 IF ( .NOT. nml%record_7%f4d ) obs_file_count = icount WRITE (obs_nudge_file,'("OBS_DOMAIN",i1,i2.2)') nml%record_2%grid_id, obs_file_count INQUIRE ( EXIST = exist , FILE = obs_nudge_file ) - CALL output_obs ( obs , 2 , trim(obs_nudge_file), number_of_obs , 1 , & - .TRUE., exist, nml%record_2%remove_data_above_qc_flag, & - nml%record_2%remove_unverified_data, kbu_alloc, pressure, & - filename ) + !BPR BEGIN + use_grid_relative_winds = .TRUE. + CALL output_obs ( obs , 2 , trim(obs_nudge_file), number_of_obs , 1 , & + .TRUE., exist, nml%record_2%remove_data_above_qc_flag, & + nml%record_2%remove_unverified_data, kbu_alloc, pressure , & + filename, use_grid_relative_winds ) + ! filename ) + !BPR END ELSEIF ( fdda_loop.EQ.2 .AND. fdda_loop.NE.fdda_loop_max) THEN obs_file_count = (icount-1)*2 WRITE (obs_nudge_file,'("OBS_DOMAIN",i1,i2.2)') nml%record_2%grid_id, obs_file_count INQUIRE ( EXIST = exist , FILE = obs_nudge_file ) - CALL output_obs ( obs , 2 , trim(obs_nudge_file), number_of_obs , 1 , & - .TRUE., exist, nml%record_2%remove_data_above_qc_flag, & + !BPR BEGIN + use_grid_relative_winds = .TRUE. + CALL output_obs ( obs , 2 , trim(obs_nudge_file), number_of_obs , 1 , & + .TRUE., exist, nml%record_2%remove_data_above_qc_flag, & nml%record_2%remove_unverified_data, kbu_alloc, pressure, & - filename ) + filename, use_grid_relative_winds ) + ! filename ) + !BPR END ELSEIF ( fdda_loop.GT.2 .AND. fdda_loop.NE.fdda_loop_max) THEN obs_file_count = (icount-1)*2 WRITE (obs_nudge_file,'("OBS_DOMAIN",i1,i2.2)') nml%record_2%grid_id, obs_file_count - CALL output_obs ( obs , 2 , trim(obs_nudge_file), number_of_obs , 1 , & - .TRUE., .FALSE., nml%record_2%remove_data_above_qc_flag, & + !BPR BEGIN + use_grid_relative_winds = .TRUE. + CALL output_obs ( obs , 2 , trim(obs_nudge_file), number_of_obs , 1 , & + .TRUE., .FALSE., nml%record_2%remove_data_above_qc_flag, & nml%record_2%remove_unverified_data, kbu_alloc, pressure, & - filename ) + filename, use_grid_relative_winds ) + ! filename ) + !BPR END ENDIF - ! If this is the time that we have observations, we can de-allocate the space. IF ( dummy_filename(1:4) .NE. 'null' .AND. fdda_loop .NE. fdda_loop_max ) THEN Index: src/final_analysis_module.F90 =================================================================== --- src/final_analysis_module.F90 (revision 374) +++ src/final_analysis_module.F90 (working copy) @@ -11,7 +11,10 @@ latitude_x , longitude_x , latitude_d , longitude_d , & slp_x , slp_C , sst , tobbox , odis , & iew_alloc , jns_alloc , kbu_alloc , iewd , jnsd , date_char, print_analysis , & -oa_type , oa_3D_type , oa_3D_option , rad_influence ) +!BPR BEGIN +!oa_type , oa_3D_type , oa_3D_option , rad_influence ) +oa_type , oa_3D_type , oa_3D_option , rad_influence , oa_psfc ) +!BPR END ! This routine assembles the correct data fields and outputs them with the ! appropriate small header. @@ -41,6 +44,9 @@ CHARACTER *(*) :: oa_type , oa_3D_type INTEGER , DIMENSION(10) :: rad_influence INTEGER :: oa_3D_option + !BPR BEGIN + LOGICAL :: oa_psfc + !BPR END INTERFACE INCLUDE 'error.int' @@ -349,13 +355,24 @@ ALLOCATE ( d3_data1 (jns_alloc, iew_alloc, kbu_alloc) ) d3_data1 = pres - DO imet = 1,iew_alloc - DO jmet = 1,jns_alloc + !BPR BEGIN + !The surface level of the pressure field can be adjusted based on changes + !made to the sea level pressure by the objective analysis. + !However, only do this if you did not create a surface pressure objective + !analysis. + IF(.NOT.oa_psfc) THEN + !BPR END + DO imet = 1,iew_alloc + DO jmet = 1,jns_alloc d3_data1(jmet,imet,1) = pres(jmet,imet,1) + & ( pres(jmet,imet,1)/slp_C(jmet,imet) ) * & ( slp_x(jmet,imet)-slp_C(jmet,imet) ) - ENDDO - ENDDO + ENDDO + ENDDO + !BPR BEGIN + ENDIF + !BPR END + ALLOCATE ( met_em_3d (iewd-1, jnsd-1, kbu_alloc) ) ALLOCATE ( met_em_2d (iewd-1, jnsd-1 ) ) CALL unexpand3 ( d3_data1, iew_alloc, jns_alloc, kbu_alloc, dum3d, iewd, jnsd ) @@ -462,6 +479,10 @@ integer :: i , j #endif +!BPR BEGIN +integer :: varid +!BPR END + ! We need to keep track of where we are sticking the data for the FDDA option. IF ( initial_time ) THEN @@ -566,7 +587,19 @@ DO i = 1, 23 IF ( i == 1 ) THEN rcode = nf_inq_var(oa_met, i, cval, itype, idm, ishape, natt) - rcode = nf_def_var(sfc_ncid, cval, itype, idm, ishape, i) + !BPR BEGIN + if(rcode.ne.0) then + PRINT *, 'ERROR: Attempt to read variable from metoa file failed: ',NF_STRERROR(rcode) + end if + !The last argument of nf_def_var is an OUTPUT of nf_def_var + !If this variable is our loop variable than this call can change the + !value of our loop value + !rcode = nf_def_var(sfc_ncid, cval, itype, idm, ishape, i) + rcode = nf_def_var(sfc_ncid, cval, itype, idm, ishape, varid) + if(rcode.ne.0) then + PRINT *, 'ERROR: Attempt to define a variable in the surface analysis nudging file failed: ',NF_STRERROR(rcode) + end if + !BPR END ELSE ishape(1) = iwe ishape(2) = isn @@ -669,14 +702,29 @@ units = 'kg m-2' description = 'WATER EQUIVALENT SNOW DEPTH' ENDIF - rcode = nf_def_var(sfc_ncid, cval, NF_FLOAT, 3, ishape, i) - rcode = nf_put_att_int(sfc_ncid, i, 'FieldType', NF_INT, 1, 104 ) - rcode = nf_put_att_text(sfc_ncid, i, 'MemoryOrder', 3, 'XY ' ) + !BPR BEGIN + !The final argument to nf_dev_var is an OUTPUT and gives the variable id of the + !variable you are inquiring about + !If this is your loop variable you can redefine your loop variable and get stuck + !in an infinite loop! + !Therefore, instead of using the loop variable "i", we now use "varid" + !rcode = nf_def_var(sfc_ncid, cval, NF_FLOAT, 3, ishape, i) + ! rcode = nf_put_att_int(sfc_ncid, i, 'FieldType', NF_INT, 1, 104 ) + ! rcode = nf_put_att_text(sfc_ncid, i, 'MemoryOrder', 3, 'XY ' ) + ! ilen = len_trim(units) + ! rcode = nf_put_att_text(sfc_ncid, i, 'units', ilen, trim(units) ) + ! ilen = len_trim(description) + ! rcode = nf_put_att_text(sfc_ncid, i, 'description', ilen, trim(description) ) + ! rcode = nf_put_att_text(sfc_ncid, i, 'stagger', 1, stagger ) + rcode = nf_def_var(sfc_ncid, cval, NF_FLOAT, 3, ishape, varid) + rcode = nf_put_att_int(sfc_ncid, varid, 'FieldType', NF_INT, 1, 104 ) + rcode = nf_put_att_text(sfc_ncid, varid, 'MemoryOrder', 3, 'XY ' ) ilen = len_trim(units) - rcode = nf_put_att_text(sfc_ncid, i, 'units', ilen, trim(units) ) + rcode = nf_put_att_text(sfc_ncid, varid, 'units', ilen, trim(units) ) ilen = len_trim(description) - rcode = nf_put_att_text(sfc_ncid, i, 'description', ilen, trim(description) ) - rcode = nf_put_att_text(sfc_ncid, i, 'stagger', 1, stagger ) + rcode = nf_put_att_text(sfc_ncid, varid, 'description', ilen, trim(description) ) + rcode = nf_put_att_text(sfc_ncid, varid, 'stagger', 1, stagger ) + !BPR END ENDIF ENDDO rcode = nf_enddef(sfc_ncid) Index: src/get_fg_at_point_module.F90 =================================================================== --- src/get_fg_at_point_module.F90 (revision 0) +++ src/get_fg_at_point_module.F90 (working copy) @@ -0,0 +1,120 @@ +!BPR BEGIN +MODULE get_fg_at_point + +CONTAINS + +!Retrieve first guess temperature at chosen point +SUBROUTINE GET_FG_T_AT_POINT(fg_3d_t,fg_3d_h,latitude,longitude,height_msl,fg_t_at_point) + + USE map_utils + USE map_utils_helper + + IMPLICIT NONE + + REAL, DIMENSION(:,:,:),INTENT(IN) :: fg_3d_t !First guess 3D temperature + REAL, DIMENSION(:,:,:),INTENT(IN) :: fg_3d_h !First guess 3D height + REAL, INTENT(IN) :: latitude !Latitude where we want temperature + REAL, INTENT(IN) :: longitude !Longitude where we want temperature + REAL, INTENT(IN) :: height_msl !Height (m MSL) where want want temperature + REAL, INTENT(OUT) :: fg_t_at_point !First guess temperature at chosen point + + REAL :: x_location, y_location + INTEGER :: jns_alloc,iew_alloc,kbu_alloc + INTEGER, DIMENSION(4) :: nearest_x_locations, nearest_y_locations + INTEGER :: cur_point + LOGICAL :: found + REAL :: h_below, h_above + REAL :: t_below, t_above + REAL :: vinterp_weight + REAL, DIMENSION(4) :: fg_t_at_near_point + INTEGER :: k + REAL :: x_location_nonint, y_location_nonint + + INCLUDE 'missing.inc' + INCLUDE 'constants.inc' + + !Get location of ob in domain grid cell space + CALL latlon_to_ij(projd , latitude , longitude , x_location , y_location ) + + jns_alloc = size(fg_3d_t,1) + iew_alloc = size(fg_3d_t,2) + kbu_alloc = size(fg_3d_t,3) + + !If ob is outside domain then mark as missing and return + IF( (x_location .LT. 1) .OR. (x_location .GT. iew_alloc-1 ) .OR. & + (y_location .LT. 1) .OR. (y_location .GT. jns_alloc-1 ) ) THEN + fg_t_at_point = missing_r + RETURN + ENDIF + + !Create four surrounding points + nearest_x_locations(1)=floor(x_location) + nearest_y_locations(1)=floor(y_location) + nearest_x_locations(2)=floor(x_location) + nearest_y_locations(2)=ceiling(y_location) + nearest_x_locations(3)=ceiling(x_location) + nearest_y_locations(3)=floor(y_location) + nearest_x_locations(4)=ceiling(x_location) + nearest_y_locations(4)=ceiling(y_location) + + + !Vertically interpolate temperature to desired height at four surrounding grid points + !If desired point is below bottom of first guess temperature, extrapolate using standard + !atmosphere + DO cur_point = 1,4 + found = .FALSE. + find_trapping : DO k = 2 , kbu_alloc-1 + IF ( ( height_msl .GE. fg_3d_h(nearest_y_locations(cur_point), & + nearest_x_locations(cur_point),k ) ) .AND. & + ( height_msl .LE. fg_3d_h(nearest_y_locations(cur_point), & + nearest_x_locations(cur_point),k+1 ) ) ) THEN + found = .TRUE. + h_below = fg_3d_h(nearest_y_locations(cur_point),nearest_x_locations(cur_point),k ) + h_above = fg_3d_h(nearest_y_locations(cur_point),nearest_x_locations(cur_point),k+1 ) + t_below = fg_3d_t(nearest_y_locations(cur_point),nearest_x_locations(cur_point),k ) + t_above = fg_3d_t(nearest_y_locations(cur_point),nearest_x_locations(cur_point),k+1 ) + EXIT find_trapping + END IF + END DO find_trapping + !If we have first guess levels sandwiching the desired height + IF( found ) THEN + !Vertically interpolate temperature to desired height + !height_msl + vinterp_weight = (h_above - height_msl ) / ( h_above - h_below ) + fg_t_at_near_point(cur_point) = vinterp_weight * t_below + (1 - vinterp_weight) * t_above + + !The desired height is below the lowest non-surface first-guess level + ELSEIF( height_msl .LT. fg_3d_h(nearest_y_locations(cur_point),nearest_x_locations(cur_point),2 ) ) THEN + fg_t_at_near_point(cur_point) = fg_3d_t(nearest_y_locations(cur_point),nearest_x_locations(cur_point),2 ) - & + ( dTdz_std_atm_lt_11000m * & + ( fg_3d_h(nearest_y_locations(cur_point),nearest_x_locations(cur_point),2 ) - height_msl ) ) + + !The desired height is above the top first-guess level + ELSEIF( height_msl .GT. & + fg_3d_h(nearest_y_locations(cur_point),nearest_x_locations(cur_point),kbu_alloc ) ) THEN + fg_t_at_point = missing_r + RETURN + !We should never get to this branch so error out + ELSE + PRINT *,'ERROR: GET_FG_T_AT_POINT: Entered branch that should never be entered: height_msl = ',height_msl + STOP + + ENDIF + + ENDDO + + !Horizontally interpolate to desired location + x_location_nonint = x_location - floor(x_location) + y_location_nonint = y_location - floor(y_location) + + fg_t_at_point = ( 1.-x_location_nonint ) * ( ( 1.-y_location_nonint ) * fg_t_at_near_point(1) & + + y_location_nonint * fg_t_at_near_point(2) ) & + + x_location_nonint * ( ( 1.-y_location_nonint ) * fg_t_at_near_point(3) & + + y_location_nonint * fg_t_at_near_point(4) ) + + +END SUBROUTINE GET_FG_T_AT_POINT + +!------------------------------------------------------------------------------ +END MODULE get_fg_at_point +!BPR END Property changes on: src/get_fg_at_point_module.F90 ___________________________________________________________________ Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: src/main.F90 =================================================================== --- src/main.F90 (revision 374) +++ src/main.F90 (working copy) @@ -14,6 +14,10 @@ INCLUDE 'error.inc' INCLUDE 'big_header.inc' + !BPR BEGIN + INCLUDE 'netcdf.inc' + INTEGER rcode, OLD_FORMAT + !BPR END TYPE ( all_nml ) :: nml ! all NAMELIST information from all NAMELIST records CHARACTER ( LEN = 132 ) :: filename, filename_out @@ -81,12 +85,24 @@ WRITE ( UNIT = * , FMT = '(" ")' ) WRITE ( UNIT = * , FMT = '("################################ ")' ) WRITE ( UNIT = * , FMT = '(" WRF OBSGRID ")' ) - WRITE ( UNIT = * , FMT = '(" Version 3.7.0 ")' ) - WRITE ( UNIT = * , FMT = '(" July 2014 ")' ) + WRITE ( UNIT = * , FMT = '(" Version 3.7.0 ")' ) + WRITE ( UNIT = * , FMT = '(" July 2014 ")' ) !!WRITE ( UNIT = * , FMT = '(" pre-release - 02/-9/10 ")' ) WRITE ( UNIT = * , FMT = '("################################ ")' ) WRITE ( UNIT = * , FMT = '(" ")' ) + !BPR BEGIN + !Set the default format for netCDF files that are written to deal with larger + !file sizes + rcode = NF_SET_DEFAULT_FORMAT(nf_format_64bit, OLD_FORMAT) + IF (rcode .NE. NF_NOERR) THEN + PRINT *, 'Attempt to set default netCDF format failed with: ',NF_STRERROR(rcode) + STOP 'Stop: Error setting default netCDF format' + ELSE + PRINT *,'Changing default netCDF format from ',OLD_FORMAT,' to nf_format_64bit' + ENDIF + !BPR END + ! Read in the NAMELIST information. This file is connected to the given ! compile-time name. All error processing on the NAMELIST data are ! handled by this routine (all simple errors that allow consistency Index: src/map_utils_module.F90 =================================================================== --- src/map_utils_module.F90 (revision 374) +++ src/map_utils_module.F90 (working copy) @@ -563,7 +563,8 @@ ! Compute cone factor CALL lc_cone(proj%truelat1, proj%truelat2, proj%cone) - PRINT '(A,F8.6)', 'Computed cone factor: ', proj%cone + !mchen PRINT '(A,F8.6)', 'Computed cone factor: ', proj%cone + PRINT '(A,F8.3)', 'Computed cone factor: ', proj%cone ! Compute longitude differences and ensure we stay out of the ! forbidden "cut zone" deltalon1 = proj%lon1 - proj%stdlon Index: src/module_report.F90 =================================================================== --- src/module_report.F90 (revision 374) +++ src/module_report.F90 (working copy) @@ -38,10 +38,14 @@ end type report - character(len=84), parameter :: rpt_format =& + !BPR BEGIN + !rpt_format was dimensioned one character too small + !character(len=84), parameter :: rpt_format =& + character(len=85), parameter :: rpt_format =& ' ( 2f20.5 , 2a40 , ' & // ' 2a40 , 1f20.5 , 5i10 , 3L10 , ' & // ' 2i10 , a20 , 13( f13.5 , i7 ) ) ' + !BPR END character(len=22), parameter :: meas_format =& ' ( 10( f13.5 , i7 ) ) ' Index: src/module_skewt.F90 =================================================================== --- src/module_skewt.F90 (revision 374) +++ src/module_skewt.F90 (working copy) @@ -705,9 +705,15 @@ END FUNCTION SKEWT_TSA FUNCTION SKEWT_W(T, P) -! W(GRAMS WATER VAPOR/KILOGRAM DRY AIR ), P(MILLIBAR ) +!BPR BEGIN +!! W(GRAMS WATER VAPOR/KILOGRAM DRY AIR ), P(MILLIBAR ) +! SKEWT_W(GRAMS WATER VAPOR/KILOGRAM DRY AIR ), P(MILLIBAR ) +!BPR END IF (T .GE. 999.) THEN - W = 0.0 +!BPR BEGIN +! W = 0.0 + SKEWT_W = 0.0 +!BPR END ELSE X = SKEWT_ESAT(T) SKEWT_W = 621.97 * X / (P - X) Index: src/namelist.common =================================================================== --- src/namelist.common (revision 374) +++ src/namelist.common (working copy) @@ -34,21 +34,35 @@ COMMON /record4_nmll/ qc_test_error_max , & qc_test_buddy , & qc_test_vert_consistency , & - qc_test_convective_adj +!BPR BEGIN +! qc_test_convective_adj + qc_test_convective_adj , & + qc_psfc , & + use_p_tolerance_one_lev +!BPR END COMMON /record4_nmlr/ max_error_t , & max_error_uv , & max_error_z , & max_error_rh , & +!BPR BEGIN + max_error_dewpoint , & +!BPR END max_error_p , & max_buddy_t , & max_buddy_uv , & max_buddy_z , & max_buddy_rh , & +!BPR BEGIN + max_buddy_dewpoint , & +!BPR END max_buddy_p , & buddy_weight , & max_p_extend_t , & - max_p_extend_w + max_p_extend_w +!BPR BEGIN + COMMON /record4_nmli/ max_p_tolerance_one_lev_qc +!BPR END ! Record 5 NAMELIST variables. @@ -91,4 +105,14 @@ COMMON /record9_nmll/ oa_min_switch , & oa_max_switch , & - oa_3D_option +!BPR BEGIN +! oa_3D_option & + oa_3D_option , & + oa_psfc , & + scale_cressman_rh_decreases +!BPR END + +!BPR BEGIN + COMMON /record9_nmlr/ radius_influence_sfc_mult + COMMON /record9_nmli/ max_p_tolerance_one_lev_oa +!BPR END Index: src/namelist.inc =================================================================== --- src/namelist.inc (revision 374) +++ src/namelist.inc (working copy) @@ -32,22 +32,37 @@ LOGICAL :: qc_test_error_max , & qc_test_buddy , & qc_test_vert_consistency , & - qc_test_convective_adj + ! BPR BEGIN + ! qc_test_convective_adj + qc_test_convective_adj , & + qc_psfc + ! BPR END REAL :: max_error_t , & max_error_uv , & max_error_z , & max_error_rh , & + ! BPR BEGIN + max_error_dewpoint , & + ! BPR END max_error_p , & max_buddy_t , & max_buddy_uv , & max_buddy_z , & max_buddy_rh , & + ! BPR BEGIN + max_buddy_dewpoint , & + ! BPR END max_buddy_p , & buddy_weight , & max_p_extend_t , & max_p_extend_w +!BPR BEGIN + LOGICAL :: use_p_tolerance_one_lev + INTEGER :: max_p_tolerance_one_lev_qc +!BPR END + ! Record 5 NAMELIST variables. LOGICAL :: print_found_obs , & @@ -90,3 +105,9 @@ LOGICAL :: oa_min_switch , & oa_max_switch INTEGER :: oa_3D_option + !BPR BEGIN + LOGICAL :: oa_psfc + LOGICAL :: scale_cressman_rh_decreases + REAL :: radius_influence_sfc_mult + INTEGER :: max_p_tolerance_one_lev_oa + !BPR END Index: src/namelist.nml =================================================================== --- src/namelist.nml (revision 374) +++ src/namelist.nml (working copy) @@ -31,20 +31,35 @@ qc_test_buddy , & qc_test_vert_consistency , & qc_test_convective_adj , & +! BPR BEGIN + qc_psfc , & +! BPR END max_error_t , & max_error_uv , & max_error_z , & max_error_rh , & +! BPR BEGIN + max_error_dewpoint , & +! BPR END max_error_p , & max_buddy_t , & max_buddy_uv , & max_buddy_z , & max_buddy_rh , & +! BPR BEGIN + max_buddy_dewpoint , & +! BPR END max_buddy_p , & buddy_weight , & max_p_extend_t , & - max_p_extend_w +! BPR BEGIN +! max_p_extend_w + max_p_extend_w , & + use_p_tolerance_one_lev , & + max_p_tolerance_one_lev_qc +! BPR END + ! Record 5 NAMELIST variables. NAMELIST /record5/ print_found_obs , & @@ -84,4 +99,10 @@ oa_min_switch , & oa_max_switch , & oa_3D_option , & - oa_3D_type + oa_3D_type , & +!BPR BEGIN + oa_psfc , & + scale_cressman_rh_decreases , & + radius_influence_sfc_mult , & + max_p_tolerance_one_lev_oa +!BPR END Index: src/namelist_module.F90 =================================================================== --- src/namelist_module.F90 (revision 374) +++ src/namelist_module.F90 (working copy) @@ -37,21 +37,39 @@ LOGICAL :: qc_test_error_max , & ! perform error-max check qc_test_buddy , & ! perform buddy check qc_test_vert_consistency , & ! T, u, v vertical consistency - qc_test_convective_adj ! lapse rate tests +!BPR BEGIN +! qc_test_convective_adj ! lapse rate tests + qc_test_convective_adj , & ! lapse rate tests + qc_psfc ! quality control surface pressure +!BPR END REAL :: max_error_t , & max_error_uv , & max_error_z , & max_error_rh , & +!BPR BEGIN + max_error_dewpoint , & +!BPR END max_error_p , & max_buddy_t , & max_buddy_uv , & max_buddy_z , & max_buddy_rh , & +!BPR BEGIN + max_buddy_dewpoint , & +!BPR END max_buddy_p , & buddy_weight , & max_p_extend_t , & max_p_extend_w +!BPR BEGIN + LOGICAL :: use_p_tolerance_one_lev ! T/F Allow single level obs vertically placed within + ! a pressure tolerance of the pressure level for which + ! obs are being sought to be used + INTEGER :: max_p_tolerance_one_lev_qc ! Pressure tolerance between ob and the first guess + ! pressure level that can be used to QC it (Pa) + ! use_p_tolerance_one_lev must be set to true to use this +!BPR END END TYPE nml_record_4 TYPE nml_record_5 @@ -97,6 +115,21 @@ LOGICAL :: oa_min_switch , & ! if MQD has less than number of obs, T/F use Cressman oa_max_switch ! if MQD has more than number of obs, T/F use Cressman INTEGER :: oa_3D_option +!BPR BEGIN + LOGICAL :: oa_psfc ! Should surface pressure be objectively analyzed + LOGICAL :: scale_cressman_rh_decreases ! Should Cressman-caused decreases to RH at locations + ! with lower RH than at the ob location be scaled based + ! on the relationship between RH at the ob location and + ! at the application location + REAL :: radius_influence_sfc_mult ! To get the radius of influence used for surface obs + ! what should the radius of influence used by non-surface + ! obs be multiplied by + INTEGER :: max_p_tolerance_one_lev_oa ! Maximum pressure difference allowed between a single + ! level ob and the pressure level it is analyzed on + ! User must enable use_p_tolerance_one_lev for this to + ! have an effect and must be less than max_p_tolerance_one_lev_qc + +!BPR END END TYPE nml_record_9 TYPE all_nml @@ -286,6 +319,73 @@ call error_handler ( error_number , error_message , fatal , listing ) END IF test_405 + !BPR BEGIN + ! 406: Perform quality control on surface pressure? + ! Tell user whether or not they are quality controlling surface pressure + ! If they are not, but they are objectively analyzing that field, then throw + ! a fatal error given the danger of objectively analyzing non-QC'd obs + + error_message = ' ' + error_number = 00013406 + error_message(1:31) = 'check_namelist ' + listing = .false. + + test_406 : IF ( .NOT. nml%record_4%qc_psfc ) THEN + IF ( nml%record_9%oa_psfc ) THEN + fatal = .true. + error_message(32:) = ' Surface pressure is NOT being QCed but is being objectively analyzed.' + ELSE + fatal = .false. + error_message(32:) = ' Quality control checks will NOT be performed on surface pressure.' + END IF + ELSE + fatal = .false. + error_message(32:) = ' Quality control checks WILL be performed on surface pressure.' + END IF test_406 + call error_handler ( error_number , error_message , fatal , listing ) + + ! 407: Allow pressure tolerance when doing QC of single level above surface obs? + + error_message = ' ' + error_number = 00013407 + error_message(1:31) = 'check_namelist ' + listing = .false. + fatal = .false. + + test_407 : IF ( nml%record_4%use_p_tolerance_one_lev ) THEN + WRITE(error_message(32:),'(A,I9,A)') ' Single-level above-surface obs will be QCed & + &against the nearest pressure level in the analysis as long as the pressure difference & + &is no more than',& + nml%record_4%max_p_tolerance_one_lev_qc, & + ' Pa. No obs will be extended to the nearest pressure level.' + ELSE + error_message(32:) = ' Single-level above-surface obs will generally only be QCed if they fall& + & on an analysis pressure level. However, will attempt to adjust obs marked as FM-88 SATOB or FM-97 AIREP& + & to the nearest pressure level.' + END IF test_407 + call error_handler ( error_number , error_message , fatal , listing ) + + ! 408: Make user aware that a pressure tolerance of 1 hPa for using a single-level above-surface + ! ob on the objective analysis at a pressure level, is like having no + ! tolerance at all. + + test_408 : IF ( ( nml%record_4%use_p_tolerance_one_lev ) .AND. & + ( nml%record_4%max_p_tolerance_one_lev_qc .EQ. 1 ) ) THEN + error_message = ' ' + error_number = 00013408 + error_message(1:31) = 'check_namelist ' + listing = .false. + fatal = .false. + + WRITE(error_message(32:),'(A)') ' Pressure tolerance for single-level above-surface obs being QCed & + &against an analysis is 1 hPa. This is equivalent to having no pressure tolerance!' + + call error_handler ( error_number , error_message , fatal , listing ) + END IF test_408 + + + !BPR END + ! 501: Print out message for each observation found? test_501 : IF ( .NOT. nml%record_5%print_found_obs ) THEN @@ -537,8 +637,14 @@ ! 901: Which objective analysis technique is to be used? + !BPR BEGIN + !Allow "None" option for cases in which we only want to QC data + !test_901 : IF ( ( nml%record_9%oa_type .NE. 'MQD' ) .AND. & + ! ( nml%record_9%oa_type .NE. 'Cressman' ) ) THEN test_901 : IF ( ( nml%record_9%oa_type .NE. 'MQD' ) .AND. & - ( nml%record_9%oa_type .NE. 'Cressman' ) ) THEN + ( nml%record_9%oa_type .NE. 'Cressman' ) .AND. & + ( nml%record_9%oa_type .NE. 'None' ) ) THEN + !BPR END nml%record_9%oa_type = 'Cressman' END IF test_901 @@ -654,6 +760,89 @@ call error_handler ( error_number , error_message , fatal , listing ) END IF test_906 + + !BPR BEGIN + ! 907: What is the minimum number of observations required to still do MQD? + + error_message = ' ' + error_number = 00013907 + error_message(1:31) = 'check_namelist ' + fatal = .false. + listing = .false. + + test_907 : IF ( nml%record_9%oa_psfc ) THEN + error_message(32:) = ' Will objectively analyze surface pressure [psfc will NOT be affected by SLP analysis].' + ELSE + error_message(32:) = ' Will NOT objectively analyze surface pressure [psfc will be affected by SLP analysis].' + END IF test_907 + call error_handler ( error_number , error_message , fatal , listing ) + + ! 908: Allow pressure tolerance when doing OA of single level above + ! surface obs? + + error_message = ' ' + error_number = 00013908 + error_message(1:31) = 'check_namelist ' + listing = .false. + fatal = .false. + + test_908 : IF ( nml%record_4%use_p_tolerance_one_lev ) THEN + WRITE(error_message(32:),'(A,I9,A)') ' Single-level above-surface obs will be included in & + &the objective analysis on the nearest pressure level as long as the pressure difference & + &is no more than',& + nml%record_9%max_p_tolerance_one_lev_oa, & + ' Pa. No obs will be extended to the nearest pressure level.' + ELSE + error_message(32:) = ' Single-level above-surface obs will generally only be included in & + &the objective analysis if they fall on an analysis pressure level. However, will attempt & + &to adjust obs marked as FM-88 SATOB or FM-97 AIREP to the nearest pressure level.' + END IF test_908 + call error_handler ( error_number , error_message , fatal , listing ) + + ! 909: Make sure pressure tolerance allowed between single-level above-surface obs and the + ! pressure level they are analyzed on is less than or equal to the pressure tolerance used to + ! do the QC + + test_909 : IF ( ( nml%record_4%use_p_tolerance_one_lev ) .AND. & + ( nml%record_9%max_p_tolerance_one_lev_oa .GT. nml%record_4%max_p_tolerance_one_lev_qc ) ) THEN + error_message = ' ' + error_number = 00013909 + error_message(1:31) = 'check_namelist ' + listing = .false. + fatal = .true. + + WRITE(error_message(32:),'(A,I9,A,I9,A)') ' Pressure tolerance for single-level above-surface obs being included in & + &the objective analysis is greater than the pressure tolerance for QCing the data. For QC tolerance is: ',& + nml%record_4%max_p_tolerance_one_lev_qc, & + ' Pa, for OA tolerance is: ',& + nml%record_9%max_p_tolerance_one_lev_oa,& + ' Pa' + + call error_handler ( error_number , error_message , fatal , listing ) + END IF test_909 + + ! 910: Make user aware that a pressure tolerance of 1 hPa for using a single-level above-surface + ! ob on the objective analysis at a pressure level, is like having no + ! tolerance at all. + + test_910 : IF ( ( nml%record_4%use_p_tolerance_one_lev ) .AND. & + ( nml%record_9%max_p_tolerance_one_lev_oa .EQ. 1 ) ) THEN + error_message = ' ' + error_number = 00013910 + error_message(1:31) = 'check_namelist ' + listing = .false. + fatal = .false. + + WRITE(error_message(32:),'(A)') ' Pressure tolerance for single-level above-surface obs being included in & + &the objective analysis is 1 hPa. This is equivalent to having no pressure tolerance!' + + call error_handler ( error_number , error_message , fatal , listing ) + END IF test_910 + !BPR END + + + + END SUBROUTINE check_namelist !------------------------------------------------------------------------------ @@ -745,11 +934,17 @@ nml%record_4%qc_test_buddy = qc_test_buddy nml%record_4%qc_test_vert_consistency = qc_test_vert_consistency nml%record_4%qc_test_convective_adj = qc_test_convective_adj + !BPR BEGIN + nml%record_4%qc_psfc = qc_psfc + !BPR END nml%record_4%max_error_t = max_error_t nml%record_4%max_error_uv = max_error_uv nml%record_4%max_error_z = max_error_z nml%record_4%max_error_rh = max_error_rh + !BPR BEGIN + nml%record_4%max_error_dewpoint = max_error_dewpoint + !BPR END nml%record_4%max_error_p = max_error_p nml%record_4%max_buddy_t = max_buddy_t @@ -756,6 +951,9 @@ nml%record_4%max_buddy_uv = max_buddy_uv nml%record_4%max_buddy_z = max_buddy_z nml%record_4%max_buddy_rh = max_buddy_rh + !BPR BEGIN + nml%record_4%max_buddy_dewpoint = max_buddy_dewpoint + !BPR END nml%record_4%max_buddy_p = max_buddy_p nml%record_4%buddy_weight = buddy_weight @@ -762,6 +960,10 @@ nml%record_4%max_p_extend_t = max_p_extend_t nml%record_4%max_p_extend_w = max_p_extend_w + !BPR BEGIN + nml%record_4%use_p_tolerance_one_lev = use_p_tolerance_one_lev + nml%record_4%max_p_tolerance_one_lev_qc = max_p_tolerance_one_lev_qc + !BPR END ! Record 5 NAMELIST values: @@ -803,9 +1005,94 @@ nml%record_9%oa_min_switch = oa_min_switch nml%record_9%oa_max_switch = oa_max_switch nml%record_9%oa_3D_option = oa_3D_option + !BPR BEGIN + nml%record_9%oa_psfc = oa_psfc + nml%record_9%scale_cressman_rh_decreases = scale_cressman_rh_decreases + nml%record_9%radius_influence_sfc_mult = radius_influence_sfc_mult + nml%record_9%max_p_tolerance_one_lev_oa = max_p_tolerance_one_lev_oa + !BPR END END SUBROUTINE store_namelist !------------------------------------------------------------------------------ +!Find pressure ranges where observations cannot be QC'ed against the first-guess +!field or cannot be included in the objective analysis +SUBROUTINE one_lev_coverage_check ( kbu_alloc , pressure , & + max_p_tolerance_one_lev_qc, max_p_tolerance_one_lev_oa ) + + IMPLICIT NONE + + !Number of vertical levels in first guess field + INTEGER, INTENT(IN) :: kbu_alloc + !Pressure of each vertical level (Pa). The first level, the surface just + !has 100100 since in reality the pressure of the surface varies across the + !domain + REAL , INTENT(IN) , DIMENSION(kbu_alloc) :: pressure + !Maximum pressure difference (Pa) between ob and the pressure level it can + !be QCed against + INTEGER , INTENT(IN) :: max_p_tolerance_one_lev_qc + !Maximum pressure difference (Pa) between ob and the pressure level on which + !it can be used in an objective analysis + INTEGER , INTENT(IN) :: max_p_tolerance_one_lev_oa + + INTEGER :: curk + INTEGER , DIMENSION(kbu_alloc) :: top_of_p_range_qc, bot_of_p_range_qc + INTEGER , DIMENSION(kbu_alloc) :: top_of_p_range_oa, bot_of_p_range_oa + + !Find the top and bottom of pressure ranges around each pressure level in + !the first-guess field + top_of_p_range_qc(1) = 0 + bot_of_p_range_qc(1) = 0 + DO curk = 2,kbu_alloc + top_of_p_range_qc(curk) = pressure(curk) - max_p_tolerance_one_lev_qc + bot_of_p_range_qc(curk) = pressure(curk) + max_p_tolerance_one_lev_qc + top_of_p_range_oa(curk) = pressure(curk) - max_p_tolerance_one_lev_oa + bot_of_p_range_oa(curk) = pressure(curk) + max_p_tolerance_one_lev_oa + ENDDO + + IF( max_p_tolerance_one_lev_qc .GT. 1 ) THEN + PRINT *,' ' + PRINT *,'*******************************************************************************' + PRINT *,'Quality control of single-level above-surface observations can use first-guess ' + PRINT *,' fields within ',max_p_tolerance_one_lev_qc,' Pa of observation' + PRINT *,'Single-level above-surface obs will not be quality controlled against ' + PRINT *,' the first guess field if their pressure falls in the following ranges: ' + WRITE(*,'(A3,I6,A3)') '<= ',top_of_p_range_qc(kbu_alloc),' Pa' + DO curk = kbu_alloc,3,-1 + !Use ".GE." here because pressure difference must be LESS than the + !tolerance, not equal to the tolerance + IF( top_of_p_range_qc(curk-1) .GE. bot_of_p_range_qc(curk) ) THEN + WRITE(*,'(I6,A3,I6,A3)') top_of_p_range_qc(curk-1),' - ',bot_of_p_range_qc(curk),' Pa' + END IF + ENDDO + WRITE(*,'(A3,I6,A3)') '>= ',bot_of_p_range_qc(2),' Pa' + PRINT *,'*******************************************************************************' + PRINT *,' ' + END IF + + IF( max_p_tolerance_one_lev_oa .GT. 1 ) THEN + PRINT *,'*******************************************************************************' + PRINT *,' ' + PRINT *,'Objective analysis can use single-level above-surface observations within ' + PRINT *,max_p_tolerance_one_lev_oa,' Pa of observation' + PRINT *,'Single-level above-surface obs will not be included in the objective analysis' + PRINT *,' at any level if their pressure falls in the following ranges: ' + WRITE(*,'(A3,I6,A3)') '<= ',top_of_p_range_oa(kbu_alloc),' Pa' + DO curk = kbu_alloc,3,-1 + !Use ".GE." here because pressure difference must be LESS than the + !tolerance, not equal to the tolerance + IF( top_of_p_range_oa(curk-1) .GE. bot_of_p_range_oa(curk) ) THEN + WRITE(*,'(I6,A3,I6,A3)') top_of_p_range_oa(curk-1),' - ',bot_of_p_range_oa(curk),' Pa' + END IF + ENDDO + WRITE(*,'(A3,I6,A3)') '>= ',bot_of_p_range_oa(2),' Pa' + PRINT *,'*******************************************************************************' + PRINT *,' ' + END IF + + +END SUBROUTINE one_lev_coverage_check +!------------------------------------------------------------------------------ + END MODULE namelist Index: src/oa_module.F90 =================================================================== --- src/oa_module.F90 (revision 374) +++ src/oa_module.F90 (working copy) @@ -26,11 +26,17 @@ gridded , iew , jns , & crsdot , name , radius_influence , & dxd , u_banana , v_banana , pressure , passes , smooth_type , lat_center , & -use_first_guess ) +!BPR BEGIN +!use_first_guess ) +use_first_guess , fg_value, scale_cressman_rh_decreases ) +!BPR END ! arguments: ! obs_value: difference between 1st guess and obs ! ( need any duplicated obs removed from this array ) +!BPR BEGIN + ! ACTUALLY obs_value appears to be a value NOT a difference +!BPR END ! obs: difference of obs_value and the interpolated first_guess value ! numobs: number of station obs ! xob,yob: x, y locations of station obs @@ -42,6 +48,9 @@ ! v_banana first guess field of V for banana scheme ! pressure pressure of this level (Pa, 1001 signifies surface) ! use_first_guess T/F to include first-guess in the objective analysis +!BPR BEGIN + ! fg_value: Model value at the ob location +!BPR END IMPLICIT NONE @@ -48,6 +57,9 @@ INTEGER, INTENT ( IN ) :: numobs INTEGER, INTENT ( IN ) :: iew, jns REAL, INTENT ( IN ), DIMENSION ( numobs ) :: obs_value , obs, xob, yob +!BPR BEGIN + REAL, INTENT ( IN ), DIMENSION ( numobs ) :: fg_value +!BPR END REAL, INTENT ( INOUT ), DIMENSION( iew , jns ) :: gridded INTEGER, INTENT ( IN ) :: radius_influence CHARACTER (8), INTENT ( IN ) :: name @@ -70,8 +82,16 @@ jns_min , & iew_max , & jns_max - CHARACTER (LEN=8) :: choice + !BPR BEGIN + !CHARACTER (LEN=8) :: choice + REAL :: scale_factor + REAL :: real_crsdot_divide_2, real_i, real_j + REAL :: distance_squared + INTEGER :: radius_influence_squared + LOGICAL :: scale_cressman_rh_decreases + !BPR END + INCLUDE 'error.inc' INTERFACE @@ -85,6 +105,12 @@ denominator = 0 perturbation = 0 + !BPR BEGIN + !Calculate this here so we do not have to do it many, many times below + real_crsdot_divide_2 = REAL(crsdot)/2. + radius_influence_squared = radius_influence * radius_influence + !BPR END + ! Loop over each observation. obs_loop : DO n = 1 , numobs @@ -92,16 +118,24 @@ ! Every observation has a first guess wind field location that will be specified ! as an objective analysis that is circular, elliptical or banana-shaped. - choice = 'undefine' + !BPR BEGIN + !choice = 'undefine' + !BPR END ! To use the observation, it must be far enough inside our domain ! (using the cross/dot configuration for the staggered variables). If ! it is not, we simply zoom back to the top of this loop. - IF ( ( xob(n) .LE. 1 + REAL(crsdot)/2. ) .OR. & - ( yob(n) .LE. 1 + REAL(crsdot)/2. ) .OR. & - ( xob(n) .GE. iew - REAL(crsdot)/2. ) .OR. & - ( yob(n) .GE. jns - REAL(crsdot)/2. ) ) THEN + !BPR BEGIN + !IF ( ( xob(n) .LE. 1 + REAL(crsdot)/2. ) .OR. & + ! ( yob(n) .LE. 1 + REAL(crsdot)/2. ) .OR. & + ! ( xob(n) .GE. iew - REAL(crsdot)/2. ) .OR. & + ! ( yob(n) .GE. jns - REAL(crsdot)/2. ) ) THEN + IF ( ( xob(n) .LE. 1 + real_crsdot_divide_2 ) .OR. & + ( yob(n) .LE. 1 + real_crsdot_divide_2 ) .OR. & + ( xob(n) .GE. iew - real_crsdot_divide_2 ) .OR. & + ( yob(n) .GE. jns - real_crsdot_divide_2 ) ) THEN + !BPR END CYCLE obs_loop END IF @@ -109,47 +143,49 @@ ! some points inside this window the observation will have zero ! influence, but this allows a direct method rather than a search. - iew_min = MAX ( NINT ( (xob(n)-REAL(crsdot)/2.) ) - 4*radius_influence -1 , 1 ) - jns_min = MAX ( NINT ( (yob(n)-REAL(crsdot)/2.) ) - 4*radius_influence -1 , 1 ) - iew_max = MIN ( NINT ( (xob(n)-REAL(crsdot)/2.) ) + 4*radius_influence +1 , iew - crsdot ) - jns_max = MIN ( NINT ( (yob(n)-REAL(crsdot)/2.) ) + 4*radius_influence +1 , jns - crsdot ) + !BPR BEGIN + !iew_min = MAX ( NINT ( (xob(n)-REAL(crsdot)/2.) ) - 4*radius_influence -1 , 1 ) + !jns_min = MAX ( NINT ( (yob(n)-REAL(crsdot)/2.) ) - 4*radius_influence -1 , 1 ) + !iew_max = MIN ( NINT ( (xob(n)-REAL(crsdot)/2.) ) + 4*radius_influence +1 , iew - crsdot ) + !jns_max = MIN ( NINT ( (yob(n)-REAL(crsdot)/2.) ) + 4*radius_influence +1 , jns - crsdot ) + iew_min = MAX ( NINT ( (xob(n)-real_crsdot_divide_2) ) - 4*radius_influence -1 , 1 ) + jns_min = MAX ( NINT ( (yob(n)-real_crsdot_divide_2) ) - 4*radius_influence -1 , 1 ) + iew_max = MIN ( NINT ( (xob(n)-real_crsdot_divide_2) ) + 4*radius_influence +1 , iew - crsdot ) + jns_max = MIN ( NINT ( (yob(n)-real_crsdot_divide_2) ) + 4*radius_influence +1 , jns - crsdot ) + !BPR END ! For the grid points in the window, compute this observation's ! influence. - DO j = jns_min , jns_max - DO i = iew_min , iew_max + !BPR BEGIN + !Significantly edited code for compuatational efficiency. The exact same + !algorithms are used, they are just restructured to speed things up since + !we want to minimize the code called on the innermost loop where we are + !looping over the x and y direction for every observation. - ! Distance from each of the windowed grid points to the observation. The distance is in - ! grid point units. If this is a horizontal wind field or a moisture field, we can use the - ! banana scheme, which elongates the acceptable distance from an observation and also is - ! able to handle first-guess wind curvature. + ! Distance from each of the windowed grid points to the observation. The distance is in + ! grid point units. If this is a horizontal wind field or a moisture field, we can use the + ! banana scheme, which elongates the acceptable distance from an observation and also is + ! able to handle first-guess wind curvature. - IF ( ( name(1:8) .EQ. 'UU ' ) .OR. & - ( name(1:8) .EQ. 'VV ' ) .OR. & - ( name(1:8) .EQ. 'RH ' ) ) THEN - CALL dist ( i , j , iew , jns , crsdot , name(1:8) , obs_value(n) , xob(n) , yob(n) , & - dxd , pressure , u_banana , v_banana , radius_influence , & - lat_center , choice , distance ) - - ELSE - distance = SQRT ( ( ( xob(n) - REAL(crsdot)/2. - REAL(i) ) * ( xob(n) - REAL(crsdot)/2. - REAL(i) ) ) + & - ( ( yob(n) - REAL(crsdot)/2. - REAL(j) ) * ( yob(n) - REAL(crsdot)/2. - REAL(j) ) ) ) + IF ( ( name(1:8) .EQ. 'UU ' ) .OR. & + ( name(1:8) .EQ. 'VV ' ) .OR. & + ( name(1:8) .EQ. 'RH ' ) ) THEN + CALL dist_uu_vv_rh ( iew_min, iew_max, jns_min, jns_max, iew , jns , real_crsdot_divide_2, name(1:8) , & + obs(n), xob(n) , yob(n) , & + dxd , pressure , u_banana , v_banana , radius_influence , & + lat_center , gridded, fg_value(n), use_first_guess, numerator , denominator, & !BPR BEGIN + scale_cressman_rh_decreases ) +! lat_center , gridded, fg_value(n), use_first_guess, numerator , denominator ) !BPR END + ELSE - END IF + CALL dist_not_uu_vv_rh ( iew_min, iew_max, jns_min, jns_max, real_crsdot_divide_2, name(1:8) , & + obs(n), xob(n) , yob(n) , radius_influence , numerator, denominator ) - ! If the observation is within a grid point's radius of influence, - ! sum its effect. + END IF + !BPR END - IF ( distance .LT. radius_influence ) THEN - weight = ( radius_influence**2 - distance**2 ) / ( radius_influence**2 + distance**2 ) - numerator(i,j) = numerator(i,j) + weight * weight * obs(n) - denominator(i,j) = denominator(i,j) + weight - END IF - END DO - END DO - END DO obs_loop ! We have looped over all of the observations, and have computed the required @@ -187,46 +223,59 @@ ! This routine computes the distance from an observation to a grid point. ! The Euclidean distance is modified to allow for elongation due to ! strong stream-wise flow, as well as for curvature. +! The distance is then used to update the denominator and numerator used to +! determine the grid point perturbations for the objective analysis -SUBROUTINE dist ( i , j , iew , jns , crsdot , name , obs_value , xob , yob , dxd , & +SUBROUTINE dist_uu_vv_rh ( iew_min, iew_max, jns_min, jns_max, iew , jns , real_crsdot_divide_2, name , & + obs_n, xob , yob , dxd , & pressure , u_banana , v_banana , radius_influence , & - lat_center , choice , distance ) + lat_center , gridded, fg_value_n, use_first_guess, numerator_pert, denominator_pert, & !BPR BEGIN + scale_cressman_rh_decreases ) +! lat_center , gridded, fg_value_n, use_first_guess, numerator_pert, denominator_pert ) !BPR END IMPLICIT NONE ! Input variables. - INTEGER , INTENT(IN) :: i , & - j , & - iew , & - jns , & - crsdot , & - radius_influence - REAL , INTENT(IN) :: obs_value , & + INTEGER , INTENT(IN) :: iew_min , & + iew_max , & + jns_min , & + jns_max + INTEGER , INTENT(IN) :: iew , & + jns + REAL , INTENT(IN) :: real_crsdot_divide_2 + CHARACTER (LEN=8), INTENT(IN) :: name + REAL , INTENT(IN) :: obs_n , & xob , & yob , & dxd , & - pressure , & - lat_center - REAL , INTENT(IN) , DIMENSION(iew,jns) :: u_banana , & + pressure + REAL , INTENT(IN) , DIMENSION(:,:) :: u_banana , & v_banana - CHARACTER (LEN=8), INTENT(IN) :: name + INTEGER , INTENT(IN) :: radius_influence + REAL , INTENT(IN) :: lat_center + REAL, INTENT (IN), DIMENSION(:,:) :: gridded + REAL, INTENT (IN) :: fg_value_n + LOGICAL , INTENT (IN) :: use_first_guess + LOGICAL , INTENT (IN) :: scale_cressman_rh_decreases !BPR BEGIN/END + ! Input/Output variables. + REAL , INTENT(INOUT), DIMENSION (:,:) :: numerator_pert, & + denominator_pert - ! Output variables. - - REAL , INTENT(OUT) :: distance - CHARACTER (LEN=8) , INTENT(INOUT) :: choice - ! Local variables. + REAL :: real_i , real_j + REAL :: radius_influence_squared + REAL :: numerator_curvature, denominator_curvature REAL :: radius_of_curvature , k REAL :: u_ob , dudx , speed_ob , dydx REAL :: v_ob , dvdx - REAL :: numerator , denominator REAL :: x1 , x2 , y1 , y2 , d12 INTEGER :: i1 , i2 , j1 , j2 INTEGER :: ia, ja, ib, jb, ic, jc, id, jd REAL :: ua, va, ub, vb, uc, vc, ud, vd, vor + INTEGER :: ii,jj + REAL :: distance_squared INTEGER :: i_center , j_center @@ -239,104 +288,110 @@ REAL :: x , y , theta_diff - SAVE + LOGICAL :: USE_CIRCLE, USE_BANANA, USE_ELLIPSE + LOGICAL :: do_rh_scaling + + REAL :: one_divided_by_pi, one_divided_by_twopi + REAL :: weight + REAL :: scale_factor - IF ( choice .EQ. 'undefine' ) THEN + USE_CIRCLE = .FALSE. + USE_BANANA = .FALSE. + USE_ELLIPSE = .FALSE. + radius_influence_squared = REAL(radius_influence * radius_influence) - ! The critical speed to decide whether or not to do a simple - ! circular Cressman or one of the anisotropic schemes. If the - ! observation speed is < the critical speed, do a circular - ! Cressman distance computation. At 1000 hPa, the critical speed - ! is 5 m/s; this goes linearly to 15 m/s at 500 hPa. + ! The critical speed to decide whether or not to do a simple + ! circular Cressman or one of the anisotropic schemes. If the + ! observation speed is < the critical speed, do a circular + ! Cressman distance computation. At 1000 hPa, the critical speed + ! is 5 m/s; this goes linearly to 15 m/s at 500 hPa. + + IF ( pressure < 500. ) THEN + vcrit = 15. + ELSE + vcrit = 25. - 0.02 * pressure + END IF + + ! Compute the radius of curvature for the streamline passing through + ! this observation location. This is the information available from + ! eqn 4.16 in Manning and Haagenson, 1994 (MH94). Note that MH94 eqn 4.16 + ! is different from that computation given below. + + ! dy/dx = slope of streamline through grid point, which is v/u + ! d^2y/dx^2 = rate of change of slope of streamline through grid point, which + ! is d(v/u)/dx if taken along stream line trajectory, which is + ! u*dv/dx - v*du/dx + ! ----------------- + ! u^2 + ! k = curvature = abs(dy/dx) / ( 1 + (d^2y/dx^2)^2 ) ^(3/2) + ! radius of curvature = 1 / k + + ! Values of u and v nearest the observation point in the first-guess data and the + ! respective speed at this point. + + u_ob = u_banana(NINT(xob),NINT(yob)) + v_ob = v_banana(NINT(xob),NINT(yob)) + speed_ob = SQRT ( u_ob**2 + v_ob**2 ) + + ! We can see if we need to do any fancy computations by just looking at the speed + ! of the observation. If the first-guess speed is less than the cutoff, we do a + ! simple circular Cressman and vamoose. + + IF ( speed_ob .GE. vcrit ) THEN - IF ( pressure < 500. ) THEN - vcrit = 15. + ! The dy/dx has to be defined, so u can't be = zero. If it is, then we end up + ! doing an ellipse, not a banana, anyways. + + IF ( ABS(u_ob) .GE. 1.E-5 ) THEN + dydx = v_ob / u_ob ELSE - vcrit = 25. - 0.02 * pressure + dydx = 1000. END IF + + ! Differences are computed from approximately 3 grid points on either side of the + ! observation location, normal and tanget to the flow at the observation location. + + x1 = xob + 3. * ( u_ob/speed_ob ) + y1 = yob + 3. * ( v_ob/speed_ob ) + x2 = xob - 3. * ( u_ob/speed_ob ) + y2 = yob - 3. * ( v_ob/speed_ob ) + + ! These values need to be inside the domain since we will use them as indices + ! for the u and v first-guess winds. + + i1 = MAX ( MIN ( NINT(x1) , iew ) , 1 ) + j1 = MAX ( MIN ( NINT(y1) , jns ) , 1 ) + i2 = MAX ( MIN ( NINT(x2) , iew ) , 1 ) + j2 = MAX ( MIN ( NINT(y2) , jns ) , 1 ) + + ! And the distance between these two points (in grid point units) is + + d12 = SQRT( (REAL(i1) - REAL(i2))**2 + (REAL(j1) - REAL(j2))**2 ) + + ! So we can now compute the du/dx and the dv/dx, where we are assumming that + ! we are on a streamline through this observation point. We extend the + ! difference calculation in approximately 3 grid points to either side of + ! the observation point on the tangent to the streamline. + + dudx = ( u_banana(i1,j1) - u_banana(i2,j2) ) / ( d12 * dxd ) + dvdx = ( v_banana(i1,j1) - v_banana(i2,j2) ) / ( d12 * dxd ) - ! Compute the radius of curvature for the streamline passing through - ! this observation location. This is the information available from - ! eqn 4.16 in Manning and Haagenson, 1994 (MH94). Note that MH94 eqn 4.16 - ! is different from that computation given below. + ! Build the numerator and denominator of the curvature equation (k from above). + ! Since the radius of curvature is the inverse of the curvature, there is no + ! real need to do a temporary computation and then an inverse, so we just do it + ! directly (denominator / numerator). The radius of curvature is now in units + ! of (m). - ! dy/dx = slope of streamline through grid point, which is v/u - ! d^2y/dx^2 = rate of change of slope of streamline through grid point, which - ! is d(v/u)/dx if taken along stream line trajectory, which is - ! u*dv/dx - v*du/dx - ! ----------------- - ! u^2 - ! k = curvature = abs(dy/dx) / ( 1 + (d^2y/dx^2)^2 ) ^(3/2) - ! radius of curvature = 1 / k - - ! Values of u and v nearest the observation point in the first-guess data and the - ! respective speed at this point. - - u_ob = u_banana(NINT(xob),NINT(yob)) - v_ob = v_banana(NINT(xob),NINT(yob)) - speed_ob = SQRT ( u_ob**2 + v_ob**2 ) - - ! We can see if we need to do any fancy computations by just looking at the speed - ! of the observation. If the first-guess speed is less than the cutoff, we do a - ! simple circular Cressman and vamoose. - - IF ( speed_ob .GE. vcrit ) THEN - - ! The dy/dx has to be defined, so u can't be = zero. If it is, then we end up - ! doing an ellipse, not a banana, anyways. - - IF ( ABS(u_ob) .GE. 1.E-5 ) THEN - dydx = v_ob / u_ob - ELSE - dydx = 1000. - END IF - - ! Differences are computed from approximately 3 grid points on either side of the - ! observation location, normal and tanget to the flow at the observation location. - - x1 = xob + 3. * ( u_ob/speed_ob ) - y1 = yob + 3. * ( v_ob/speed_ob ) - x2 = xob - 3. * ( u_ob/speed_ob ) - y2 = yob - 3. * ( v_ob/speed_ob ) - - ! These values need to be inside the domain since we will use them as indices - ! for the u and v first-guess winds. - - i1 = MAX ( MIN ( NINT(x1) , iew ) , 1 ) - j1 = MAX ( MIN ( NINT(y1) , jns ) , 1 ) - i2 = MAX ( MIN ( NINT(x2) , iew ) , 1 ) - j2 = MAX ( MIN ( NINT(y2) , jns ) , 1 ) - - ! And the distance between these two points (in grid point units) is - - d12 = SQRT( (REAL(i1) - REAL(i2))**2 + (REAL(j1) - REAL(j2))**2 ) - - ! So we can now compute the du/dx and the dv/dx, where we are assumming that - ! we are on a streamline through this observation point. We extend the - ! difference calculation in approximately 3 grid points to either side of - ! the observation point on the tangent to the streamline. - - dudx = ( u_banana(i1,j1) - u_banana(i2,j2) ) / ( d12 * dxd ) - dvdx = ( v_banana(i1,j1) - v_banana(i2,j2) ) / ( d12 * dxd ) - - ! Build the numerator and denominator of the curvature equation (k from above). - ! Since the radius of curvature is the inverse of the curvature, there is no - ! real need to do a temporary computation and then an inverse, so we just do it - ! directly (denominator / numerator). The radius of curvature is now in units - ! of (m). - - numerator = ABS ( u_ob * dvdx - v_ob * dudx ) / u_ob**2 - denominator = ( SQRT ( 1. + dydx**2 ) ) ** 3 - - IF ( ( ABS ( u_ob) .LT. 1.E-5 ) .OR. ( ABS ( v_ob) .LT. 1.E-5 ) .OR. & - ( ABS ( numerator ) .LT. 1.E-5 ) ) THEN - radius_of_curvature = 1000. - ELSE - ! k = numerator / denominator - ! radius_of_curvature = 1. / k - radius_of_curvature = denominator / numerator - END IF - + numerator_curvature = ABS ( u_ob * dvdx - v_ob * dudx ) / u_ob**2 + denominator_curvature = ( SQRT ( 1. + dydx**2 ) ) ** 3 + + IF ( ( ABS ( u_ob) .LT. 1.E-5 ) .OR. ( ABS ( v_ob) .LT. 1.E-5 ) .OR. & + ( ABS ( numerator_curvature ) .LT. 1.E-5 ) ) THEN + radius_of_curvature = 1000. + ELSE +! k = numerator / denominator +! radius_of_curvature = 1. / k + radius_of_curvature = denominator_curvature / numerator_curvature END IF END IF @@ -347,10 +402,7 @@ IF ( speed_ob .LT. vcrit ) THEN - distance = SQRT ( ( ( xob - REAL(crsdot)/2. - REAL(i) ) * ( xob - REAL(crsdot)/2. - REAL(i) ) ) + & - ( ( yob - REAL(crsdot)/2. - REAL(j) ) * ( yob - REAL(crsdot)/2. - REAL(j) ) ) ) - - choice = 'circle ' + USE_CIRCLE = .true. ! If the wind speed is large enough and the radius of curvature is small enough, then ! the banana shaped ellipse is used. @@ -357,108 +409,230 @@ ELSE IF ( ABS(radius_of_curvature) .LT. ( 3.* radius_influence ) * dxd ) THEN - IF ( choice .EQ. 'undefine' ) THEN + USE_BANANA = .true. - ! The radius of curvature needs to have a sign associated with it so that - ! from the given point on the radius, we can find the center. From the - ! observation location, find the four surrounding points that are three - ! grid units away. However, as before, they need to be in the domain. + + one_divided_by_twopi = 1. / twopi + one_divided_by_pi = 1. / pi + + ! The radius of curvature needs to have a sign associated with it so that + ! from the given point on the radius, we can find the center. From the + ! observation location, find the four surrounding points that are three + ! grid units away. However, as before, they need to be in the domain. + + ia = MAX ( NINT ( xob ) - 3 , 1 ) + ja = yob + ib = xob + jb = MAX ( NINT ( yob ) - 3 , 1 ) + ic = MIN ( NINT ( xob ) + 3 , iew ) + jc = yob + id = xob + jd = MIN ( NINT ( yob ) + 3 , jns ) + + ! With these four grid points, get the first-guess u and v at these points. + + ua = u_banana(ia,ja) + va = v_banana(ia,ja) + ub = u_banana(ib,jb) + vb = v_banana(ib,jb) + uc = u_banana(ic,jc) + vc = v_banana(ic,jc) + ud = u_banana(id,jd) + vd = v_banana(id,jd) + + ! Get the SIGN of the relative vorticity at the observation point with + ! these surrounding values. - ia = MAX ( NINT ( xob ) - 3 , 1 ) - ja = yob - ib = xob - jb = MAX ( NINT ( yob ) - 3 , 1 ) - ic = MIN ( NINT ( xob ) + 3 , iew ) - jc = yob - id = xob - jd = MIN ( NINT ( yob ) + 3 , jns ) - - ! With these four grid points, get the first-guess u and v at these points. - - ua = u_banana(ia,ja) - va = v_banana(ia,ja) - ub = u_banana(ib,jb) - vb = v_banana(ib,jb) - uc = u_banana(ic,jc) - vc = v_banana(ic,jc) - ud = u_banana(id,jd) - vd = v_banana(id,jd) - - ! Get the SIGN of the relative vorticity at the observation point with - ! these surrounding values. + vor = ( vc - va ) - ( ud - ub ) + + ! Modify the SIGN of the radius_of_curvature with the SIGN of the relative + ! vorticity. + + radius_of_curvature = radius_of_curvature * SIGN ( 1.,vor ) + + ! Compute the location of the center of curvature for the streamline + ! passing through this observation location. This is the same as eqn 4.20a + ! and 4.20b in MH94, except that i0, ic, j0, jc are not in the reversed + ! model directions (i is east west, j is north south). + + i_center = NINT(xob) - radius_of_curvature * v_ob / speed_ob + j_center = NINT(yob) + radius_of_curvature * u_ob / speed_ob - vor = ( vc - va ) - ( ud - ub ) - - ! Modify the SIGN of the radius_of_curvature with the SIGN of the relative - ! vorticity. - - radius_of_curvature = radius_of_curvature * SIGN ( 1.,vor ) - - ! Compute the location of the center of curvature for the streamline - ! passing through this observation location. This is the same as eqn 4.20a - ! and 4.20b in MH94, except that i0, ic, j0, jc are not in the reversed - ! model directions (i is east west, j is north south). - - i_center = NINT(xob) - radius_of_curvature * v_ob / speed_ob - j_center = NINT(yob) + radius_of_curvature * u_ob / speed_ob - - END IF ! Compute the angle from the positve x-direction to the line connecting the - ! center of the curvature to a couple of different points - the observation and - ! the (i,j) locations in the window. The variable theta_ob is the same as the - ! output from eqn 4.21 MH94, and theta_ij is eqn 4.22 from MH94. + ! center of the curvature to the observation. + ! The variable theta_ob is the same as the output from eqn 4.21 MH94. theta_ob = ATAN2 ( REAL(j_center)-yob , REAL(i_center)-xob ) - theta_ij = ATAN2 ( REAL(j_center)-REAL(j) , REAL(i_center)-REAL(i) ) + elongation = 1. + 0.7778 / vcrit * speed_ob + + + ! We have the fast enough speed (the observation speed is larger than the critical + ! speed), but too large radius of curvature (radius_curvature > 3X radius of influence). + ! This is, therefore, the ellipse scheme. + + ELSE + + USE_ELLIPSE = .true. + + elongation = 1. + 0.7778 / vcrit * speed_ob + + END IF + + + !Determine if RH scaling should be done + IF ( ( scale_cressman_rh_decreases ) .AND. & + ( ( name(1:8) .EQ. 'RH ' ) .AND. ( use_first_guess ) .AND. ( obs_n .lt. 0.0) ) ) THEN + do_rh_scaling = .TRUE. + ELSE + do_rh_scaling = .FALSE. + ENDIF + + !Now that we have done everything we can do outside the loop, loop over all + !relevant x, y points for whatever shape was chosen + JLOOP: DO jj = jns_min , jns_max + real_j = REAL(jj) + ILOOP: DO ii = iew_min , iew_max + real_i = REAL(ii) + + IF( USE_CIRCLE ) THEN + + distance_squared = & + ( ( xob - real_crsdot_divide_2 - real_i ) * ( xob - real_crsdot_divide_2 - real_i ) ) + & + ( ( yob - real_crsdot_divide_2 - real_j ) * ( yob - real_crsdot_divide_2 - real_j ) ) + + ELSEIF( USE_BANANA ) THEN + + ! Compute the angle from the positve x-direction to the line connecting the + ! center of the curvature to the (i,j) locations in the window. + ! theta_ij is eqn 4.22 from MH94. + + theta_ij = ATAN2 ( REAL(j_center)-real_j , REAL(i_center)-real_i ) + ! Compute the distance from the center of the curvature to the grid point ! (i,j). This corresponds to the definition given for rij on the top of page 63 ! in MH94. - rij = SQRT ( (REAL(i_center)-REAL(i))**2 + (REAL(j_center)-REAL(j))**2 ) + rij = SQRT ( (REAL(i_center)-real_i)**2 + (REAL(j_center)-real_j)**2 ) - ! All of the required pieces are assembled for use in MH94 eqns 4.15a, + ! All of the required pieces are assembled for use in MH94 eqns 4.15a, ! 4.14b and 4.14. First, make sure that the values of theta_diff are bounded ! by +- 2 pi. Then, bound the values between +- pi. theta_diff = theta_ob - theta_ij - IF ( theta_diff .GT. twopi ) theta_diff = theta_diff - REAL ( INT ( theta_diff/twopi ) ) * twopi - IF ( theta_diff .LT. -twopi ) theta_diff = theta_diff + REAL ( INT ( theta_diff/twopi ) ) * twopi + IF ( theta_diff .GT. twopi ) theta_diff = theta_diff - & + REAL ( INT ( theta_diff*one_divided_by_twopi ) ) * twopi + IF ( theta_diff .LT. -twopi ) theta_diff = theta_diff + & + REAL ( INT ( theta_diff*one_divided_by_twopi ) ) * twopi IF ( theta_diff .GT. pi ) theta_diff = theta_diff - twopi IF ( theta_diff .LT. -pi ) theta_diff = theta_diff + twopi ! x = radius_of_curvature * theta_diff - x = radius_of_curvature * theta_diff / pi ! this dividing by pi is home grown + x = radius_of_curvature * theta_diff * one_divided_by_pi ! this dividing by pi is home grown y = ABS(radius_of_curvature) - rij - elongation = 1. + 0.7778 / vcrit * speed_ob - distance = SQRT ( x*x/elongation + y*y ) + distance_squared = x*x/elongation + y*y - choice = 'banana ' + ELSEIF( USE_ELLIPSE ) THEN - ! We have the fast enough speed (the observation speed is larger than the critical - ! speed), but too large radius of curvature (radius_curvature > 3X radius of influence). - ! This is, therefore, the ellipse scheme. + x = ( u_ob * ( real_i - xob ) + v_ob * ( real_j - yob ) ) / speed_ob + y = ( v_ob * ( real_i - xob ) - u_ob * ( real_j - yob ) ) / speed_ob - ELSE + distance_squared = x*x/elongation + y*y - x = ( u_ob * ( REAL(i) - xob ) + v_ob * ( REAL(j) - yob ) ) / speed_ob - y = ( v_ob * ( REAL(i) - xob ) - u_ob * ( REAL(j) - yob ) ) / speed_ob + ELSE + PRINT *,'ERROR: dist_uu_vv_rh: It appears that no shape was chosen for radius of influence' + STOP + ENDIF - elongation = 1. + 0.7778 / vcrit * speed_ob - distance = SQRT ( x*x/elongation + y*y ) - choice = 'ellipse ' + ! If the observation is within a grid point's radius of influence, + ! sum its effect. - END IF + IF ( distance_squared .LT. radius_influence_squared ) THEN + scale_factor = 1.0 + IF ( do_rh_scaling ) THEN + !do_rh_scaling is TRUE if ALL of the following are TRUE: + !1) Field is RH, 2) use_first_guess is true, 3) the observed RH was + !lower than the first-guess RH + !Calculate a scale_factor that is: + ! first-guess RH at current location (i.e., location where ob is being applied + ! DIVIDED BY + ! First-guess RH at the location of the observation + scale_factor = min( 1.0, gridded(ii,jj)/fg_value_n) + END IF + !weight = ( radius_influence**2 - distance**2 ) / ( radius_influence**2 + distance**2 ) + weight = ( radius_influence_squared - distance_squared ) / ( radius_influence_squared + distance_squared ) + numerator_pert(ii,jj) = numerator_pert(ii,jj) + weight * weight * obs_n * scale_factor + denominator_pert(ii,jj) = denominator_pert(ii,jj) + weight + END IF + + END DO ILOOP + END DO JLOOP -END SUBROUTINE dist +END SUBROUTINE dist_uu_vv_rh + +SUBROUTINE dist_not_uu_vv_rh ( iew_min, iew_max, jns_min, jns_max, real_crsdot_divide_2, name , & + obs_n, xob , yob , radius_influence , numerator, denominator ) + + IMPLICIT NONE + + ! Input variables. + + INTEGER , INTENT(IN) :: iew_min + INTEGER , INTENT(IN) :: iew_max + INTEGER , INTENT(IN) :: jns_min + INTEGER , INTENT(IN) :: jns_max + REAL , INTENT(IN) :: real_crsdot_divide_2 + CHARACTER (LEN=8), INTENT(IN) :: name + REAL , INTENT(IN) :: obs_n + REAL , INTENT(IN) :: xob + REAL , INTENT(IN) :: yob + INTEGER , INTENT(IN) :: radius_influence + ! Input/Output variables. + REAL , INTENT(INOUT), DIMENSION (:,:) :: denominator, & + numerator + + ! Local variables. + REAL :: real_i , real_j + REAL :: distance_squared, radius_influence_squared + INTEGER :: ii,jj + REAL :: weight + + radius_influence_squared = REAL(radius_influence * radius_influence) + + + JLOOP: DO jj = jns_min , jns_max + real_j = REAL(jj) + ILOOP: DO ii = iew_min , iew_max + real_i = REAL(ii) + + distance_squared = & + ( ( xob - real_crsdot_divide_2 - real_i ) * ( xob - real_crsdot_divide_2 - real_i ) ) + & + ( ( yob - real_crsdot_divide_2 - real_j ) * ( yob - real_crsdot_divide_2 - real_j ) ) + + IF ( distance_squared .LT. radius_influence_squared ) THEN + + weight = ( radius_influence_squared - distance_squared ) / ( radius_influence_squared + distance_squared ) + numerator(ii,jj) = numerator(ii,jj) + weight * weight * obs_n + denominator(ii,jj) = denominator(ii,jj) + weight + + END IF + + END DO ILOOP + END DO JLOOP + +END SUBROUTINE dist_not_uu_vv_rh + ! !--------------------------------------------------------------------------- -SUBROUTINE get_background_for_oa ( t , u , v , rh , slp_x , & +!BPR BEGIN +!SUBROUTINE get_background_for_oa ( t , u , v , rh , slp_x , & +SUBROUTINE get_background_for_oa ( t , u , v , rh , slp_x , pres , & +!BPR END dum2d , kp , name , & iew_alloc , jns_alloc , kbu_alloc ) @@ -487,6 +661,16 @@ CASE ( 'PMSL ' ) which_var_to_get CALL yx2xy ( slp_x , jns_alloc , iew_alloc , dum2d ) +!BPR BEGIN + CASE ( 'PSFC ' ) which_var_to_get + IF(kp.ne.1) THEN + PRINT *,'ERROR: get_background_for_oa: When get_background_for_oa is called' + PRINT *,' for PSFC the vertical level must be set to 1 but kp = ',kp + STOP + END IF + CALL yx2xy ( pres(1,1,kp) , jns_alloc , iew_alloc , dum2d ) +!BPR END + END SELECT which_var_to_get END SUBROUTINE get_background_for_oa @@ -568,6 +752,9 @@ CASE ( 'UU ' ) ; errm = 2. CASE ( 'VV ' ) ; errm = 2. CASE ( 'PMSL ' ) ; errm = 1. + !BPR BEGIN + CASE ( 'PSFC ' ) ; errm = 1. + !BPR END CASE ( 'TT ' ) ; errm = 0.5 CASE ( 'RH ' ) ; errm = 10. END SELECT @@ -636,7 +823,10 @@ ! !--------------------------------------------------------------------------- -SUBROUTINE put_background_from_oa ( t , u , v , rh , slp_x , & +!BPR BEGIN +!SUBROUTINE put_background_from_oa ( t , u , v , rh , slp_x , & +SUBROUTINE put_background_from_oa ( t , u , v , rh , slp_x , pres , & +!BPR END dum2d , kp , name , & iew_alloc , jns_alloc , kbu_alloc ) @@ -665,6 +855,16 @@ CASE ( 'PMSL ' ) which_var_to_get CALL yx2xy ( dum2d , iew_alloc , jns_alloc , slp_x ) + !BPR BEGIN + CASE ( 'PSFC ' ) which_var_to_get + IF(kp.ne.1) THEN + PRINT *,'ERROR: pot_background_from_oa: When put_background_from_oa is called' + PRINT *,' for PSFC the vertical level must be set to 1 but kp = ',kp + STOP + ENDIF + CALL yx2xy ( dum2d , iew_alloc , jns_alloc , pres(1,1,kp) ) + !BPR END + END SELECT which_var_to_get END SUBROUTINE put_background_from_oa @@ -806,5 +1006,146 @@ END SUBROUTINE smooth_5 !--------------------------------------------------------------------------- +!NOTE THAT THIS IS AN ELEMENTAL FUNCTION SO WE CAN PASS IT SCALARS OR ARRAYS +!AND THE RESULT WILL BE CONFORMABLE TO THE INPUT ARRAYS +!Find water vapor mixing ratio given relative humidity, temperature, and +!pressure +!Algorithm from subroutine mixing_ratio in final_analysis_module.F90 +!If error then returns -99999, since we cannot "STOP" because this is an +!elemental function +ELEMENTAL FUNCTION qv_from_rh(rh,tk,pressure) + IMPLICIT NONE + + REAL, INTENT(IN) :: rh !Relative humidity (%) + REAL, INTENT(IN) :: tk !Temperature (Kelvin) + REAL, INTENT(IN) :: pressure !Pressure (hPa) + REAL :: qv_from_rh !Water vapor mixing ratio (kg/kg) + + REAL :: es, qs + REAL :: rh_use + + REAL, PARAMETER :: EPS = 0.622 + REAL, PARAMETER :: SVP1 = 0.6112 + REAL, PARAMETER :: SVP2 = 17.67 + REAL, PARAMETER :: SVP3 = 29.65 + REAL, PARAMETER :: SVPT0 = 273.15 + + REAL, PARAMETER :: small_diff = 0.00001 + LOGICAL :: error_found + + error_found = .false. + + !Do some error checking on the inputs + + !RH + rh_use = rh + IF((rh.lt.0.0).or.(rh.gt.100.0)) THEN + !If its out of range by more than a little bit then error out + IF((rh.lt.(0.0-small_diff)).OR.(rh.gt.(100.0+small_diff))) THEN + !PRINT *,'ERROR: qv_from_rh: RH outside expected ranges: ',rh + error_found = .true. + qv_from_rh = -99999.0 + ELSE + !If its just slightly outside of range then set it to edge of range + rh_use=max(min(rh,100.0),0.0) + ENDIF + ENDIF + + !Temperature + IF(tk<50.0) THEN + !PRINT *,'ERROR: qv_from_rh: Temperature of less than 50 K passed in: ',tk + error_found = .true. + qv_from_rh = -99998.0 + ENDIF + + !Pressure + IF((pressure<0.0).or.(pressure>1200.0)) THEN + !PRINT *,'ERROR: qv_from_rh: Pressure outside expected range passed in: ',pressure,' (hPa aka mb)' + error_found = .true. + qv_from_rh = -99997.0 + ENDIF + + IF(error_found) THEN + !qv_from_rh = -99999.0 + ELSE + es = svp1 * 10. * EXP(svp2 * (tk - svpt0) / (tk - svp3)) + qs = eps * es / (pressure - es) + qv_from_rh = MAX(0.01 * rh_use * qs,0.0) + ENDIF + +END FUNCTION qv_from_rh + +!--------------------------------------------------------------------------- +!NOTE THAT THIS IS AN ELEMENTAL FUNCTION SO WE CAN PASS IT SCALARS OR ARRAYS +!AND THE RESULT WILL BE CONFORMABLE TO THE INPUT ARRAYS +!Find relative humidity given water vapor mixing ratio, temperature, and +!pressure +!Based on algorithm in qv_from_rh which is in turned based on algorithm from +!subroutine mixing_ratio in final_analysis_module.F90 +!If error then returns -99999, since we cannot "STOP" because this is an +!elemental function +ELEMENTAL FUNCTION rh_from_qv(qv,tk,pressure) + + IMPLICIT NONE + + REAL, INTENT(IN) :: qv !Water vapor mixing ratio (kg/kg) + REAL, INTENT(IN) :: tk !Temperature (Kelvin) + REAL, INTENT(IN) :: pressure !Pressure (hPa) + REAL :: rh_from_qv !Relative humidity (%) + + REAL :: es, qs + REAL :: qv_use + + REAL, PARAMETER :: EPS = 0.622 + REAL, PARAMETER :: SVP1 = 0.6112 + REAL, PARAMETER :: SVP2 = 17.67 + REAL, PARAMETER :: SVP3 = 29.65 + REAL, PARAMETER :: SVPT0 = 273.15 + + REAL, PARAMETER :: small_diff = 0.00001 + LOGICAL :: error_found + + error_found = .false. + + qv_use = qv + !Do some error checking on the inputs + + !Water vapor mixing ratio + IF(qv.lt.0.0) THEN + !If its out of range by more than a little bit then error out + IF(qv.lt.(0-small_diff)) THEN + !PRINT *,'ERROR: rh_from_qv: QV outside expected ranges: ',qv + error_found = .true. + ELSE + !If its just slightly outside of range then set it to edge of range + qv_use=0.0 + ENDIF + ENDIF + + !Temperature + IF(tk<50.0) THEN + !PRINT *,'ERROR: rh_from_qv: Temperature of less than 50 K passed in: ',tk + error_found = .true. + ENDIF + + !Pressure + IF((pressure<0.0).or.(pressure>1200.0)) THEN + !PRINT *,'ERROR: rh_from_qv: Pressure outside expected range passed in: ',pressure,' (hPa aka mb)' + error_found = .true. + ENDIF + + IF(error_found) THEN + rh_from_qv = -99999.0 + ELSE + es = svp1 * 10. * EXP(svp2 * (tk - svpt0) / (tk - svp3)) + qs = eps * es / (pressure - es) + rh_from_qv = MIN ( MAX( ( 100.0 * qv_use ) / qs,0.0), 100.0 ) + ENDIF + +END FUNCTION rh_from_qv + + + + END MODULE obj_analysis Index: src/ob_access_module.F90 =================================================================== --- src/ob_access_module.F90 (revision 374) +++ src/ob_access_module.F90 (working copy) @@ -1,14 +1,19 @@ MODULE ob_access USE observation +!BPR BEGIN +USE get_fg_at_point +!BPR END CONTAINS !------------------------------------------------------------------------------ - SUBROUTINE query_ob ( obs , date , time , & request_variable , request_level , request_qc_max , request_p_diff , & -value , qc ) +!BPR BEGIN +!value , qc ) +value , qc , fg_3d_t , fg_3d_h ) +!BPR END USE date_pack @@ -22,6 +27,11 @@ request_p_diff REAL :: value INTEGER :: qc +!BPR BEGIN + REAL , INTENT(IN) , OPTIONAL , DIMENSION(:,:,:) :: fg_3d_t, & + fg_3d_h + LOGICAL :: mark_as_false_extend +!BPR END TYPE ( measurement ) , POINTER :: next_one @@ -33,6 +43,20 @@ CHARACTER (LEN=19) :: obs_date , analysis_date , analysis_p1_date , analysis_m1_date LOGICAL :: close +!BPR BEGIN + REAL :: curpdiff, minpdiff + INTEGER :: curlev, uselev, maxlev + LOGICAL :: isin + +!BPR END +!BPR BEGIN + REAL :: rh_value, t_value, rh_qc, t_qc + REAL :: t_layer_average, t_value_at_sea_level + REAL :: psfc_value, psfc_qc + REAL :: elevation_value + INCLUDE 'constants.inc' +!BPR END + INCLUDE 'error.inc' INCLUDE 'missing.inc' INTERFACE @@ -41,6 +65,10 @@ not_missing ( r ) = ( ABS ( r - missing_r ) .GT. 1. ) +!BPR BEGIN + mark_as_false_extend = .false. +!BPR END + ! Compute the difference in times bewteen the requested ! data and time, and the observation's data and time. @@ -108,7 +136,7 @@ value = obs%ground%slp%data qc = obs%ground%slp%qc END IF - + CASE DEFAULT slp_vs_others ! All of the data that is not in the "ground" section is stored in @@ -168,7 +196,96 @@ value = next_one%meas%rh%data qc = next_one%meas%rh%qc END IF - + + !BPR BEGIN + !Added so that we can perform QC on dewpoint to improve QC of RH + CASE ( 'DEWPOINT' ) which_sfc_var + IF ( ( ( not_missing ( next_one%meas%rh%data ) ) .AND. & + ( next_one%meas%rh%qc .LT. request_qc_max ) ) .AND. & + ( ( not_missing ( next_one%meas%temperature%data ) ) .AND. & + ( next_one%meas%temperature%qc .LT. request_qc_max ) ) ) THEN + rh_value = next_one%meas%rh%data + rh_qc = next_one%meas%rh%qc + t_value = next_one%meas%temperature%data + t_qc = next_one%meas%temperature%qc + !Calculate dewpoint using formula in + !obs_sort_module.F90 + !Since log(0) is undefined, if RH is near zero then set it to + !some small value + IF( rh_value .LT. very_small_rh ) THEN + value = 1. / ( 1./t_value - Rv_over_L * LOG ( very_small_rh / 100. ) ) + ELSE + value = 1. / ( 1./t_value - Rv_over_L * LOG ( rh_value / 100. ) ) + END IF + + !Since purpose is to improve QC of RH, pull RH QC values + qc = rh_qc + END IF + + CASE ( 'PSFC ' ) which_sfc_var + !Although ob header stores surface pressure in obs%ground%psfc%data + !Obsgrid seems to use the pressure in the surface level of the ob + IF ( ( not_missing ( next_one%meas%pressure%data ) ) .AND. & + ( next_one%meas%pressure%qc .LT. request_qc_max ) ) THEN + value = next_one%meas%pressure%data + qc = next_one%meas%pressure%qc + END IF + + CASE ( 'PMSLPSFC' ) which_sfc_var + !This case means that we want to obtain the sea level pressure + !based on the observed surface pressure + !Although ob header stores surface pressure in obs%ground%psfc%data + !Obsgrid seems to use the pressure in the surface level of the ob + IF ( ( not_missing ( next_one%meas%pressure%data ) ) .AND. & + ( next_one%meas%pressure%qc .LT. request_qc_max ) .AND. & + ( .NOT. eps_equal ( obs%info%elevation , missing_r , 1. ) ) ) THEN + psfc_value = next_one%meas%pressure%data + psfc_qc = next_one%meas%pressure%qc + elevation_value = obs%info%elevation + IF ( ( not_missing ( next_one%meas%temperature%data ) ) .AND. & + ( next_one%meas%temperature%qc .LT. request_qc_max ) ) THEN + !If the ob has a surface temperature ob that is not bad + !then use that to define temperature at level of ob and + !standard atmosphere to get to sea level + t_value = next_one%meas%temperature%data + t_value_at_sea_level = t_value - dTdz_std_atm_lt_11000m * elevation_value + t_layer_average = 0.5 * (t_value_at_sea_level + t_value ) + ELSE + !If the ob does NOT have a good surface temperature ob + !Then use first guess field to determine temperature + IF( PRESENT(fg_3d_t) .AND. PRESENT(fg_3d_h) ) THEN + CALL get_fg_t_at_point(fg_3d_t,fg_3d_h,obs%location%latitude, & + obs%location%longitude, 0.5*elevation_value, t_layer_average) + !If average temperature of layer is missing and the height of the observation is + !below 11000m (the top of the level where the lapse rate is applicable) then use + !standard atmosphere to find temperature + IF ( ( ABS ( t_layer_average - missing_r ) .LT. 1. ) .AND. & + ( elevation_value .LE. 11000.0 ) ) THEN + t_value_at_sea_level = t_at_sea_level_std_atm + t_value = t_value_at_sea_level + dTdz_std_atm_lt_11000m * elevation_value + t_layer_average = 0.5 * (t_value_at_sea_level + t_value ) + END IF + + ELSE + !If somehow we do not have the first guess field then error out since + !this should not happen. + !standard atmosphere + PRINT *,'ERROR: query_ob: First guess temperature and/or height are not available' + STOP + !In theory we could use standard atmosphere to get temperature + t_value_at_sea_level = t_at_sea_level_std_atm + t_value = t_value_at_sea_level + dTdz_std_atm_lt_11000m * elevation_value + t_layer_average = 0.5 * (t_value_at_sea_level + t_value ) + END IF + END IF + IF ( not_missing ( t_layer_average ) ) THEN + !Value is the sea level pressure calculated from the + !surface pressure, elevation, and temperatures + value = psfc_value / (exp( ( -elevation_value * g ) / (gasr * t_layer_average ) ) ) + qc = psfc_qc + END IF + END IF + !BPR END END SELECT which_sfc_var ! We found the surface data, and we either have the data or not. @@ -206,12 +323,108 @@ next_one => obs%surface still_more = ASSOCIATED ( next_one ) - + + !BPR BEGIN + !Original code assumes that a given ob will be close enough to + !at most one of the pressure levels in the background field + !If one expands the tolerance (request_p_diff) enough then one + !must determine if current pressure level of analysis is the + !one closest to this ob + curlev=0 + uselev=-1 + minpdiff=9999999 + search_for_levela : DO WHILE ( still_more ) + curlev=curlev+1 + curpdiff=ABS ( next_one%meas%pressure%data - request_level * 100 ) + found_levela : IF ( curpdiff .LT. request_p_diff ) THEN + IF(curpdiff.LT.minpdiff) THEN + minpdiff=curpdiff + uselev=curlev + ENDIF + END IF found_levela + next_one => next_one%next + still_more = ASSOCIATED ( next_one ) + END DO search_for_levela + + next_one => obs%surface + still_more = ASSOCIATED ( next_one ) + maxlev=curlev + !BPR END + ! Does the level exist? Loop through the linked list of the data ! to find a correct pressure. + !BPR BEGIN + curlev=0 + !BPR END search_for_level : DO WHILE ( still_more ) - found_level : IF ( ABS ( next_one%meas%pressure%data - request_level * 100 ) .LT. request_p_diff ) THEN + !BPR BEGIN + curlev=curlev+1 + curpdiff=ABS ( next_one%meas%pressure%data - request_level * 100 ) + !found_level : IF ( ABS ( next_one%meas%pressure%data - request_level * 100 ) .LT. request_p_diff ) THEN + !found_level : IF ( curlev .EQ. uselev ) THEN + !If single-level ob (e.g., ACARS) and the pressure is + !within request_p_diff Pa of analysis level we are working on + !OR + !If multiple-level ob (e.g., sonde) and the pressure is + !very close to the analysis level (within 1 Pa) + !AND + !It is not a surface ob (excluded by ensuring that the height + !of the ob is not identical to the station elevation, + !unless the height of the ob is missing) + ! + !Single-level obs need larger tolerances in order to be + !included since they cannot be interpolated to analysis pressures + !Note that some single-level obs (FM-97 AIREP and FM-88 SATOB) will + !be extended to the closest analysis level if possible by + !extend_to_closest if use_p_tolerance_one_lev = .FALSE. + !If this has been done they will have + !the original level and the new level and so they will be + !treated with the multiple-level obs + + !Multiple-level obs cannot use the broader tolerance, because they + !have been pre-interpolated to the analysis levels. If we + !use the broader tolerance we may include both the + !non-interpolated levels close to the analysis level and + !the new interpolated value at the analysis level + + !Surface obs are excluded here since they should be + !qc'd/analyzed at the surface level (which is dealt with in + !the not_slp_which_level switch under the case where + !request_level=1001). Here we are under the DEFAULT case. + !In the original code, surface obs were basically excluded + !from this section of code because surface obs would only be + !included if their pressure was within 0Pa of analysis + !level. With the pressure tolerance expanded for + !single-level obs, surface obs would be much more likely to + !be included in another level if surface obs were not + !specifically excluded. + found_level : IF (((( maxlev .EQ. 1 ).AND.( curlev .EQ. uselev )).OR.& + (( maxlev .NE. 1 ).AND.( curpdiff .LT. 1 ))).AND.& + ((.NOT. eps_equal ( next_one%meas%height%data , obs%info%elevation , 1. ) ) .OR. & + ( eps_equal ( next_one%meas%height%data , missing_r , 1. ) ))) THEN + !If ob qualifies because it is a single-level non-surface + !ob close enough to desired level + IF (( maxlev .EQ. 1 ).AND.( curlev .EQ. uselev )) THEN + !Note that this ob should be marked with the qc flag + !usually used for single-level obs that are assigned + !the closest analysis pressure level + !This will allow ob to pass following remove_unverified test on + !output in obs_sort_module.F90: + !IF ( is_sounding .AND. next%meas%pressure%qc .gt. 4 ) keep_data = .TRUE. + mark_as_false_extend = .true. + ENDIF + IF ( mark_as_false_extend ) THEN + !Check if pressure qc flag already contains extend_influence flag + isin = contains_2n ( next_one%meas%pressure%qc, extend_influence ) + !If pressure qc flag does NOT already contain the + !extend_influence flag, then add it + IF ( .NOT. isin ) THEN + next_one%meas%pressure%qc = next_one%meas%pressure%qc + extend_influence + ENDIF + ENDIF + + !BPR END ! Does the data exist on this level? We found the right ! pressure level, now check the existence of the @@ -247,6 +460,32 @@ value = next_one%meas%rh%data qc = next_one%meas%rh%qc END IF + + !BPR BEGIN + !Added so that we can perform QC on dewpoint to improve QC of RH + CASE ( 'DEWPOINT' ) which_var + IF ( ( ( not_missing ( next_one%meas%rh%data ) ) .AND. & + ( next_one%meas%rh%qc .LT. request_qc_max ) ) .AND. & + ( ( not_missing ( next_one%meas%temperature%data ) ) .AND. & + ( next_one%meas%temperature%qc .LT. request_qc_max ) ) ) THEN + rh_value = next_one%meas%rh%data + rh_qc = next_one%meas%rh%qc + t_value = next_one%meas%temperature%data + t_qc = next_one%meas%temperature%qc + !Calculate dewpoint using formula in + !obs_sort_module.F90 + !Since log(0) is undefined, if RH is near zero then set it to + !some small value + IF( rh_value .LT. very_small_rh ) THEN + value = 1. / ( 1./t_value - Rv_over_L * LOG ( very_small_rh / 100. ) ) + ELSE + value = 1. / ( 1./t_value - Rv_over_L * LOG ( rh_value / 100. ) ) + END IF + !Since the purpose is to improve the QC of RH, + !pull RH QC values + qc = rh_qc + END IF + !BPR END END SELECT which_var @@ -284,15 +523,15 @@ error_message(1:31) = 'query_ob ' error_message(32:) = ' Wrong time for observation ' // & TRIM ( obs%location%id ) // ', at time = ' // obs%valid_time%date_char(1:12) // & - ' (ccyymmddhhmm).' + ' (ccyymmddhhmm).' fatal = .false. listing = .false. ! foo -! CALL error_handler ( error_number , error_message , & -! fatal , listing ) +! CALL error_handler ( error_number , error_message , & +! fatal , listing ) ! print*,obs%surface%meas%temperature%data, obs%surface%meas%temperature%qc END IF right_time - + END SUBROUTINE query_ob !------------------------------------------------------------------------------ @@ -316,6 +555,12 @@ CHARACTER ( LEN = 4 ) :: plevel_char + !BPR BEGIN + TYPE ( measurement ) , POINTER :: next_next_one + LOGICAL :: multi_level_ob + INTEGER :: request_p_diff_use + !BPR END + INCLUDE 'error.inc' INTERFACE INCLUDE 'error.int' @@ -379,6 +624,36 @@ CASE ( 'RH ' ) which_sfc_var next_one%meas%rh%data = value next_one%meas%rh%qc =qc + + !BPR BEGIN + CASE ( 'PMSLPSFC' ) which_sfc_var + !Do not change the value because this case means that the program + !requested a sea level pressure based on the surface pressure + !so that quality control of surface pressure could be done since + !directly QC'ing surface pressure is difficult due to it's elevation + !dependence. + !So although we want to update the surface pressure QC flag, we do not + !want to put a derived sea level pressure (which is what is in value) + !into surface pressure + !Surface pressure is in the header but is also the + !pressure of the surface level so we need to adjust both + obs%ground%psfc%qc = qc + next_one%meas%pressure%qc = qc + !BPR END + + !BPR BEGIN + !Dewpoint is the moisture variable read in from the + !littler file, but RH is used internally and output to + !the OBS_DOMAIN file. However, the qc_obs* files use + !dewpoint so dewpoint is also important + !When we do QC on dewpoint we are trying to check for + !bad RH. So if we are here we are trying to update + !the RH QC flags based on the dewpoint QC + CASE ( 'DEWPOINT' ) which_sfc_var + next_one%meas%dew_point%data = value + next_one%meas%dew_point%qc =qc + next_one%meas%rh%qc =qc + !BPR END END SELECT which_sfc_var @@ -432,10 +707,38 @@ next_one => obs%surface still_more = ASSOCIATED ( next_one ) + !BPR BEGIN + !Determine the pressure tolerance to use when storing observations + !If this is a multi-level observation then one should effectively + ! use no tolerance since the ob has already been interpolated to + ! the pressure of the first-guess field so it should exactly match + !If this is a single-level (above-surface) observation then we + ! need to use the same tolerance used in pulling the ob in the + ! first place. However, we do not need to worry about determining + ! whether the current pressure level is the one closest to the ob + ! since we would not have retrieved this ob in the first place if + ! the current request_level was not the closest pressure level + !If the first level exists... + IF ( still_more ) THEN + !Make a pointer to the next level + next_next_one => next_one%next + !If this pointer is valid then this ob has at least two levels + multi_level_ob = ASSOCIATED ( next_next_one) + IF(multi_level_ob) THEN + request_p_diff_use = 1 + ELSE + request_p_diff_use = request_p_diff + END IF + END IF + !BPR END + ! Loop through the linked list of the data to find a correct pressure. search_for_level : DO WHILE ( still_more ) - found_level : IF ( ABS ( next_one%meas%pressure%data - request_level * 100 ) .LT. request_p_diff ) THEN + !BPR BEGIN + !found_level : IF ( ABS ( next_one%meas%pressure%data - request_level * 100 ) .LT. request_p_diff ) THEN + found_level : IF ( ABS ( next_one%meas%pressure%data - request_level * 100 ) .LT. request_p_diff_use ) THEN + !BPR END ! We found the right pressure level, so just shove the data in. @@ -457,6 +760,20 @@ next_one%meas%rh%data = value next_one%meas%rh%qc = qc + !BPR BEGIN + !Dewpoint is the moisture variable read in from the + !littler file, but RH is used internally and output to + !the OBS_DOMAIN file. However, the qc_obs* files use + !dewpoint so dewpoint is also important + !When we do QC on dewpoint we are trying to check for + !bad RH. So if we are here we are trying to update + !the RH QC flags based on the dewpoint QC + CASE ( 'DEWPOINT' ) which_var + next_one%meas%dew_point%data = value + next_one%meas%dew_point%qc = qc + next_one%meas%rh%qc = qc + !BPR END + END SELECT which_var ! Inside this IF, we found the level. We have stored the data, @@ -500,5 +817,5 @@ END SUBROUTINE store_ob !------------------------------------------------------------------------------ - + END MODULE ob_access Index: src/obs_sort_module.F90 =================================================================== --- src/obs_sort_module.F90 (revision 374) +++ src/obs_sort_module.F90 (working copy) @@ -426,7 +426,10 @@ ! ! ------------------------------------------------------------------------- -LOGICAL FUNCTION time_eq ( a , b , date , time ) +!BPR BEGIN +!LOGICAL FUNCTION time_eq ( a , b , date , time ) +LOGICAL FUNCTION time_eq ( a , b , date , time, ob_closest_to_target_time ) +!BPR END ! This defines operator .EQ. for 'time_info' data type @@ -436,9 +439,24 @@ TYPE ( time_info ) , INTENT ( INOUT ) :: a , b INTEGER , INTENT ( IN ) :: date , time + !BPR BEGIN + INTEGER , INTENT ( OUT ) :: ob_closest_to_target_time + !BPR END ! Local variables. + !BPR BEGIN + !This defines how close two times need to be to be equal + !In the standard version of OBSGRID this is 30 minutes (1800 seconds), but + !this results in 15-minute surface obs being merged and this may not be + !desirable + !The assumption that the first of two obs to be merged is correct (all else + !being equal) could conceivably introduce a bias, especially since the time + !of the merged ob appears to be assigned the current "target time"). + !INTEGER, PARAMETER ::time_equal_tolerance_seconds = 300 + INTEGER, PARAMETER ::time_equal_tolerance_seconds = 1800 + !BPR END + CHARACTER (LEN=19) :: target_date , a_date , b_date INTEGER :: diff_seconds , a_diff_seconds , b_diff_seconds @@ -472,35 +490,80 @@ CALL geth_idts ( a_date , b_date , diff_seconds ) + + !BPR BEGIN + !BEGIN OLD COMMENT ! If the times (a and b) are within half an hour of each other, we say that they ! are the same time. + !END OLD COMMENT + ! If the times (a and b) are within time_equal_tolerance_seconds of each other, we say that they + ! are the same time. + !IF ( ABS ( diff_seconds ) .LT. 1800 ) THEN + IF ( ABS ( diff_seconds ) .LT. time_equal_tolerance_seconds ) THEN + !BPR END - IF ( ABS ( diff_seconds ) .LT. 1800 ) THEN - + !BPR BEGIN OLD COMMENT ! Now that we know a and b are the same time, the important question is now ! are they the time that we want? If they are the same (which means either ! a or b is within an hour of the target time), we set both of these times ! to the target time. + !BPR END OLD COMMENT + CALL geth_idts ( target_date , a_date , a_diff_seconds ) CALL geth_idts ( target_date , b_date , b_diff_seconds ) + !BPR BEGIN + !Determine which ob of "duplicates" is closer to the target time + IF( ABS ( a_diff_seconds ) .LT. ABS ( b_diff_seconds ) ) THEN + ob_closest_to_target_time = 1 + ELSE + ob_closest_to_target_time = 2 + ENDIF + + + !Check if either ob is within 1 hour of target time, if so then use + !the time of the ob closest to the target time + !Previously it just used the target time. + !Advantages of the new method include: + ! 1) It keeps the actual time so subsequent "duplicate" obs can compare + ! their time against the previously merged ob to see if the new + ! duplicate is closer to the target time + ! 2) An original ob time is kept. so that the merging of say a 1224 UTC + ! 1226 UTC ob does not result in a merged 1200 UTC ob, but rather a + ! 1224 UTC ob. + !Disadvantage of the new method: + ! 1) In cases where one merged ob is prior to the ob time and one merged + ! ob is after the ob time, the new ob time will not represent any type + ! of rough average of the obs times. For example, previously a 1145 + ! UTC and 1214 UTC ob would be merged to form a 1200 UTC ob, whereas + ! now they would be merged to form a 1214 UTC bill. + !BPR END + IF ( ( ABS ( a_diff_seconds ) .LT. 3600 ) .OR. & ( ABS ( b_diff_seconds ) .LT. 3600 ) ) THEN - a%date_char( 1: 4) = target_date( 1: 4) - a%date_char( 5: 6) = target_date( 6: 7) - a%date_char( 7: 8) = target_date( 9:10) - a%date_char( 9:10) = target_date(12:13) - a%date_char(11:12) = target_date(15:16) - a%date_char(13:14) = target_date(18:19) + !BPR BEGIN + IF ( ob_closest_to_target_time == 1 ) THEN + b%date_char(1:14) = a%date_char(1:14) + ELSEIF ( ob_closest_to_target_time == 2 ) THEN + a%date_char(1:14) = b%date_char(1:14) + ENDIF - b%date_char( 1: 4) = target_date( 1: 4) - b%date_char( 5: 6) = target_date( 6: 7) - b%date_char( 7: 8) = target_date( 9:10) - b%date_char( 9:10) = target_date(12:13) - b%date_char(11:12) = target_date(15:16) - b%date_char(13:14) = target_date(18:19) + !a%date_char( 1: 4) = target_date( 1: 4) + !a%date_char( 5: 6) = target_date( 6: 7) + !a%date_char( 7: 8) = target_date( 9:10) + !a%date_char( 9:10) = target_date(12:13) + !a%date_char(11:12) = target_date(15:16) + !a%date_char(13:14) = target_date(18:19) + + !b%date_char( 1: 4) = target_date( 1: 4) + !b%date_char( 5: 6) = target_date( 6: 7) + !b%date_char( 7: 8) = target_date( 9:10) + !b%date_char( 9:10) = target_date(12:13) + !b%date_char(11:12) = target_date(15:16) + !b%date_char(13:14) = target_date(18:19) + !BPR END END IF time_eq = .TRUE. @@ -508,6 +571,9 @@ ELSE time_eq = .FALSE. + !BPR BEGIN + ob_closest_to_target_time = 0 + !BPR END END IF @@ -540,6 +606,9 @@ INTEGER , INTENT ( OUT ) :: total_dups INTEGER , INTENT ( IN ) :: date , & time + !BPR BEGIN + INTEGER :: ob_closest_to_target_time + !BPR END INCLUDE 'error.inc' INTERFACE @@ -585,10 +654,20 @@ ! If time fields are not completely identical, go to next observation. ! Sort is by location ONLY, not by time; so next+1 may be identical ! even though next has different time. + !BPR BEGIN + ! time_eq does NOT check that the "time fields are not completely + ! identical" as the above comment states. Rather it checks if they are + ! close enough, which in the default version was 30 minutes. This + ! tolerance is now a parameter in time_eq + !BPR END ! foo ! IF ( .NOT. ( obs(first)%valid_time .EQ. obs(second)%valid_time ) ) THEN - IF ( .NOT. time_eq ( obs(first)%valid_time , obs(second)%valid_time , date , time ) ) THEN +!BPR BEGIN +! IF ( .NOT. time_eq ( obs(first)%valid_time , obs(second)%valid_time , date , time ) ) THEN + IF ( .NOT. time_eq ( obs(first)%valid_time , obs(second)%valid_time , & + date , time , ob_closest_to_target_time ) ) THEN +!BPR END error_number = 0332001 error_message(1:31) = 'check_duplicate_ob ' error_message(32:) = ' Found multiple times for ' & @@ -605,11 +684,16 @@ ! Observations are from same location and time, so merge them. - CALL merge_obs ( obs(first) , obs(second) ) +!BPR BEGIN +! CALL merge_obs ( obs(first) , obs(second) ) + CALL merge_obs ( obs(first) , obs(second) , ob_closest_to_target_time ) +!BPR END ! Mark second of pair as discarded; data is put in 'first'. Note that ! a duplicate has been found by incrementing the counter. + + obs(second)%info%discard = .true. obs(first)%info%num_dups = obs(first)%info%num_dups + 1 total_dups = total_dups + 1 @@ -873,7 +957,10 @@ ! !------------------------------------------------------------------------------ -SUBROUTINE height_to_pres ( height , pressure , iew , jns , kbu , & +!BPR BEGIN +!SUBROUTINE height_to_pres ( height , pressure , iew , jns , kbu , & +SUBROUTINE height_to_pres ( height , pressure , slp_x , temperature , iew , jns , kbu , & +!BPR END map_projection , latitude , longitude , value ) ! This routine is used to compute the pressure of an observation @@ -885,6 +972,9 @@ USE map_utils USE map_utils_helper +!BPR BEGIN + USE get_fg_at_point +!BPR END IMPLICIT NONE @@ -892,6 +982,10 @@ map_projection REAL , DIMENSION ( jns , iew , kbu ) :: height REAL , DIMENSION ( kbu ) :: pressure + !BPR BEGIN + REAL , INTENT ( IN ) , DIMENSION ( jns, iew ) :: slp_x + REAL , INTENT ( IN ) , DIMENSION ( jns , iew , kbu ) :: temperature + !BPR END REAL :: latitude , & longitude TYPE ( meas_data ) :: value @@ -908,6 +1002,22 @@ LOGICAL :: found + !BPR BEGIN + REAL :: h_start, h_diff + REAL :: p_start + REAL :: t_start, t_ob, t_layer_average + INTEGER :: request_qc_max + LOGICAL :: not_missing + LOGICAL :: is_missing + REAL :: r + + INCLUDE 'constants.inc' + + INCLUDE 'missing.inc' + not_missing ( r ) = ( ABS ( r - missing_r ) .GT. 1. ) + is_missing ( r ) = ( ABS ( r - missing_r ) .LE. 1. ) + !BPR END + ! We can zoom directly out of here if the height with which we are ! going to do the approximation is a flag (meaning it is missing). @@ -922,9 +1032,18 @@ ! exit out of here without modifying the input data set (value). inside_domain : IF ( ( x_location .GT. 1 ) .AND. & - ( x_location .LT. iew - 2 ) .AND. & +!BPR BEGIN +!Cannot see reason why this cannot be -1 instead of -2. This will allow obs +!near the edge of the domain to use the first guess fields in estimating surface +!pressure and thus not need to rely on less accurate standard atmosphere method +! ( x_location .LT. iew - 2 ) .AND. & + ( x_location .LT. iew - 1 ) .AND. & +!BPR END ( y_location .GT. 1 ) .AND. & - ( y_location .LT. jns - 2 ) ) THEN +!BPR BEGIN +! ( y_location .LT. jns - 2 ) ) THEN + ( y_location .LT. jns - 1 ) ) THEN +!BPR END i = NINT ( x_location ) j = NINT ( y_location ) @@ -950,18 +1069,98 @@ END IF END DO find_trapping + !BPR BEGIN + !Original comment ! If we found the trapping levels, great. If not, we drop through this ! block and out of the routine. This would happen if an observation was ! reported very high in the atmosphere, beyond the analysis levels. We ! compute the new pressure from linear interpolation. The QC flag for the ! pressure reflects that this data was interpolated. + !New comment + ! If we found the trapping levels, interpolate the pressure. + ! If not, we now consider three cases: + ! 1) Ob height is between sea level and height of lowest non-surface + ! pressure level. + ! -Interpolate between these two layers + ! -Pressure QC flag marked to indicate pressure from first guess + ! 2) Ob height is below sea level + ! -Use temperature and height of ob with standard atmosphere lapse + ! rate to obtain temperature at sea level + ! -Use sea level pressure/height and estimated temperature along + ! with ob height/temperature to calculate pressure at ob + ! -Pressure QC flag marked to indicate pressure from first guess + ! since we only use std atmosphere lapse rate for temperature + ! 3) Ob height is above top level + ! -Do nothing here and let standard atmosphere routine deal with this + ! case + !BPR END found_level : IF ( found ) THEN value%pressure%data = EXP ( ( LOG ( p_below ) * ( h_ob - h_above) + LOG ( p_above ) * ( h_below - h_ob ) ) / & ( h_below - h_above ) ) value%pressure%qc = p_from_h_first_guess + !BPR BEGIN + ELSE + + !DO NOT USE THE SURFACE PRESSURE / HEIGHT BECAUSE THIS IS SENSITIVE TO + !THE HORIZONTAL INTERPOLATION FROM THE SOURCE OF THE FIRST GUESS TO + !THE CURRENT GRID + + !If the ob is between sea level and the 1st layer above the surface + IF ( ( h_ob .GE. 0.0 ) .AND. ( h_ob .LT. height(j,i,2) ) ) THEN + h_below = 0.0 + h_above = height(j,i,2) + p_below = slp_x(j,i) * 100. + p_above = pressure(2) * 100. + !If distance to interpolate pressure over is at least 1 meter + !then do the interpolation + IF(ABS(h_below-h_above).GE.1.0) THEN + value%pressure%data = EXP ( ( LOG ( p_below ) * ( h_ob - h_above) + LOG ( p_above ) * ( h_below - h_ob ) ) / & + ( h_below - h_above ) ) + ELSE !If distance to interpolate pressure over is less than 1 meter + !then simply average the two pressures. The interpolation equation may fail to + !produce acceptable results if interpolating over too small a vertical separation + value%pressure%data = 0.5 * ( p_below + p_above ) + END IF + value%pressure%qc = p_from_h_first_guess + + !If the ob is below sea level (and below the second layer) + ELSE IF ( h_ob .LT. 0.0 ) THEN + h_start = 0.0 + p_start = slp_x(j,i) * 100. + h_diff = h_start - h_ob + + !Set the QC flag value below which the temperature must be QC'd at + !in order to use it here to calculate pressure + !Note this appears to be BEFORE we do the error_max or buddy check so + !we do not know if the temperature ob will pass these QC checks but we + !check the QC flag anyway in case it indicates we should not use the + !temperature ob + request_qc_max = MIN ( fails_error_max , fails_buddy_check ) + IF ( ( not_missing ( value%temperature%data ) ) .AND. & + !If temperature ob is available, use temperature at ob height + ( value%temperature%qc .LT. request_qc_max ) ) THEN + t_ob = value%temperature%data + ELSE + !If temperature ob is missing or QC indicates problems get + !temperature from first guess field + CALL GET_FG_T_AT_POINT(temperature,height,latitude,longitude,h_ob,t_ob) + !If this fails then use standard atmosphere + IF ( is_missing ( t_ob ) ) THEN + t_ob = t_at_sea_level_std_atm + dTdz_std_atm_lt_11000m * h_ob + END IF + END IF + !Find temperature at lowest level to estimate using standard atmosphere lapse rate the + !temperature at the ob + t_start = t_ob + dTdz_std_atm_lt_11000m * h_diff + t_layer_average = 0.5 * (t_start + t_ob ) + value%pressure%data = p_start / (exp( ( -h_diff * g ) / (gasr * t_layer_average ) ) ) + value%pressure%qc = p_from_h_first_guess + END IF + + !BPR END END IF found_level - + END IF inside_domain END IF missing_height @@ -1990,6 +2189,7 @@ ELSE IF ( field1%qc .EQ. field2%qc ) THEN + ! Second, if they have the same quality control values, use data ! that was chosen 'best' @@ -2300,7 +2500,10 @@ ! ! -------------------------------------------------------------------------- -SUBROUTINE merge_obs ( first , second ) +!BPR BEGIN +!SUBROUTINE merge_obs ( first , second ) +SUBROUTINE merge_obs ( first , second , ob_closest_to_target_time ) +!BPR END ! Reports 'first' and 'second' have been found to have same location and ! time, therefore they must be merged and one of them discarded. @@ -2312,6 +2515,9 @@ TYPE ( report ) , INTENT ( INOUT ) :: first , & second + !BPR BEGIN + INTEGER , INTENT ( IN ) :: ob_closest_to_target_time + !BPR END INTEGER :: best INCLUDE 'error.inc' @@ -2336,15 +2542,28 @@ ELSE IF ( first%info%seq_num .LT. second%info%seq_num ) THEN best = 2 ELSE - best = 1 - error_number =3321001 - error_message(1:31) = 'merge_obs ' - error_message(32:) = ' Arbitrarily assuming "first" obs & - &is better than "second" for ' // & - TRIM ( first%location%name ) // ' ' // & - TRIM ( first%location%id ) // '.' - fatal = .false. - listing = .false. +!BPR BEGIN + !Previously this first branch is all that existed. + !Now, if we have no reason to pick one ob over the other, choose the ob + !that is closest to the target date. I believe that the target date is the date/time + !whose analysis you are QC'ing the obs against + !Only revert to old code if ob_closest_to_target_time was set to zero (which I + !do not think should ever be true at this point in the code) + IF( (ob_closest_to_target_time .ne. 1) .and. (ob_closest_to_target_time .ne. 2) ) THEN + best = 1 + error_number =3321001 + error_message(1:31) = 'merge_obs ' + error_message(32:) = ' Arbitrarily assuming "first" obs & + &is better than "second" for ' // & + TRIM ( first%location%name ) // ' ' // & + TRIM ( first%location%id ) // '.' + fatal = .false. + listing = .false. + ELSE + best = ob_closest_to_target_time + ENDIF +!BPR END + ! CALL error_handler ( error_number , error_message , & ! fatal , listing ) END IF @@ -2500,9 +2719,12 @@ ! ! ------------------------------------------------------------------------- +!BPR BEGIN SUBROUTINE output_obs ( obs , unit , file_name , num_obs , out_opt, forinput, & new_file, qc_flag_keep, remove_unverified, num_klvls, pressure, & - met_file) + met_file, use_grid_relative_winds ) +! met_file ) +!BPR END ! Take the array of observations and write them including measurements ! at all levels. The two options (out_opt and forinput) are described @@ -2525,6 +2747,9 @@ INTEGER , INTENT ( IN ) :: unit CHARACTER ( LEN = * ) , INTENT ( IN ) :: file_name LOGICAL , INTENT ( IN ) :: forinput, new_file + !BPR BEGIN + LOGICAL , INTENT ( IN ) :: use_grid_relative_winds + !BPR END INTEGER :: i , iout, num_klvls INTEGER :: qc_flag_keep @@ -2536,7 +2761,7 @@ LOGICAL :: keep_data LOGICAL :: OBS_data=.FALSE. LOGICAL :: is_sounding - LOGICAL :: no_qc_done + LOGICAL :: no_qc_done INTEGER :: true_num_obs INTEGER :: track_surface_data @@ -2545,10 +2770,10 @@ INTEGER :: ilev CHARACTER (LEN=132) :: nfile REAL, ALLOCATABLE, DIMENSION(:) :: press_nc, z_nc, temp_nc, td_nc - REAL, ALLOCATABLE, DIMENSION(:) :: spd_nc, dir_nc, u_nc, v_nc, rh_nc + REAL, ALLOCATABLE, DIMENSION(:) :: spd_nc, dir_nc, u_nc, v_nc, rh_nc INTEGER, ALLOCATABLE, DIMENSION(:) :: press_qc_nc, z_qc_nc, temp_qc_nc, td_qc_nc - INTEGER, ALLOCATABLE, DIMENSION(:) :: spd_qc_nc, dir_qc_nc, u_qc_nc, v_qc_nc, rh_qc_nc + INTEGER, ALLOCATABLE, DIMENSION(:) :: spd_qc_nc, dir_qc_nc, u_qc_nc, v_qc_nc, rh_qc_nc INTEGER :: int_sounding, int_bogus, int_discard, iout_nc CHARACTER ( LEN = 132 ) :: met_file @@ -2555,8 +2780,14 @@ INTEGER :: met_ncid INTEGER :: idummy REAL :: rdummy - + !BPR BEGIN + TYPE ( meas_data ) :: meas_earth_winds, meas_write + !BPR END + !BPR BEGIN + LOGICAL :: qc_done_against_nearby_level + !BPR END + end_meas%pressure%data = end_data_r end_meas%height%data = end_data_r end_meas%temperature%data = end_data_r @@ -2837,7 +3068,14 @@ ENDIF !!!!IF ( obs(i)%info%num_vld_fld == 1 .AND. (obs(i)%info%elevation .ne. missing) .AND. & IF ( obs(i)%info%num_vld_fld == 1 .AND. & - ( next%meas%height%data .eq. obs(i)%info%elevation ) ) THEN + !BPR BEGIN + !Without inclusion of check to see if elevation is missing, + !single-level above-surface observations with the height and + !elevation set to missing will be converted to being not a sounding + !( next%meas%height%data .eq. obs(i)%info%elevation ) ) THEN + ( next%meas%height%data .eq. obs(i)%info%elevation ) .AND. & + (obs(i)%info%elevation .ne. missing) ) THEN + !BPR END is_sounding = .FALSE. obs(i)%info%is_sound = .FALSE. ENDIF @@ -2857,9 +3095,25 @@ ENDIF IF ( remove_unverified ) THEN keep_data = any ( ( pres_hPA .eq. next%meas%pressure%data ) ) - IF ( (obs(i)%info%elevation .ne. missing) .AND. & + IF ( (obs(i)%info%elevation .ne. missing) .AND. & (next%meas%height%data .eq. obs(i)%info%elevation) ) keep_data = .TRUE. !!IF ( is_sounding .AND. next%meas%pressure%qc .gt. 4 ) keep_data = .TRUE. + !BPR BEGIN + !keep_data is initially defined above by being set to true for + !any observations whose pressure matched one of the pressure + !levels in the first-guess field we are using to QC against + !Then it is set to true for all surface observations + !Now, check to see if the QC flag is set to the value (extend_influence) that + !indicates that either: + !1) The ob was assigned a nearby pressure where there is a level in the first guess + !2) The ob was QC'd against a nearby pressure where there is a level in the first guess + !If the first is true then the pressure would have been changed + !in the data so keep_data should already be true, but if the + !second is true then keep_data is currently set to false but + !should be set to true + qc_done_against_nearby_level = contains_2n ( next%meas%pressure%qc, extend_influence ) + IF ( is_sounding .AND. qc_done_against_nearby_level ) keep_data = .TRUE. + !BPR END IF ( keep_data ) true_num_obs = true_num_obs + 1 ELSE true_num_obs = true_num_obs + 1 @@ -2872,7 +3126,7 @@ !! 2012-12-17 cB - If surface and we don't want data above a set qc value - discard !! Also discard soundings with no data - if ( .NOT. is_sounding .and. (obs(i)%surface%meas%pressure%qc >= qc_flag_keep) ) & + if ( .NOT. is_sounding .and. (obs(i)%surface%meas%pressure%qc >= qc_flag_keep) ) & obs(i)%info%discard = .TRUE. if ( is_sounding .and. true_num_obs.eq.0) obs(i)%info%discard = .TRUE. @@ -2885,97 +3139,130 @@ endif IF ( OBS_data ) THEN - IF ( .NOT. obs(i)%info%discard ) THEN - IF ( is_sounding ) THEN - WRITE ( UNIT = unit , FMT='(1x,A14)' ) obs(i)%valid_time%date_char(1:14) - WRITE ( UNIT = unit , FMT='(2x,2(F9.4,1x))' ) obs(i)%location%latitude, obs(i)%location%longitude - IF ( obs(i)%location%id(1:5) == "US un" ) THEN - WRITE ( UNIT = unit , FMT='(2x,2(A40,3x))' ) & - "----- " , & - "Unknown Station " - ELSE - WRITE ( UNIT = unit , FMT='(2x,2(A40,3x))' ) obs(i)%location%id, obs(i)%location%name - ENDIF - WRITE ( UNIT = unit , FMT='(2x,2(A16,2x),F8.0,2x,2(L4,2x),I5)' ) & - obs(i)%info%platform, obs(i)%info%source, obs(i)%info%elevation, & - is_sounding, obs(i)%info%bogus, true_num_obs + IF ( .NOT. obs(i)%info%discard ) THEN + WRITE ( UNIT = unit , FMT='(1x,A14)' ) obs(i)%valid_time%date_char(1:14) + WRITE ( UNIT = unit , FMT='(2x,2(F9.4,1x))' ) obs(i)%location%latitude, obs(i)%location%longitude + IF ( obs(i)%location%id(1:5) == "US un" ) THEN + WRITE ( UNIT = unit , FMT='(2x,2(A40,3x))' ) & + "----- " , & + "Unknown Station " + ELSE + WRITE ( UNIT = unit , FMT='(2x,2(A40,3x))' ) obs(i)%location%id, obs(i)%location%name ENDIF - ENDIF + WRITE ( UNIT = unit , FMT='(2x,2(A16,2x),F8.0,2x,2(L4,2x),I5)' ) & + obs(i)%info%platform, obs(i)%info%source, obs(i)%info%elevation, & + is_sounding, obs(i)%info%bogus, true_num_obs + ENDIF ELSE - IF ( .NOT. obs(i)%info%discard ) THEN - WRITE ( UNIT = unit , FMT = rpt_format ) & - obs(i)%location , obs(i)%info , obs(i)%valid_time , obs(i)%ground + !BPR BEGIN + !Only write the header if this ob has not been marked as discard=TRUE + !or the output file we are writing is not supposed to throw out obs + !IF ( .NOT. obs(i)%info%discard ) THEN + IF ( out_opt .EQ. 0 .OR. & + ( out_opt .GT. 0 .AND. .NOT. obs(i)%info%discard ) .OR. & + ( out_opt .LT. 0 .AND. obs(i)%info%discard ) ) THEN + !BPR END + WRITE ( UNIT = unit , FMT = rpt_format ) & + obs(i)%location , obs(i)%info , obs(i)%valid_time , obs(i)%ground ENDIF ENDIF - + track_surface_data = 0 next => obs(i)%surface - IF ( .not. OBS_data ) THEN - ilev = 0 - if(allocated(press_nc)) deallocate(press_nc) - allocate(press_nc(true_num_obs)) - if(allocated(press_qc_nc)) deallocate(press_qc_nc) - allocate(press_qc_nc(true_num_obs)) - if(allocated(z_nc)) deallocate(z_nc) - allocate(z_nc(true_num_obs)) - if(allocated(z_qc_nc)) deallocate(z_qc_nc) - allocate(z_qc_nc(true_num_obs)) - if(allocated(temp_nc)) deallocate(temp_nc) - allocate(temp_nc(true_num_obs)) - if(allocated(temp_qc_nc)) deallocate(temp_qc_nc) - allocate(temp_qc_nc(true_num_obs)) - if(allocated(td_nc)) deallocate(td_nc) - allocate(td_nc(true_num_obs)) - if(allocated(td_qc_nc)) deallocate(td_qc_nc) - allocate(td_qc_nc(true_num_obs)) - if(allocated(spd_nc)) deallocate(spd_nc) - allocate(spd_nc(true_num_obs)) - if(allocated(spd_qc_nc)) deallocate(spd_qc_nc) - allocate(spd_qc_nc(true_num_obs)) - if(allocated(dir_nc)) deallocate(dir_nc) - allocate(dir_nc(true_num_obs)) - if(allocated(dir_qc_nc)) deallocate(dir_qc_nc) - allocate(dir_qc_nc(true_num_obs)) - if(allocated(u_nc)) deallocate(u_nc) - allocate(u_nc(true_num_obs)) - if(allocated(u_qc_nc)) deallocate(u_qc_nc) - allocate(u_qc_nc(true_num_obs)) - if(allocated(v_nc)) deallocate(v_nc) - allocate(v_nc(true_num_obs)) - if(allocated(v_qc_nc)) deallocate(v_qc_nc) - allocate(v_qc_nc(true_num_obs)) - if(allocated(rh_nc)) deallocate(rh_nc) - allocate(rh_nc(true_num_obs)) - if(allocated(rh_qc_nc)) deallocate(rh_qc_nc) - allocate(rh_qc_nc(true_num_obs)) + IF ( .not. OBS_data ) THEN + ilev = 0 + if(allocated(press_nc)) deallocate(press_nc) + allocate(press_nc(true_num_obs)) + if(allocated(press_qc_nc)) deallocate(press_qc_nc) + allocate(press_qc_nc(true_num_obs)) + if(allocated(z_nc)) deallocate(z_nc) + allocate(z_nc(true_num_obs)) + if(allocated(z_qc_nc)) deallocate(z_qc_nc) + allocate(z_qc_nc(true_num_obs)) + if(allocated(temp_nc)) deallocate(temp_nc) + allocate(temp_nc(true_num_obs)) + if(allocated(temp_qc_nc)) deallocate(temp_qc_nc) + allocate(temp_qc_nc(true_num_obs)) + if(allocated(td_nc)) deallocate(td_nc) + allocate(td_nc(true_num_obs)) + if(allocated(td_qc_nc)) deallocate(td_qc_nc) + allocate(td_qc_nc(true_num_obs)) + if(allocated(spd_nc)) deallocate(spd_nc) + allocate(spd_nc(true_num_obs)) + if(allocated(spd_qc_nc)) deallocate(spd_qc_nc) + allocate(spd_qc_nc(true_num_obs)) + if(allocated(dir_nc)) deallocate(dir_nc) + allocate(dir_nc(true_num_obs)) + if(allocated(dir_qc_nc)) deallocate(dir_qc_nc) + allocate(dir_qc_nc(true_num_obs)) + if(allocated(u_nc)) deallocate(u_nc) + allocate(u_nc(true_num_obs)) + if(allocated(u_qc_nc)) deallocate(u_qc_nc) + allocate(u_qc_nc(true_num_obs)) + if(allocated(v_nc)) deallocate(v_nc) + allocate(v_nc(true_num_obs)) + if(allocated(v_qc_nc)) deallocate(v_qc_nc) + allocate(v_qc_nc(true_num_obs)) + if(allocated(rh_nc)) deallocate(rh_nc) + allocate(rh_nc(true_num_obs)) + if(allocated(rh_qc_nc)) deallocate(rh_qc_nc) + allocate(rh_qc_nc(true_num_obs)) + + press_nc = missing_r + press_qc_nc = missing + z_nc = missing_r + z_qc_nc = missing + temp_nc = missing_r + temp_qc_nc = missing + td_nc = missing_r + td_qc_nc = missing + spd_nc = missing_r + spd_qc_nc = missing + dir_nc = missing_r + dir_qc_nc = missing + u_nc = missing_r + u_qc_nc = missing + v_nc = missing_r + v_qc_nc = missing + rh_nc = missing_r + rh_qc_nc = missing + ENDIF - press_nc = missing_r - press_qc_nc = missing - z_nc = missing_r - z_qc_nc = missing - temp_nc = missing_r - temp_qc_nc = missing - td_nc = missing_r - td_qc_nc = missing - spd_nc = missing_r - spd_qc_nc = missing - dir_nc = missing_r - dir_qc_nc = missing - u_nc = missing_r - u_qc_nc = missing - v_nc = missing_r - v_qc_nc = missing - rh_nc = missing_r - rh_qc_nc = missing - ENDIF - DO WHILE ( ASSOCIATED ( next ) ) - if ( obs(i)%info%discard ) exit + if ( obs(i)%info%discard ) exit keep_data = any ( ( pres_hPA .eq. next%meas%pressure%data ) ) IF ( (obs(i)%info%elevation .ne. missing) .AND. (next%meas%height%data .eq. obs(i)%info%elevation) ) keep_data = .TRUE. !!IF ( is_sounding .AND. next%meas%pressure%qc .gt. 4 ) keep_data = .TRUE. + !BPR BEGIN + !keep_data is initially defined above by being set to true for + !any observations whose pressure matched one of the pressure + !levels in the first-guess field we are using to QC against + !Then it is attempted to be set to true for all surface observations + ! However if a surface observation has a slightly different height + ! than the elevation of the station keep_data will NOT be set to + ! true. Sometimes these may differ by ~1m. HOWEVER, keep_data is + ! no longer used in determining whether to output surface data so + ! this does not matter + ! So that someone in the future does not rely on keep_data for + ! surface data and have it behave unexpectedly one could use the + ! following line to put a tolerance on the difference between station elevation and + ! height of ob BUT this might break cases where one has a sounding + ! with a value very close to the surface so will not use this line + ! at this time. + ! IF ( (obs(i)%info%elevation .ne. missing) .AND. & + ! eps_equal( next%meas%height%data, obs(i)%info%elevation, 3. ) ) keep_data = .TRUE. + !Now, check to see if the QC flag is set to the value (extend_influence) that + !indicates that either: + !1) The ob was assigned a nearby pressure where there is a level in the first guess + !2) The ob was QC'd against a nearby pressure where there is a level in the first guess + !If the first is true then the pressure would have been changed + !in the data so keep_data should already be true, but if the + !second is true then keep_data is currently set to false but + qc_done_against_nearby_level = contains_2n ( next%meas%pressure%qc, extend_influence ) + IF ( is_sounding .AND. qc_done_against_nearby_level ) keep_data = .TRUE. + !BPR END !!! Make sure no surface ob has more than one entry track_surface_data = track_surface_data + 1 @@ -3062,185 +3349,210 @@ IF (obs(i)%ground%precip%data == missing_r) obs(i)%ground%precip%qc = missing IF (next%meas%thickness%data == missing_r) next%meas%thickness%qc = missing - IF ( is_sounding ) THEN - IF ( (keep_data .AND. remove_unverified) .OR. (.NOT. remove_unverified) ) THEN - IF ( OBS_data ) THEN - WRITE ( UNIT = unit , FMT='(1x,6(F11.3,1x,F11.3,1x))' ) & - next%meas%pressure%data, real(next%meas%pressure%qc), & - next%meas%height%data, real(next%meas%height%qc), & - next%meas%temperature%data, real(next%meas%temperature%qc), & - next%meas%u%data, real(next%meas%u%qc), & - next%meas%v%data, real(next%meas%v%qc), & - next%meas%rh%data, real(next%meas%rh%qc) - ELSE - WRITE ( UNIT = unit , FMT = meas_format ) next%meas - ENDIF - IF ( file_name(1:10) == "qc_obs_use" ) THEN - ilev = ilev + 1 - press_nc(ilev) = next%meas%pressure%data - press_qc_nc(ilev) = next%meas%pressure%qc - z_nc(ilev) = next%meas%height%data - z_qc_nc(ilev) = next%meas%height%qc - temp_nc(ilev) = next%meas%temperature%data - temp_qc_nc(ilev) = next%meas%temperature%qc - td_nc(ilev) = next%meas%dew_point%data - td_qc_nc(ilev) = next%meas%dew_point%qc - spd_nc(ilev) = next%meas%speed%data - spd_qc_nc(ilev) = next%meas%speed%qc - dir_nc(ilev) = next%meas%direction%data - dir_qc_nc(ilev) = next%meas%direction%qc - u_nc(ilev) = next%meas%u%data - u_qc_nc(ilev) = next%meas%u%qc - v_nc(ilev) = next%meas%v%data - v_qc_nc(ilev) = next%meas%v%qc - rh_nc(ilev) = next%meas%rh%data - rh_qc_nc(ilev) = next%meas%rh%qc - ENDIF - ENDIF + ! BPR BEGIN + ! Rotate winds to earth-relative (keep non wind fields the same) and place + ! in meas_earth_winds + ! Also fills in speed/direction in grid-relative winds (next%meas) since the + ! direction is not necessarily correct at this point + CALL grid_to_earth_winds ( next%meas, meas_earth_winds, obs(i)%location%longitude ) + IF( use_grid_relative_winds ) THEN + meas_write = next%meas ELSE - IF ( (keep_data .AND. remove_unverified) .OR. (.NOT. remove_unverified) ) THEN - IF ( OBS_data ) THEN - WRITE ( UNIT = unit , FMT='(1x,A14)' ) obs(i)%valid_time%date_char(1:14) - WRITE ( UNIT = unit , FMT='(2x,2(F9.4,1x))' ) obs(i)%location%latitude, obs(i)%location%longitude - IF ( obs(i)%location%id(1:5) == "US un" ) THEN - WRITE ( UNIT = unit , FMT='(2x,2(A40,3x))' ) & - "----- " , & - "Unknown Station " - ELSE - WRITE ( UNIT = unit , FMT='(2x,2(A40,3x))' ) obs(i)%location%id, obs(i)%location%name - ENDIF - WRITE ( UNIT = unit , FMT='(2x,2(A16,2x),F8.0,2x,2(L4,2x),I5)' ) & - obs(i)%info%platform, obs(i)%info%source, obs(i)%info%elevation, & - is_sounding, obs(i)%info%bogus, true_num_obs - WRITE ( UNIT = unit , FMT='(1x,9(F11.3,1x,F11.3,1x))' ) & - obs(i)%ground%slp%data, real(obs(i)%ground%slp%qc), & - obs(i)%ground%ref_pres%data, real(obs(i)%ground%ref_pres%qc),& - next%meas%height%data, real(next%meas%height%qc), & - next%meas%temperature%data, real(next%meas%temperature%qc), & - next%meas%u%data, real(next%meas%u%qc), & - next%meas%v%data, real(next%meas%v%qc), & - next%meas%rh%data, real(next%meas%rh%qc), & - next%meas%pressure%data, real(next%meas%pressure%qc), & - obs(i)%ground%precip%data, real(obs(i)%ground%precip%qc) - ELSE - WRITE ( UNIT = unit , FMT = meas_format ) next%meas - ENDIF - IF ( file_name(1:10) == "qc_obs_use" ) THEN - ilev = ilev + 1 - press_nc(ilev) = next%meas%pressure%data - press_qc_nc(ilev) = next%meas%pressure%qc - z_nc(ilev) = next%meas%height%data - z_qc_nc(ilev) = next%meas%height%qc - temp_nc(ilev) = next%meas%temperature%data - temp_qc_nc(ilev) = next%meas%temperature%qc - td_nc(ilev) = next%meas%dew_point%data - td_qc_nc(ilev) = next%meas%dew_point%qc - spd_nc(ilev) = next%meas%speed%data - spd_qc_nc(ilev) = next%meas%speed%qc - dir_nc(ilev) = next%meas%direction%data - dir_qc_nc(ilev) = next%meas%direction%qc - u_nc(ilev) = next%meas%u%data - u_qc_nc(ilev) = next%meas%u%qc - v_nc(ilev) = next%meas%v%data - v_qc_nc(ilev) = next%meas%v%qc - rh_nc(ilev) = next%meas%rh%data - rh_qc_nc(ilev) = next%meas%rh%qc - ENDIF - ENDIF + meas_write = meas_earth_winds ENDIF - + ! BPR END - IF ( file_name(1:10) == "qc_obs_raw" ) THEN - ilev = ilev + 1 - press_nc(ilev) = next%meas%pressure%data - press_qc_nc(ilev) = next%meas%pressure%qc - z_nc(ilev) = next%meas%height%data - z_qc_nc(ilev) = next%meas%height%qc - temp_nc(ilev) = next%meas%temperature%data - temp_qc_nc(ilev) = next%meas%temperature%qc - td_nc(ilev) = next%meas%dew_point%data - td_qc_nc(ilev) = next%meas%dew_point%qc - spd_nc(ilev) = next%meas%speed%data - spd_qc_nc(ilev) = next%meas%speed%qc - dir_nc(ilev) = next%meas%direction%data - dir_qc_nc(ilev) = next%meas%direction%qc - u_nc(ilev) = next%meas%u%data - u_qc_nc(ilev) = next%meas%u%qc - v_nc(ilev) = next%meas%v%data - v_qc_nc(ilev) = next%meas%v%qc - rh_nc(ilev) = next%meas%rh%data - rh_qc_nc(ilev) = next%meas%rh%qc - ENDIF - next => next%next - END DO + !BPR BEGIN REORDERING + !The following code is reordered from the standard version to be less repetitive + !and hopefully clearer to read + !All of the old code is not included in commented out fashion because it would be confusing + !The reordering involved... + !-The keep_data/remove_unverified conditional is placed on the outside of the is_sounding + ! conditional since it is applied to both + !-The is_sounding conditional is placed within the OBS_data conditional since it only matters + ! for the obs nudging files + !-The qc_obs_raw* netcdf file writing is included with the qc_obs_used* netcdf file + ! writing because the code to write them is identical. While the keep_data/remove_unverified + ! conditional was previously not applied to qc_obs_raw it can be applied for simplicity since + ! remove_unverified will be FALSE for the qc_obs_raw file + IF ( (keep_data .AND. remove_unverified) .OR. (.NOT. remove_unverified) ) THEN + IF ( OBS_data ) THEN + IF ( is_sounding ) THEN + WRITE ( UNIT = unit , FMT='(1x,6(F11.3,1x,F11.3,1x))' ) & + !BPR BEGIN + !next%meas%pressure%data, real(next%meas%pressure%qc), & + !next%meas%height%data, real(next%meas%height%qc), & + !next%meas%temperature%data, real(next%meas%temperature%qc), & + !next%meas%u%data, real(next%meas%u%qc), & + !next%meas%v%data, real(next%meas%v%qc), & + !next%meas%rh%data, real(next%meas%rh%qc) + meas_write%pressure%data, real(meas_write%pressure%qc), & + meas_write%height%data, real(meas_write%height%qc), & + meas_write%temperature%data, real(meas_write%temperature%qc), & + meas_write%u%data, real(meas_write%u%qc), & + meas_write%v%data, real(meas_write%v%qc), & + meas_write%rh%data, real(meas_write%rh%qc) + !BPR END + ELSE !If not marked as a sounding + WRITE ( UNIT = unit , FMT='(1x,9(F11.3,1x,F11.3,1x))' ) & + !BPR BEGIN + !obs(i)%ground%slp%data, real(obs(i)%ground%slp%qc), & + !obs(i)%ground%ref_pres%data, real(obs(i)%ground%ref_pres%qc),& + !next%meas%height%data, real(next%meas%height%qc), & + !next%meas%temperature%data, real(next%meas%temperature%qc), & + !next%meas%u%data, real(next%meas%u%qc), & + !next%meas%v%data, real(next%meas%v%qc), & + !next%meas%rh%data, real(next%meas%rh%qc), & + !next%meas%pressure%data, real(next%meas%pressure%qc), & + !obs(i)%ground%precip%data, real(obs(i)%ground%precip%qc) + obs(i)%ground%slp%data, real(obs(i)%ground%slp%qc), & + obs(i)%ground%ref_pres%data, real(obs(i)%ground%ref_pres%qc),& + meas_write%height%data, real(meas_write%height%qc), & + meas_write%temperature%data, real(meas_write%temperature%qc), & + meas_write%u%data, real(meas_write%u%qc), & + meas_write%v%data, real(meas_write%v%qc), & + meas_write%rh%data, real(meas_write%rh%qc), & + meas_write%pressure%data, real(meas_write%pressure%qc), & + obs(i)%ground%precip%data, real(obs(i)%ground%precip%qc) + !BPR END + ENDIF !If marked as sounding or not + ELSE !Not an OBS_DOMAIN file + !BPR BEGIN + !WRITE ( UNIT = unit , FMT = meas_format ) next%meas + WRITE ( UNIT = unit , FMT = meas_format ) meas_write + !BPR END + ENDIF !If OBS_DOMAIN file or not - IF ( .not. OBS_data .and. .not. obs(i)%info%discard ) THEN + !Prepare data for netCDF output files + IF ( ( file_name(1:10) == "qc_obs_use" ) .OR. ( file_name(1:10) == "qc_obs_raw" ) ) THEN + ilev = ilev + 1 + !BPR BEGIN + !Change next%meas to meas_write so we can choose whether to use + !grid- or earth-relative winds - iout_nc = iout_nc + 1 - int_sounding = 0 - IF (is_sounding) int_sounding = 1 - int_bogus = 0 - IF (obs(i)%info%bogus) int_bogus = 1 - int_discard = 0 - IF (obs(i)%info%discard) int_discard = 1 + !press_nc(ilev) = next%meas%pressure%data + !press_qc_nc(ilev) = next%meas%pressure%qc + !z_nc(ilev) = next%meas%height%data + !z_qc_nc(ilev) = next%meas%height%qc + !temp_nc(ilev) = next%meas%temperature%data + !temp_qc_nc(ilev) = next%meas%temperature%qc + !td_nc(ilev) = next%meas%dew_point%data + !td_qc_nc(ilev) = next%meas%dew_point%qc + !spd_nc(ilev) = next%meas%speed%data + !spd_qc_nc(ilev) = next%meas%speed%qc + !dir_nc(ilev) = next%meas%direction%data + !dir_qc_nc(ilev) = next%meas%direction%qc + !u_nc(ilev) = next%meas%u%data + !u_qc_nc(ilev) = next%meas%u%qc + !v_nc(ilev) = next%meas%v%data + !v_qc_nc(ilev) = next%meas%v%qc + !rh_nc(ilev) = next%meas%rh%data + !rh_qc_nc(ilev) = next%meas%rh%qc + press_nc(ilev) = meas_write%pressure%data + press_qc_nc(ilev) = meas_write%pressure%qc + z_nc(ilev) = meas_write%height%data + z_qc_nc(ilev) = meas_write%height%qc + temp_nc(ilev) = meas_write%temperature%data + temp_qc_nc(ilev) = meas_write%temperature%qc + td_nc(ilev) = meas_write%dew_point%data + td_qc_nc(ilev) = meas_write%dew_point%qc + spd_nc(ilev) = meas_write%speed%data + spd_qc_nc(ilev) = meas_write%speed%qc + dir_nc(ilev) = meas_write%direction%data + dir_qc_nc(ilev) = meas_write%direction%qc + u_nc(ilev) = meas_write%u%data + u_qc_nc(ilev) = meas_write%u%qc + v_nc(ilev) = meas_write%v%data + v_qc_nc(ilev) = meas_write%v%qc + rh_nc(ilev) = meas_write%rh%data + rh_qc_nc(ilev) = meas_write%rh%qc + !BPR END + ENDIF !If writing a file for which a netCDF version should be written + ENDIF !If current ob data is to be written to this file + !BPR END REORDERING - start = 1 - count = 1 - start(2) = iout_nc - count(1) = 14 - CALL check (nf_put_vara_text(ncid, 1,start,count,obs(i)%valid_time%date_char(1:14))) - count(1) = 40 - CALL check (nf_put_vara_text(ncid, 2,start,count,obs(i)%location%id)) - CALL check (nf_put_vara_text(ncid, 3,start,count,obs(i)%location%name)) - CALL check (nf_put_vara_text(ncid, 4,start,count,obs(i)%info%platform)) - CALL check (nf_put_vara_text(ncid, 5,start,count,obs(i)%info%source)) + next => next%next + END DO + !BPR BEGIN + !Default use of discard does not appear to account for out_opt variable + !IF ( .not. OBS_data .and. .not. obs(i)%info%discard ) THEN + IF ( .not. OBS_data ) THEN + !Only write the footer if this ob has not been marked as discard=TRUE + !or the output file we are writing is not supposed to throw out obs + IF ( out_opt .EQ. 0 .OR. & + ( out_opt .GT. 0 .AND. .NOT. obs(i)%info%discard ) .OR. & + ( out_opt .LT. 0 .AND. obs(i)%info%discard ) ) THEN + !BPR END - start = 1 - count = 1 - start(1) = iout_nc - CALL check (nf_put_vara_real(ncid, 6,start,count,obs(i)%location%longitude)) - CALL check (nf_put_vara_real(ncid, 7,start,count,obs(i)%location%latitude)) - CALL check (nf_put_vara_real(ncid, 8,start,count,obs(i)%info%elevation)) - CALL check (nf_put_vara_real(ncid, 9,start,count,obs(i)%ground%slp%data)) - CALL check (nf_put_vara_int(ncid,10,start,count,obs(i)%ground%slp%qc)) - CALL check (nf_put_vara_real(ncid,11,start,count,obs(i)%ground%ref_pres%data)) - CALL check (nf_put_vara_int(ncid,12,start,count,obs(i)%ground%ref_pres%qc)) - CALL check (nf_put_vara_int(ncid,13,start,count,true_num_obs)) - CALL check (nf_put_vara_int(ncid, 14,start,count,int_sounding)) - CALL check (nf_put_vara_int(ncid, 15,start,count,int_bogus)) - CALL check (nf_put_vara_int(ncid, 16,start,count,int_discard)) - start = 1 - count = 1 - start(2) = iout_nc - count(1) = true_num_obs - CALL check (nf_put_vara_real(ncid,17,start,count,press_nc)) - CALL check (nf_put_vara_int(ncid,18,start,count,press_qc_nc)) - CALL check (nf_put_vara_real(ncid,19,start,count,z_nc)) - CALL check (nf_put_vara_int(ncid,20,start,count,z_qc_nc)) - CALL check (nf_put_vara_real(ncid,21,start,count,temp_nc)) - CALL check (nf_put_vara_int(ncid,22,start,count,temp_qc_nc)) - CALL check (nf_put_vara_real(ncid,23,start,count,td_nc)) - CALL check (nf_put_vara_int(ncid,24,start,count,td_qc_nc)) - CALL check (nf_put_vara_real(ncid,25,start,count,spd_nc)) - CALL check (nf_put_vara_int(ncid,26,start,count,spd_qc_nc)) - CALL check (nf_put_vara_real(ncid,27,start,count,dir_nc)) - CALL check (nf_put_vara_int(ncid,28,start,count,dir_qc_nc)) - CALL check (nf_put_vara_real(ncid,29,start,count,u_nc)) - CALL check (nf_put_vara_int(ncid,30,start,count,u_qc_nc)) - CALL check (nf_put_vara_real(ncid,31,start,count,v_nc)) - CALL check (nf_put_vara_int(ncid,32,start,count,v_qc_nc)) - CALL check (nf_put_vara_real(ncid,33,start,count,rh_nc)) - CALL check (nf_put_vara_int(ncid,34,start,count,rh_qc_nc)) + iout_nc = iout_nc + 1 + int_sounding = 0 + IF (is_sounding) int_sounding = 1 + int_bogus = 0 + IF (obs(i)%info%bogus) int_bogus = 1 + int_discard = 0 + IF (obs(i)%info%discard) int_discard = 1 + + start = 1 + count = 1 + start(2) = iout_nc + count(1) = 14 + CALL check (nf_put_vara_text(ncid, 1,start,count,obs(i)%valid_time%date_char(1:14))) + count(1) = 40 + CALL check (nf_put_vara_text(ncid, 2,start,count,obs(i)%location%id)) + CALL check (nf_put_vara_text(ncid, 3,start,count,obs(i)%location%name)) + CALL check (nf_put_vara_text(ncid, 4,start,count,obs(i)%info%platform)) + CALL check (nf_put_vara_text(ncid, 5,start,count,obs(i)%info%source)) + + start = 1 + count = 1 + start(1) = iout_nc + CALL check (nf_put_vara_real(ncid, 6,start,count,obs(i)%location%longitude)) + CALL check (nf_put_vara_real(ncid, 7,start,count,obs(i)%location%latitude)) + CALL check (nf_put_vara_real(ncid, 8,start,count,obs(i)%info%elevation)) + CALL check (nf_put_vara_real(ncid, 9,start,count,obs(i)%ground%slp%data)) + CALL check (nf_put_vara_int(ncid,10,start,count,obs(i)%ground%slp%qc)) + CALL check (nf_put_vara_real(ncid,11,start,count,obs(i)%ground%ref_pres%data)) + CALL check (nf_put_vara_int(ncid,12,start,count,obs(i)%ground%ref_pres%qc)) + CALL check (nf_put_vara_int(ncid,13,start,count,true_num_obs)) + CALL check (nf_put_vara_int(ncid, 14,start,count,int_sounding)) + CALL check (nf_put_vara_int(ncid, 15,start,count,int_bogus)) + CALL check (nf_put_vara_int(ncid, 16,start,count,int_discard)) + + start = 1 + count = 1 + start(2) = iout_nc + count(1) = true_num_obs + CALL check (nf_put_vara_real(ncid,17,start,count,press_nc)) + CALL check (nf_put_vara_int(ncid,18,start,count,press_qc_nc)) + CALL check (nf_put_vara_real(ncid,19,start,count,z_nc)) + CALL check (nf_put_vara_int(ncid,20,start,count,z_qc_nc)) + CALL check (nf_put_vara_real(ncid,21,start,count,temp_nc)) + CALL check (nf_put_vara_int(ncid,22,start,count,temp_qc_nc)) + CALL check (nf_put_vara_real(ncid,23,start,count,td_nc)) + CALL check (nf_put_vara_int(ncid,24,start,count,td_qc_nc)) + CALL check (nf_put_vara_real(ncid,25,start,count,spd_nc)) + CALL check (nf_put_vara_int(ncid,26,start,count,spd_qc_nc)) + CALL check (nf_put_vara_real(ncid,27,start,count,dir_nc)) + CALL check (nf_put_vara_int(ncid,28,start,count,dir_qc_nc)) + CALL check (nf_put_vara_real(ncid,29,start,count,u_nc)) + CALL check (nf_put_vara_int(ncid,30,start,count,u_qc_nc)) + CALL check (nf_put_vara_real(ncid,31,start,count,v_nc)) + CALL check (nf_put_vara_int(ncid,32,start,count,v_qc_nc)) + CALL check (nf_put_vara_real(ncid,33,start,count,rh_nc)) + CALL check (nf_put_vara_int(ncid,34,start,count,rh_qc_nc)) - WRITE ( UNIT = unit , FMT = meas_format ) end_meas - WRITE ( UNIT = unit , FMT = end_format ) obs(i)%info%num_vld_fld, & + WRITE ( UNIT = unit , FMT = meas_format ) end_meas + WRITE ( UNIT = unit , FMT = end_format ) obs(i)%info%num_vld_fld, & obs(i)%info%num_error, obs(i)%info%num_warning + !BPR BEGIN + ENDIF !If regarding whether data is to be discarded + !BPR END ENDIF IF ( .NOT. forinput ) & write(unit,*) 'End of measurements for observation ' , i - END IF + END IF END DO @@ -3263,7 +3575,10 @@ !--------------------------------------------------------------------------- SUBROUTINE read_measurements ( file_num , surface , location , bad_data , error , & -height , pressure , iew , jns , levels , map_projection , elevation ) +!BPR BEGIN +!height , pressure , iew , jns , levels , map_projection , elevation ) +height , pressure , slp_x , temperature , iew , jns , levels , map_projection , elevation ) +!BPR END ! This routine reads in 'measurements' at as many levels as there are in ! the report, then stops and returns when an end-of-measurements flag is @@ -3281,6 +3596,10 @@ INTEGER , INTENT ( IN ) :: levels REAL , INTENT ( IN ) , DIMENSION ( levels ) :: pressure INTEGER :: iew , jns , map_projection + !BPR BEGIN + REAL , INTENT ( IN ) , DIMENSION ( jns, iew ) :: slp_x + REAL , INTENT ( IN ) , DIMENSION ( jns , iew , levels ) :: temperature + !BPR END REAL , DIMENSION ( jns , iew , levels ) :: height CHARACTER ( LEN = 32 ) , PARAMETER :: sub_name = 'read_measurements' @@ -3291,6 +3610,10 @@ CHARACTER ( LEN = 40 ) :: location_id , & location_name REAL , INTENT(IN) :: elevation + !BPR BEGIN + LOGICAL :: used_std_atmos_so_recalc + TYPE ( meas_data ) :: meas_earth_winds + !BPR END INCLUDE 'error.inc' INTERFACE @@ -3435,9 +3758,29 @@ ! 1) With the available geopotential height field, we can compute a ! fairly accurate pressure if the height is available from the observation. + !BPR BEGIN + !Check if pressure QC flag indicates that the pressure was calculated + !using standard atmosphere + height + !Only check if the QC flag is not the missing flag + IF ( current%meas%pressure%qc - missing .ne. 0 ) THEN + used_std_atmos_so_recalc = contains_2n ( current%meas%pressure%qc, p_std_atm_and_height ) + ELSE + used_std_atmos_so_recalc = .FALSE. + ENDIF + + !IF ( ( eps_equal ( current%meas%pressure%data , missing_r , 1. ) ) .OR. & + ! ( eps_equal ( current%meas%pressure%data , 0. , 1. ) ) ) THEN + !If pressure is missing, set to zero, or was diagnosed using standard + !atmosphere + height IF ( ( eps_equal ( current%meas%pressure%data , missing_r , 1. ) ) .OR. & - ( eps_equal ( current%meas%pressure%data , 0. , 1. ) ) ) THEN - CALL height_to_pres ( height , pressure , iew , jns , levels , map_projection , & + ( eps_equal ( current%meas%pressure%data , 0. , 1. ) ) .OR. & + ( used_std_atmos_so_recalc ) ) THEN + !BPR END + !BPR BEGIN + ! CALL height_to_pres ( height , pressure , iew , jns , levels , map_projection , & + CALL height_to_pres ( height , pressure , slp_x , temperature , iew , jns , & + levels , map_projection , & + !BPR END location%latitude , location%longitude , current%meas ) END IF @@ -3479,6 +3822,13 @@ CALL diagnostics ( current%meas , location%longitude ) + ! BPR BEGIN + ! Rotate winds to earth-relative (keep non wind fields the same) and place + ! in meas_earth_winds + !TYPE ( meas_data ) :: meas_earth_winds + !CALL grid_to_earth_winds ( current%meas, meas_earth_winds, location%longitude ) + ! BPR END + ! Increment count of measurements correctly read in (these are the levels). meas_count = meas_count + 1 @@ -3549,7 +3899,10 @@ SUBROUTINE read_observations ( file_name , file_num , obs , n_obs , & total_number_of_obs , fatal_if_exceed_max_obs , print_out_found_obs , & -height , pressure , iew , jns , levels , map_projection ) +!BPR BEGIN +!height , pressure , iew , jns , levels , map_projection ) +height , pressure , slp_x , temperature , iew , jns , levels , map_projection ) +!BPR END ! This routine opens file 'file_name' and reads all observations, and ! measurements at all levels, from the file into the 'obs' observation array. @@ -3572,6 +3925,10 @@ INTEGER , INTENT ( IN ) :: levels REAL , INTENT ( IN ) , DIMENSION ( levels ) :: pressure INTEGER :: iew , jns , map_projection + !BPR BEGIN + REAL , INTENT ( IN ) , DIMENSION ( jns, iew ) :: slp_x + REAL , INTENT ( IN ) , DIMENSION ( jns , iew , levels ) :: temperature + !BPR END REAL , DIMENSION ( jns , iew , levels ) :: height CHARACTER ( LEN = 32 ) , PARAMETER :: proc_name = 'read_observations ' @@ -3814,7 +4171,11 @@ NULLIFY ( obs(obs_num)%surface ) CALL read_measurements( file_num , obs(obs_num)%surface , & obs(obs_num)%location , outside_window , error_ret , & - height , pressure , iew , jns , levels , map_projection , obs(obs_num)%info%elevation ) +!BPR BEGIN +! height , pressure , iew , jns , levels , map_projection , obs(obs_num)%info%elevation ) + height , pressure , slp_x , temperature , iew , jns , levels , map_projection , & + obs(obs_num)%info%elevation ) +!BPR END ! An error in the measurements read is handled in a couple of ways. A ! flat out error in the read requires the process to start again (cycle @@ -4102,6 +4463,277 @@ ! ! ---------------------------------------------------------------------------- + +!------------------------------------------------------------------------------ +!BPR BEGIN +! +!Check if sumtocheck decomposed includes numtofind +!If it does contains_2n = .true., elsewise contains_2n = .false. +!sumtocheck is a sum of the form: +!2**x(1)+2**x(2)+... +!where x(n) is composed of integers +LOGICAL FUNCTION contains_2n ( sumtocheck, numtofind ) + INTEGER, INTENT(IN) :: sumtocheck + INTEGER, INTENT(IN) :: numtofind + + INTEGER :: maxn + INTEGER :: curn + INTEGER :: cur2n + INTEGER :: sumnow + + !Determine maximum exponent of n to check + maxn=18 + !In theory could take floor of log base 2 of huge(sumtocheck) + !to determine this but this is overkill here (and may not work due + !to integer/real issues?) + curn=maxn + + sumnow=sumtocheck + sumnowpos: DO WHILE( sumnow > 0 ) + + IF( sumnow < numtofind ) THEN + contains_2n = .FALSE. + RETURN + ELSEIF( sumnow == numtofind ) THEN + contains_2n = .TRUE. + RETURN + ENDIF + + cur2n=2**curn + !PRINT *,'cur2n = ',cur2n,' curn=',curn + !PRINT *,'sumnow =',sumnow + + DO WHILE( cur2n > sumnow ) + curn=curn-1 + cur2n=2**curn + !PRINT *,'LOOP cur2n = ',cur2n,' curn=',curn + END DO + IF( cur2n == numtofind ) THEN + contains_2n = .TRUE. + RETURN + ENDIF + sumnow=sumnow-cur2n + END DO sumnowpos + !Have completely decomposed sumnow and did not find numtofind + contains_2n = .FALSE. + +END FUNCTION contains_2n + +!------------------------------------------------------------------------------ + +SUBROUTINE grid_to_earth_winds ( meas_grid_winds, meas_earth_winds , longitude ) + +! Copy meas_grid_winds to meas_earth_winds, but rotate winds from grid to earth +! Also set meas_earth_winds wind speed and direction since they are not necessarily +! consistent with u/v + + IMPLICIT NONE + + TYPE ( meas_data ), INTENT(INOUT) :: meas_grid_winds + TYPE ( meas_data ), INTENT(OUT) :: meas_earth_winds + REAL, INTENT(IN) :: longitude + + REAL :: dir_map , lon_dif + + INTEGER :: iew_map , jns_map + + REAL , PARAMETER :: piover180 = 3.14159265358 / 180. + REAL , PARAMETER :: one80overpi = 180. / 3.14159265358 + INTEGER :: uv_to_dirspd_status + + INCLUDE 'proc_get_info_header.inc' + INCLUDE 'map.inc' + INCLUDE 'error.inc' + INTERFACE + INCLUDE 'error.int' + END INTERFACE + + meas_earth_winds = meas_grid_winds + + ! The wind components are computed in two steps. First we need the + ! meteorological u and v from the speed and direction. Second, those + ! values of u and v are rotated to the model grid. + + IF ( ( eps_equal ( meas_grid_winds%u%data , missing_r , 1. ) ) .OR. & + ( eps_equal ( meas_grid_winds%v%data , missing_r , 1. ) ) ) THEN + meas_earth_winds%u%data = missing_r + meas_earth_winds%u%qc = missing + meas_earth_winds%v%data = missing_r + meas_earth_winds%v%qc = missing + meas_earth_winds%speed%data = missing_r + meas_earth_winds%speed%qc = missing + meas_earth_winds%direction%data = missing_r + meas_earth_winds%direction%qc = missing + ELSE + !Calculate the grid-relative wind direction / speed + !dir_map = one80overpi * acos(-1.0*meas_grid_winds%v%data/meas_earth_winds%speed%data) + CALL uv_to_dirspd(meas_grid_winds%u%data,meas_grid_winds%v%data,meas_grid_winds%direction%data,& + meas_grid_winds%speed%data,uv_to_dirspd_status) + IF(uv_to_dirspd_status.lt.0) THEN + PRINT *,'ERROR: Error in uv_to_dirspd' + STOP "ERROR in uv_to_dirspd" + ENDIF + meas_grid_winds%direction%qc = MAX( meas_earth_winds%u%qc, meas_earth_winds%v%qc ) + meas_grid_winds%speed%qc = meas_grid_winds%direction%qc + + !Set the earth-relative wind speed equal to the grid-relative value since + !they are equal + meas_earth_winds%speed%data = meas_grid_winds%speed%data + meas_earth_winds%speed%qc = meas_grid_winds%speed%qc + + !If the grid-relative wind speed is basically zero then mark direction as + !missing + !Also mark earth-relative u/v as zero + IF( eps_equal ( meas_grid_winds%speed%data, 0.0 , 0.00001 ) ) THEN + meas_grid_winds%direction%data = missing_r + meas_grid_winds%direction%qc= missing + meas_earth_winds%direction%data = missing_r + meas_earth_winds%direction%qc= missing + meas_earth_winds%u%data = 0.0 + meas_earth_winds%v%data = 0.0 + ELSE + !Calculate the earth-relative wind direction + meas_earth_winds%direction%data = meas_grid_winds%direction%data + (longitude - lon_center ) * & + cone_factor * sign(1.0, lat_center) + meas_earth_winds%direction%qc = meas_earth_winds%speed%qc + !Calculate the earth-relative u/v + meas_earth_winds%u%data = -1. * meas_earth_winds%speed%data * & + SIN ( meas_earth_winds%direction%data * piover180 ) + meas_earth_winds%v%data = -1. * meas_earth_winds%speed%data * & + COS ( meas_earth_winds%direction%data * piover180 ) + ENDIF + meas_earth_winds%u%qc = meas_earth_winds%speed%qc + meas_earth_winds%v%qc = meas_earth_winds%speed%qc + + + !lon_dif = longitude - lon_center + !if ( lon_dif .gt. 180. ) lon_dif = lon_dif - 360. + !if ( lon_dif .lt. -180. ) lon_dif = lon_dif + 360. + !dir_map = new%direction%data - ( lon_dif ) * cone_factor * SIGN ( 1. , lat_center ) + !if ( dir_map .GT. 360. ) dir_map = dir_map - 360. + !if ( dir_map .LT. 0. ) dir_map = 360 + dir_map + !new%u%data = -1. * new%speed%data * SIN ( dir_map * piover180 ) + !new%v%data = -1. * new%speed%data * COS ( dir_map * piover180 ) + !new%u%qc = new%direction%qc + !new%v%qc = new%direction%qc + END IF + +END SUBROUTINE grid_to_earth_winds + +!Calculate wind direction and wind speed from u and v +SUBROUTINE uv_to_dirspd(u,v,direction,speed,uv_to_dirspd_status) + + IMPLICIT NONE + + REAL, INTENT(IN) :: u,v + REAL, INTENT(OUT) :: direction,speed + INTEGER, INTENT(OUT) :: uv_to_dirspd_status + REAL :: pi + REAL :: absu,absv !Absolute value of wind components + + uv_to_dirspd_status = 0 + + pi = acos(-1.0) + + speed = (u*u+v*v)**0.5 + absu=abs(u) + absv=abs(v) + + if((u.ge.0).and.(v.gt.0)) then !Quadrant 1 + direction = atan(absu/absv) + elseif((u.gt.0).and.(v.le.0)) then !Quadrant 2 + direction = atan(absv/absu) + pi/2.0 + elseif((u.le.0).and.(v.lt.0)) then !Quadrant 3 + direction = atan(absu/absv) + pi + elseif((u.lt.0).and.(v.ge.0)) then !Quadrant 4 + direction = atan(absv/absu) + 3*pi/2.0 + elseif((u.eq.0).and.(v.eq.0)) then !Calm + direction = 0 + endif + + direction = 360.0*direction/(2*pi) + + !Now convert from real vector to meteorological wind vector + if((u.eq.0).and.(v.eq.0)) then !Calm + direction=0.0 + else + direction=direction+180.0 + endif + + if(direction.ge.360.0) then + direction=direction-360.0 + endif + +END SUBROUTINE uv_to_dirspd + +!BPR END + +!BPR BEGIN +!Sort list using comb sort algorithm +SUBROUTINE COMBSORT_INT(array_to_sort,order_index) + + IMPLICIT NONE + + INTEGER, PARAMETER :: INT_14DIGITS = selected_int_kind ( 14 ) + !The array to be sorted + INTEGER (kind=INT_14DIGITS) , INTENT(INOUT), DIMENSION(:) :: array_to_sort + !An array holding indices to translate between the array_to_sort + !passed in, and the one passed out + !In other words, if order_index(5)=8 it means that what was in + !array_to_sort(8) on input is now in array_to_sort(5) + INTEGER, INTENT(INOUT), DIMENSION(:) :: order_index + + + INTEGER :: num_vals_to_sort + INTEGER :: gap + !Gap shrink factor + REAL, PARAMETER :: shrink = 1.3 + INTEGER :: counter + LOGICAL :: swapped + INTEGER :: element_to_sort_temp + INTEGER :: order_index_temp + + !Determine how many values we have to sort + num_vals_to_sort = size(array_to_sort) + + !Initialize gap to number of elements to sort + gap = num_vals_to_sort + + !Initialize order index so that order_index(n) = n + !This indicates that before sorting the n'th element in array_to_sort + !when it was passed in is now the n'th element in array_to_sort + order_index(:) = (/ (counter,counter=1,num_vals_to_sort) /) + + !Set swapped to TRUE so that we can enter the loop + swapped = .true. + DO WHILE( ( gap .ne. 1 ) .OR. ( swapped ) ) + swapped = .false. + !Update the gap value + gap = gap / shrink + IF( gap .lt. 1 ) THEN + gap = 1 + END IF + + DO counter = 1, num_vals_to_sort - gap + IF( array_to_sort(counter) .gt. array_to_sort(counter + gap) ) THEN + !Swap [counter] and [counter + gap] + element_to_sort_temp = array_to_sort( counter ) + order_index_temp = order_index ( counter ) + array_to_sort( counter ) = array_to_sort( counter + gap ) + order_index ( counter ) = order_index ( counter + gap ) + array_to_sort( counter + gap ) = element_to_sort_temp + order_index ( counter + gap ) = order_index_temp + swapped = .TRUE. + END IF + END DO + END DO + + +END SUBROUTINE COMBSORT_INT +!BPR END + +!^L +! ---------------------------------------------------------------------------- subroutine check(status) INCLUDE 'netcdf.inc' Index: src/plot_soundings.F90 =================================================================== --- src/plot_soundings.F90 (revision 374) +++ src/plot_soundings.F90 (working copy) @@ -95,7 +95,10 @@ n_times = idiff / interval ! Build metoa file - WRITE(oa_file, '("metoa_em.d",I2.2".",A19,".nc")') grid_id, start_date + !BPR BEGIN + !WRITE(oa_file, '("metoa_em.d",I2.2".",A19,".nc")') grid_id, start_date + WRITE(oa_file, '("metoa_em.d",I2.2,".",A19,".nc")') grid_id, start_date + !BPR END IF ( read_metoa ) THEN print*," Attempting to open file: ", trim(oa_file) Index: src/proc_final_analysis.F90 =================================================================== --- src/proc_final_analysis.F90 (revision 374) +++ src/proc_final_analysis.F90 (working copy) @@ -11,7 +11,10 @@ max_error_t , max_error_uv , & max_error_z , max_error_p , & buddy_weight , date_char , root_filename , oa_3D_option , & -intf4d , lagtem , oa_type , oa_3D_type , rad_influence ) +!BPR BEGIN +!intf4d , lagtem , oa_type , oa_3D_type , rad_influence ) +intf4d , lagtem , oa_type , oa_3D_type , rad_influence, oa_psfc ) +!BPR END ! This routine is a driver for the required utilities to output the ! final analysis of this program. The input values are the objectively @@ -54,7 +57,11 @@ INTEGER :: oa_3D_option CHARACTER *(*) :: oa_type , oa_3D_type INTEGER , DIMENSION(10) :: rad_influence + !BPR BEGIN + LOGICAL :: oa_psfc + !BPR END + ! Temporary holding arrays for the header information. What is in the header ! on input is the data from the first guess file. To output the data, we need ! to store the final analysis header in the bhi, bhr, arrays. So that @@ -122,7 +129,10 @@ latitude_x , longitude_x , latitude_d , longitude_d , & slp_x , slp_C , sst , tobbox , odis , & iew_alloc , jns_alloc , kbu_alloc , iewd , jnsd , date_char, print_analysis , & - oa_type , oa_3D_type , oa_3D_option , rad_influence ) +!BPR BEGIN +! oa_type , oa_3D_type , oa_3D_option , rad_influence ) + oa_type , oa_3D_type , oa_3D_option , rad_influence , oa_psfc ) +!BPR END ELSE CALL make_date ( current_date_8 , current_time_6 , fdda_date_24 ) slp_x(:jns_alloc-1,:iew_alloc-1) = slp_x(:jns_alloc-1,:iew_alloc-1) * 100. Index: src/proc_final_analysis.int =================================================================== --- src/proc_final_analysis.int (revision 374) +++ src/proc_final_analysis.int (working copy) @@ -10,7 +10,10 @@ max_error_t , max_error_uv , & max_error_z , max_error_p , & buddy_weight , date_char , root_filename , oa_3D_option , & -intf4d , lagtem , oa_type , oa_3D_type , rad_influence ) +!BPR BEGIN +!intf4d , lagtem , oa_type , oa_3D_type , rad_influence ) +intf4d , lagtem , oa_type , oa_3D_type , rad_influence , oa_psfc ) +!BPR END CHARACTER *(*) , INTENT ( IN ) :: filename CHARACTER *(*) , INTENT ( INOUT ) :: filename_out @@ -39,6 +42,9 @@ INTEGER :: oa_3D_option CHARACTER *(*) :: oa_type , oa_3D_type INTEGER , DIMENSION(10) :: rad_influence + !BPR BEGIN + LOGICAL :: oa_psfc + !BPR END INCLUDE 'proc_get_info_header.inc' END SUBROUTINE proc_final_analysis Index: src/proc_namelist.F90 =================================================================== --- src/proc_namelist.F90 (revision 374) +++ src/proc_namelist.F90 (working copy) @@ -62,11 +62,43 @@ remove_data_above_qc_flag = 200000 remove_unverified_data = .FALSE. + !BPR BEGIN + ! Default to using very large tolerances for dewpoint QC checks so as to + ! default to the the original Obsgrid methodology of not doing a separate QC + ! on dewpoint in addition to the QC already performed for RH + max_error_dewpoint = 99999.9 + max_buddy_dewpoint = 99999.9 + + ! Default to original Obsgrid behavior of not doing QC on surface pressure + qc_psfc = .FALSE. + + ! Default to original Obsgrid behavior of generally only dealing with + ! single-level above-surface observations if they fall on a pressure level + ! with gridded data. + use_p_tolerance_one_lev = .FALSE. + max_p_tolerance_one_lev_qc = 1 + max_p_tolerance_one_lev_oa = 1 + + !BPR END + ! Initialize the array of radius of influence scans. This permits an easy way ! to determine the number of requested scans without an additional input value. radius_influence = -1 + !BPR BEGIN + ! Default to the original OBSGRID setup of not scaling RH decreases in the + ! Cressman analysis + scale_cressman_rh_decreases = .FALSE. + + ! Default to the original OBSGRID setup of using the same radii of influence + ! for applying surface obs to the analysis as are used for non-surface obs + radius_influence_sfc_mult = 1.0 + + ! Default to original Obsgrid behavior of not doing objective analysis on surface pressure + oa_psfc = .FALSE. + !BPR END + ! Initialize the array of observation file names to 'null'. obs_filename = 'null & Index: src/proc_oa.F90 =================================================================== --- src/proc_oa.F90 (revision 374) +++ src/proc_oa.F90 (working copy) @@ -1,11 +1,19 @@ ! ------------------------------------------------------------------------------ SUBROUTINE proc_oa ( t , u , v , rh , slp_x , & -pressure , & +!BPR BEGIN +!Add 3D pressure of first guess to fields passed in +!pressure is 2D pressure with surface pressure just set to 1001 +!pressure , & +pressure , pres , & +!BPR END iew_alloc , jns_alloc , kbu_alloc , & date , time , fdda_loop , mqd_count , mqd_abs_min , & total_numobs , num_obs_found , total_dups , & map_projection , obs , dxd , lat_center , & +!BPR BEGIN +use_p_tolerance_one_lev , & +!BPR END print_oa , print_found_obs , print_obs_files , use_first_guess , & smooth_type , smooth_sfc_wind , & smooth_sfc_temp , smooth_sfc_rh , & @@ -15,12 +23,19 @@ mqd_minimum_num_obs , mqd_maximum_num_obs , & oa_max_switch , radius_influence , & oa_min_switch , oa_3D_option , & -grid_id ) +!BPR BEGIN +!grid_id ) +grid_id, terrain, h, scale_cressman_rh_decreases, radius_influence_sfc_mult, & +oa_psfc, max_p_tolerance_one_lev_oa ) +!BPR END ! This routine is a driver routine for objective analysis. USE obj_analysis USE observation +!BPR BEGIN + USE final_analysis, only : mixing_ratio +!BPR END IMPLICIT NONE @@ -39,7 +54,15 @@ print_obs_files , & use_first_guess , & oa_min_switch , & - oa_max_switch + !BPR BEGIN + ! oa_max_switch + oa_max_switch , & + scale_cressman_rh_decreases , & + oa_psfc + REAL , INTENT ( IN ) :: radius_influence_sfc_mult + INTEGER , INTENT ( IN ) :: max_p_tolerance_one_lev_oa + LOGICAL , INTENT ( IN ) :: use_p_tolerance_one_lev + !BPR END INTEGER :: oa_3D_option INTEGER :: date , & time , & @@ -72,7 +95,10 @@ INTEGER , DIMENSION ( num_obs_found ) :: array_index , & qc_flag CHARACTER ( LEN = 8 ),DIMENSION ( num_obs_found ) :: station_id - INTEGER , PARAMETER :: num_var = 5 + !BPR BEGIN + !INTEGER , PARAMETER :: num_var = 5 + INTEGER , PARAMETER :: num_var = 6 + !BPR END CHARACTER ( LEN = 8 ) , DIMENSION ( num_var) :: name INTEGER :: passes , & num_scan @@ -82,6 +108,15 @@ passes_upper REAL , DIMENSION ( num_var) :: scale REAL , DIMENSION ( iew_alloc , jns_alloc ) :: dum2d + !BPR BEGIN + REAL , DIMENSION ( iew_alloc-1 , jns_alloc-1 ) :: dum2d_fg, t_first_guess, t_analysis, pressure_2d + + REAL :: pressure_t + REAL, DIMENSION(jns_alloc,iew_alloc,kbu_alloc) :: qv3d_first_guess_yx + REAL, DIMENSION(jns_alloc,iew_alloc) :: psfc_first_guess_model_terrain_yx + REAL, DIMENSION(iew_alloc,jns_alloc) :: psfc_first_guess_model_terrain_xy + REAL , DIMENSION ( num_obs_found ) :: fg_value + !BPR END REAL , DIMENSION ( total_numobs) :: diff REAL :: dxob , & dyob , & @@ -105,10 +140,54 @@ min_mqd , max_mqd , & mqd_count , mqd_abs_min , & test_count + !BPR BEGIN + REAL :: radius_influence_scale_factor + REAL :: radius_influence_sfc_min_cells, radius_influence_sfc_max_cells + REAL :: radius_influence_scaled, radius_influence_scaled_scan1 + LOGICAL, DIMENSION ( kbu_alloc,max_scan ) :: roi_not_printed_yet + REAL, DIMENSION(jns_alloc,iew_alloc) :: slp_x_pa - skip_to_cressman = .TRUE. + INTEGER :: curx, cury + INTEGER :: num_last_3d_variable !Index to last of variables under consideration + !that is 3D + !Maximum difference allowed between pressure of ob and pressure of background + !data it will be incorporated with for the objective analysis (Pa) + !Default is 1 + INTEGER :: request_p_diff + + !If the user wants to allow a tolerance between an obs pressure and the + !pressure level it can be used for an objective anlysis on, then use the user-specified tolerance. + !If not, then use a tolerance of 1 Pa, which is effectively no tolerance. + + IF( use_p_tolerance_one_lev ) THEN + request_p_diff = max_p_tolerance_one_lev_oa + ELSE + request_p_diff = 1 + END IF + + roi_not_printed_yet(:,:) = .TRUE. + !BPR END + + !BPR BEGIN + !It seems that we should default to NOT skipping to Cressman unless we find + !evidence that we need to + !In most cases what we initialize this to should not matter since: + !1) If the user chose MQD then the code will set skip_to_cressman as appropriate + !2) If the user chose Cressman, then skipping to Cressman will be consistent + ! with what the user chose + !However, this switch seems to be intended to be used to force a switch to + !Cressman when needed and thus it is most straightforward to only set it to + !TRUE when the situtation indicates that it is needed. Also, if there are + !any other analysis methodologies that do not explicitly set skip_to_cressman, + !those methodologies will be skipped over due to skip_to_cressman being + !defaulted to true. This situation may occur for the choice of "none" added in this + !version of Obsgrid indicating that no objective analysis should be done. + !skip_to_cressman = .TRUE. + skip_to_cressman = .FALSE. + !BPR END + num_mqd_uu = 0 num_mqd_vv= 0 num_mqd_tt = 0 @@ -120,12 +199,28 @@ ! Names of all of the variables to objectively analyze. Also, define each of the ! horizontal staggerings for the variables. - name = (/ 'UU ' , 'VV ' , 'TT ' , 'RH ' , 'PMSL ' /) - crsdot = (/ 1 , 1 , 1 , 1 , 1 /) - scale = (/ 1. , 1. , 1. , 1. , 100. /) - passes_sfc = (/ smooth_sfc_wind , smooth_sfc_wind , smooth_sfc_temp , smooth_sfc_rh , smooth_sfc_slp /) - passes_upper = (/ smooth_upper_wind , smooth_upper_wind , smooth_upper_temp , smooth_upper_rh , 0 /) + ! BPR BEGIN + ! Add surface pressure + name = (/ 'UU ' , 'VV ' , 'TT ' , 'RH ' ,& + 'PMSL ' , 'PSFC ' /) + crsdot = (/ 1 , 1 , 1 , 1 ,& + 1 , 1 /) + scale = (/ 1. , 1. , 1. , 1. ,& + 100. , 1. /) + passes_sfc = (/ smooth_sfc_wind , smooth_sfc_wind , smooth_sfc_temp , smooth_sfc_rh ,& + smooth_sfc_slp , smooth_sfc_slp /) + passes_upper = (/ smooth_upper_wind , smooth_upper_wind , smooth_upper_temp , smooth_upper_rh ,& + 0 , 0 /) + !name = (/ 'UU ' , 'VV ' , 'TT ' , 'RH ' , 'PMSL ' /) + !crsdot = (/ 1 , 1 , 1 , 1 , 1 /) + !scale = (/ 1. , 1. , 1. , 1. , 100. /) + !passes_sfc = (/ smooth_sfc_wind , smooth_sfc_wind , smooth_sfc_temp , smooth_sfc_rh , smooth_sfc_slp /) + !passes_upper = (/ smooth_upper_wind , smooth_upper_wind , smooth_upper_temp , smooth_upper_rh , 0 /) + !Index to last of the 3D variables + num_last_3d_variable = 4 + !BPR END + ! This is where we are outputting the information about the observations ! that went into this analysis, for this level, this variable, this time. ! First, check to see if this UNIT has been OPENed yet. The first thing @@ -171,12 +266,19 @@ IF ( oa_3D_type .EQ. 'MQD' .AND. fdda_loop == 1 ) THEN skip_to_cressman = .FALSE. - DO ivar = 1 , num_var-1 + !BPR BEGIN + !DO ivar = 1 , num_var-1 + !Loop over the 3D variables + DO ivar = 1 , num_last_3d_variable + !BPR END first_time = .TRUE. test_count = 0 DO kp = 2 , kbu_alloc CALL proc_ob_access ( 'use', name(ivar) , print_found_obs , & - pressure(kp) , date , time , 1 , & + !BPR BEGIN + !pressure(kp) , date , time , 1 , & + pressure(kp) , date , time , request_p_diff , & + !BPR END MIN ( fails_error_max , fails_buddy_check ) , & num_obs_found , num_obs_pass , obs , & iew_alloc , jns_alloc , kbu_alloc , & @@ -209,7 +311,7 @@ IF ( skip_to_cressman .AND. oa_3D_option == 0 .AND. oa_min_switch ) THEN write (*,*) write (*,'("###########################################################################")') - write (*,'(" Too few observations have been found to do MQD scheme for all levels ")') + write (*,'(" Too few observations have been found to do MQD scheme for all levels ")') write (*,'(" Your options are: ")') write (*,'(" Set oa_3D_option = 1 (if any upper-air level have too few observations,")') write (*,'(" revert to Cressman for all upper-air levels for that time) ")') @@ -227,6 +329,26 @@ END IF + ! BPR BEGIN + ! Calculate the 3D QV field of the first guess field + ! To calculate the surface QV, this calculates a surface pressure based on the + ! first-guess file which includes the terrain field for the target model + ! configuration. This means the surface pressure (psfc_first_guess_model_terrain_*) + ! is not merely that of the coarse grid model, but a potentially much more + ! finescale product depending on the difference in horizontal resolution + ! between the coarse grid model and the target model. + + ! Convert from hPa (mb) to Pa because that is what mixing_ratio expects + slp_x_pa = slp_x * 100.0 + ! Initialize 3D QV and surface pressure to zero + qv3d_first_guess_yx=0.0 + psfc_first_guess_model_terrain_yx=0.0 + CALL mixing_ratio ( rh , t , h , terrain , slp_x_pa , pressure , iew_alloc , jns_alloc , & + kbu_alloc , qv3d_first_guess_yx , psfc_first_guess_model_terrain_yx) + !Change from (y,x) to (x,y) + CALL yx2xy( psfc_first_guess_model_terrain_yx, jns_alloc, iew_alloc, psfc_first_guess_model_terrain_xy) + ! BPR END + ! Loop through all analysis levels (remember that level kp=1 is the ! surface value of the 3-D field). @@ -249,11 +371,17 @@ ! data going into this routine are still on the (y,x) orientation, but ! the 2-D fields on output (u_banana and v_banana) are oriented (x,y). - CALL get_background_for_oa ( t , u , v , rh , slp_x , & + !BPR BEGIN + !CALL get_background_for_oa ( t , u , v , rh , slp_x , & + CALL get_background_for_oa ( t , u , v , rh , slp_x , pres , & + !BPR END u_banana , kp , name(1) , & iew_alloc , jns_alloc , kbu_alloc ) - CALL get_background_for_oa ( t , u , v , rh , slp_x , & + !BPR BEGIN + !CALL get_background_for_oa ( t , u , v , rh , slp_x , & + CALL get_background_for_oa ( t , u , v , rh , slp_x , pres , & + !BPR END v_banana , kp , name(2) , & iew_alloc , jns_alloc , kbu_alloc ) @@ -261,11 +389,22 @@ variable_loop : DO ivar = 1 , num_var + !BPR BEGIN + !If user chose not to do objective analysis of surface pressure then go + !to the next variable + IF ( (.NOT. oa_psfc) .AND.( name(ivar)(1:8) .EQ. 'PSFC ' ) ) THEN + CYCLE variable_loop + ENDIF + !BPR END + ! For sea level pressure, the analysis is only at one level, so cycle the loop - ! (which should end this variable loop, but that is not important). - - - IF ( ( name(ivar)(1:8) .EQ. 'PMSL ' ) .AND. ( kp .GT. 1 ) ) THEN + ! BPR BEGIN + !! IGNORE THIS COMMENT AS IT IS NO LONGER TRUE -> (which should end this variable loop, but that is not important). + ! Same for surface pressure + !IF ( ( name(ivar)(1:8) .EQ. 'PMSL ' ) .AND. ( kp .GT. 1 ) ) THEN + IF ( ( ( name(ivar)(1:8) .EQ. 'PSFC ' ) .OR. & + ( name(ivar)(1:8) .EQ. 'PMSL ' ) ) .AND. ( kp .GT. 1 ) ) THEN + ! BPR END CYCLE variable_loop END IF @@ -283,7 +422,10 @@ ! Get obs for printing and plotting. CALL proc_ob_access ( 'use', name(ivar) , print_found_obs , & - pressure(kp) , date , time , 1 , & + !BPR BEGIN + !pressure(kp) , date , time , 1 , & + pressure(kp) , date , time , request_p_diff , & + !BPR END 200000 , & num_obs_found , num_obs_pass , obs , & iew_alloc , jns_alloc , kbu_alloc , & @@ -303,7 +445,7 @@ !WRITE ( UNIT =74 , FMT = '( 3x,a8,3x,i6,3x,i5,3x,a8,3x,g13.6,3x,16x,2(f7.2,3x),i7 )' ) & WRITE ( UNIT =74 , FMT = '( 3x,a8,3x,i6,3x,i5,3x,a8,3x,2(g13.6,3x),2(f7.2,3x),i7 )' ) & name(ivar) , NINT ( pressure(kp) ) , num , station_id(num) , & - obs_value(num) , diff(num) , xob(num) , yob(num) , qc_flag(num) + obs_value(num) , diff(num) , xob(num) , yob(num) , qc_flag(num) !obs_value(num) , xob(num) , yob(num) , qc_flag(num) END DO station_loop_74 END IF @@ -312,9 +454,11 @@ ! Obtain observations for kp level and for variable ivar for this ! time period. - CALL proc_ob_access ( 'use', name(ivar) , print_found_obs , & - pressure(kp) , date , time , 1 , & + !BPR BEGIN + !pressure(kp) , date , time , 1 , & + pressure(kp) , date , time , request_p_diff , & + !BPR END MIN ( fails_error_max , fails_buddy_check ) , & num_obs_found , num_obs_pass , obs , & iew_alloc , jns_alloc , kbu_alloc , & @@ -342,9 +486,22 @@ ! data going into this routine are still on the (y,x) orientation, but ! the 2-D field on output (dum2d) is oriented (x,y). - CALL get_background_for_oa ( t , u , v , rh , slp_x , & + !BPR BEGIN + !CALL get_background_for_oa ( t , u , v , rh , slp_x , & + CALL get_background_for_oa ( t , u , v , rh , slp_x , pres , & + !BPR END dum2d , kp , name(ivar) , & iew_alloc , jns_alloc , kbu_alloc ) + + !BPR BEGIN + !Save a copy of the first guess field before it is changed by the + !analysis + !Note that although dum2d is dimensioned iew_alloc,jns_alloc, for + !temperature and qv the values at (:,jns_alloc) and (iew_alloc,:) + !are meaningless since those points are just there to allow u/v + !points to fit + dum2d_fg = dum2d(1:iew_alloc-1,1:jns_alloc-1) + !BPR END ! Calculate the difference between the observations and first guess fields. ! Estimate the values of the first guess at each of the observation location. @@ -352,6 +509,10 @@ ! "x" and a linear "y" weighting based on distance. We are either using the ! first-guess to create perturbations, or we can do an analysis from the ! observations only. + +!BPR BEGIN + fg_value(:) = -9999999.0 +!BPR END IF ( use_first_guess ) THEN station_loop_1fg : DO num = 1, num_obs_pass @@ -365,6 +526,9 @@ dxob * ( ( 1.-dyob ) * dum2d ( iob+1 , job ) + & dyob * dum2d ( iob+1 , job+1 ) ) diff ( num ) = obs_value(num) / scale(ivar) - aob +!BPR BEGIN + fg_value(num) = aob +!BPR END IF ( ( name(ivar)(1:8) .EQ. 'PMSL ' ) ) THEN if ( iob >= 11 .and. iob <= 29 ) then if ( job >= 11 .and. job <= 29 ) then @@ -420,6 +584,13 @@ ! obs_value(num) , diff(num) , xob(num) , yob(num) , qc_flag(num) ! END DO station_loop_4 !END IF + + !BPR BEGIN + !Ensure that oa_type stays as None if that is what the user chose + IF ( oa_type .EQ. 'None' ) THEN + oa_type_tmp = 'None' + ENDIF + !BPR END IF ( ( oa_type_tmp .EQ. 'MQD' ) .AND. & ( pressure(kp) .EQ. 1001. ) .AND. & @@ -449,10 +620,17 @@ dum2d , iew_alloc , jns_alloc , & crsdot(ivar) , name(ivar) , passes , smooth_type , use_first_guess ) - - ELSE IF ( ( oa_type_tmp .EQ. 'Cressman' ) .OR. & - ( ( oa_min_switch ) .AND. ( num_obs_pass .LT. mqd_minimum_num_obs ) ) .OR. & - ( ( oa_max_switch ) .AND. ( num_obs_pass .GT. mqd_maximum_num_obs ) ) ) THEN + + !BPR BEGIN + !Code to ensure oa_type stays as None if that is what the user chose + !ELSE IF ( ( oa_type_tmp .EQ. 'Cressman' ) .OR. & + ! ( ( oa_min_switch ) .AND. ( num_obs_pass .LT. mqd_minimum_num_obs ) ) .OR. & + ! ( ( oa_max_switch ) .AND. ( num_obs_pass .GT. mqd_maximum_num_obs ) ) ) THEN + ELSE IF ( ( ( oa_type_tmp .EQ. 'Cressman' ) .OR. & + ( ( oa_min_switch ) .AND. ( num_obs_pass .LT. mqd_minimum_num_obs ) ) .OR. & + ( ( oa_max_switch ) .AND. ( num_obs_pass .GT. mqd_maximum_num_obs ) ) ) .AND. & + ( oa_type .NE. 'None' ) ) THEN + !BPR END if ( ivar == 1 .and. pressure(kp) .eq. 1001 ) print*," " if ( ivar == 1 .and. pressure(kp) .eq. 1001 ) print*," Doing Cressman for surface data" @@ -471,12 +649,102 @@ ! the scheme is called repeatedly, and the "final" analysis is ! used as input for the new difference field that is objectively ! analyzed. + + !BPR BEGIN + !IF (( name(ivar)(1:8) .EQ. 'RH ' ).AND.(kp.eq.9)) THEN + ! DO cury=1,jns_alloc + ! WRITE(110+num_scan,*) dum2d(:,cury) + ! ENDDO + !ENDIF + !BPR END + + !BPR BEGIN + !If we are dealing with the surface then use a smaller radius of influence + IF(abs(pressure(kp)-1001.0).lt.0.0001) THEN + radius_influence_scale_factor = radius_influence_sfc_mult + ELSE + radius_influence_scale_factor = 1.00 + ENDIF + radius_influence_scaled = radius_influence(num_scan)*radius_influence_scale_factor + radius_influence_scaled_scan1 = radius_influence(1)*radius_influence_scale_factor + + !IF the user chose to use reduced radii of influence for + !surface obs compared to other obs + IF(radius_influence_sfc_mult.lt.0.99999) THEN + !Ensure that for the surface the radius of influence in terms of model grid + !cells stays between radius_influence_sfc_min_cells and radius_influence_sfc_max_cells + radius_influence_sfc_max_cells = 100.0; + radius_influence_sfc_min_cells = 4.5; + IF(abs(pressure(kp)-1001.0).lt.0.0001) THEN + + !Ensure that the radius of influence on the first scan + !does not exceed a certain number of grid points. If it does, then adjust the radius + !on each scan by the amount necessary so that the first scan has a radius equal to + !the maximum allowed number of grid points. + !This is to avoid everything being washed out on finer + !resolution domains even when there is dense data + IF(radius_influence_scaled_scan1.gt.radius_influence_sfc_max_cells) THEN + radius_influence_scaled = radius_influence_scaled * & + ( radius_influence_sfc_max_cells / radius_influence_scaled_scan1 ) + ENDIF + !Ensure that the radius of influence on each scan does not go + !below a certain number of grid points. If it does, then + !skip the remaining scans. This is to minimize "spots" in the + !data. + IF(radius_influence_scaled.lt.radius_influence_sfc_min_cells) THEN + IF(num_scan.eq.1) THEN + PRINT *,'ERROR: The effective radius of influence for the first Cressman scan ',& + 'incorporating surface obs into the analysis was too few grid points so ',& + 'ZERO Cressman scans would have been used and thus the surface data would ',& + 'not be incorporated into the analysis. The radius of influence for the first',& + ' Cressman scan was set to ',radius_influence_scaled,' model grid points, but ',& + 'the minimum value is ',radius_influence_sfc_min_cells,'.' + STOP 'ERROR: Effective radius of influence for the first Cressman scan of surface obs too small to continue.' + ENDIF + EXIT multi_scan + ENDIF + ENDIF + ENDIF + IF( roi_not_printed_yet(kp,num_scan) ) THEN + IF(abs(pressure(kp)-1001.0).lt.0.0001) THEN +!BPR BEGIN +! WRITE(*,'(A,A,I3,A,F6.1,A,F7.2,A,F6.1,A)'), 'ROI for the surface,',& + WRITE(*,'(A,A,I3,A,F6.1,A,F7.2,A,F6.1,A)') 'ROI for the surface,',& +!BPR END + ' Cressman scan number ',num_scan,' is ',radius_influence_scaled,& + ' grid points or ',radius_influence_scaled*dxd,'km (',dxd,' km spacing).' + ELSE +!BPR BEGIN +! WRITE(*,'(A,F8.2,A,I3,A,F6.1,A,F7.2,A,F6.1,A)'), 'ROI for ',pressure(kp),& + WRITE(*,'(A,F8.2,A,I3,A,F6.1,A,F7.2,A,F6.1,A)') 'ROI for ',pressure(kp),& +!BPR END + ' hPa, Cressman scan number ',num_scan,' is ',radius_influence_scaled,& + ' grid points or ',radius_influence_scaled*dxd,'km (',dxd,' km spacing).' + ENDIF + roi_not_printed_yet(kp,num_scan) = .FALSE. + ENDIF + !BPR END + CALL cressman ( obs_value , diff , xob , yob , num_obs_pass , & dum2d , iew_alloc , jns_alloc , & - crsdot(ivar) , name(ivar) , radius_influence(num_scan) , & +!BPR BEGIN + crsdot(ivar) , name(ivar) , nint(radius_influence_scaled) , & +! crsdot(ivar) , name(ivar) , radius_influence(num_scan) , & +!BPR END dxd , u_banana , v_banana , pressure(kp) , passes , smooth_type , lat_center , & - use_first_guess ) +!BPR BEGIN + use_first_guess, fg_value, scale_cressman_rh_decreases ) +! use_first_guess ) +!BPR END + !BPR BEGIN + !IF (( name(ivar)(1:8) .EQ. 'RH ' ).AND.(kp.eq.9)) THEN + ! DO cury=1,jns_alloc + ! WRITE(130+num_scan,*) dum2d(:,cury) + ! ENDDO + !ENDIF + !BPR END + more_than_1_scan : IF ( max_scan .GT. 1 ) THEN IF ( use_first_guess ) THEN @@ -535,10 +803,24 @@ IF ( name(ivar)(1:8) .EQ. 'RH ' ) THEN CALL clean_rh ( dum2d , iew_alloc , jns_alloc , 0. , 100. ) END IF + + !BPR BEGIN + !If we are processing temperature, then save temperature for use in + !calculating mixing ratio once we start processing RH + IF ( name(ivar)(1:8) .EQ. 'TT ' ) THEN + t_first_guess = dum2d_fg + t_analysis = dum2d(1:iew_alloc-1,1:jns_alloc-1) + pressure_t = pressure(kp) + ENDIF + ! BPR END + ! Store the final analysis for kp level and for variable name(ivar). - CALL put_background_from_oa ( t , u , v , rh , slp_x , & + !BPR BEGIN + !CALL put_background_from_oa ( t , u , v , rh , slp_x , & + CALL put_background_from_oa ( t , u , v , rh , slp_x , pres , & + !BPR END dum2d , kp , name(ivar) , & iew_alloc , jns_alloc , kbu_alloc ) @@ -584,6 +866,9 @@ ELSE IF (kp .eq. kbu_alloc+1 ) print*," Doing Cressman for all upper level data" END IF + !BPR BEGIN + ELSEIF ( ( oa_type .EQ. 'None' ) ) THEN + !BPR END ELSE IF (kp .eq. kbu_alloc+1 ) print*," Doing Cressman for all upper level data" END IF Index: src/proc_oa.int =================================================================== --- src/proc_oa.int (revision 374) +++ src/proc_oa.int (working copy) @@ -1,9 +1,15 @@ SUBROUTINE proc_oa ( t , u , v , rh , slp_x , & -pressure , & +!BPR BEGIN +!pressure , & +pressure , pres, & +!BPR END iew_alloc , jns_alloc , kbu_alloc , & date , time , fdda_loop , mqd_count , mqd_abs_min , & total_numobs , num_obs_found , total_dups , & map_projection , obs , dxd , lat_center , & +!BPR BEGIN +use_p_tolerance_one_lev , & +!BPR END print_oa , print_found_obs , print_obs_files , use_first_guess , & smooth_type , smooth_sfc_wind , & smooth_sfc_temp , smooth_sfc_rh , & @@ -13,7 +19,11 @@ mqd_minimum_num_obs , mqd_maximum_num_obs , & oa_max_switch , radius_influence , & oa_min_switch , oa_3D_option , & -grid_id ) +!BPR BEGIN +!grid_id ) +grid_id, terrain, h, scale_cressman_rh_decreases, radius_influence_sfc_mult, & +oa_psfc, max_p_tolerance_one_lev_oa ) +!BPR END USE observation @@ -31,7 +41,18 @@ print_obs_files , & use_first_guess , & oa_min_switch , & - oa_max_switch +!BPR BEGIN +! oa_max_switch + oa_max_switch , & + scale_cressman_rh_decreases , & + oa_psfc +!BPR END + + !BPR BEGIN + REAL , INTENT ( IN ) :: radius_influence_sfc_mult + INTEGER , INTENT ( IN ) :: max_p_tolerance_one_lev_oa + LOGICAL , INTENT ( IN ) :: use_p_tolerance_one_lev + !BPR END INTEGER :: oa_3D_option INTEGER :: date , & time , & Index: src/proc_ob_access.F90 =================================================================== --- src/proc_ob_access.F90 (revision 374) +++ src/proc_ob_access.F90 (working copy) @@ -7,6 +7,9 @@ total_dups , map_projection , & get_value , get_x_location , get_y_location , get_longitude , & get_array_index , get_over_water , get_id , get_qc_info , & +!BPR BEGIN +get_elevation , get_fg_3d_t , get_fg_3d_h , & +!BPR END put_value , put_array_index , put_qc_info ) ! This routine accepts user level requests for modifying the observational @@ -55,6 +58,11 @@ get_qc_info LOGICAL , OPTIONAL , DIMENSION (:) :: get_over_water CHARACTER ( LEN = 8 ) , OPTIONAL , DIMENSION (:) :: get_id + !BPR BEGIN + REAL , OPTIONAL , DIMENSION (:) :: get_elevation + REAL , INTENT(IN) , OPTIONAL , DIMENSION(:,:,:) :: get_fg_3d_t, & + get_fg_3d_h + !BPR END ! Optional arguments for storing data. @@ -76,6 +84,12 @@ INCLUDE 'error.int' END INTERFACE + !BPR BEGIN + INTERFACE + INCLUDE 'query_ob.int' + END INTERFACE + !BPR END + ! There are a couple of simple errors to detect. Make sure that the surface ! is requested for sea level pressure. @@ -91,6 +105,37 @@ fatal , listing ) END IF oops_01 + ! BPR BEGIN + ! Make sure that the surface is requested for mean sea level pressure + ! derived from surface pressure. + + oops_01a : IF ( ( request_variable .EQ. 'PMSLPSFC' ) .AND. & + ( NINT ( request_level ) .NE. 1001 ) ) THEN + error_number = 99999001 + error_message(1:31) = 'proc_ob_access ' + error_message(32:) = ' For sea level pressure from surface pressure, request pressure & + &level 1001 mb.' + fatal = .true. + listing = .false. + CALL error_handler ( error_number , error_message , & + fatal , listing ) + END IF oops_01a + + ! Make sure that the surface is requested for surface pressure + + oops_01b : IF ( ( request_variable .EQ. 'PSFC ' ) .AND. & + ( NINT ( request_level ) .NE. 1001 ) ) THEN + error_number = 99999001 + error_message(1:31) = 'proc_ob_access ' + error_message(32:) = ' For surface pressure, request pressure & + &level 1001 mb.' + fatal = .true. + listing = .false. + CALL error_handler ( error_number , error_message , & + fatal , listing ) + END IF oops_01b + ! BPR END + ! For QC, we are either retrieving data from the data structure (get), or we are ! storing it there (put). For OA, we are requesting the data (use), but a bit ! differently from (get). This routine is separated into those three pieces @@ -145,7 +190,10 @@ CALL query_ob ( obs(loop_index) , date , time , & request_variable , request_level , request_qc_max , request_p_diff , & - value , qc ) +!BPR BEGIN +! value , qc ) + value , qc , fg_3d_h = get_fg_3d_h, fg_3d_t = get_fg_3d_t ) +!BPR END ! If qc is returned .NE. to missing, then this routine found a ! valid observation. @@ -166,6 +214,10 @@ IF ( PRESENT ( get_id ) ) & get_id(counter) = obs(loop_index)%location%id(1:8) + !BPR BEGIN + IF ( PRESENT ( get_elevation ) ) get_elevation(counter) = obs(loop_index)%info%elevation + !BPR END + ! Is this observation located in our domain? We can ! check easily if this has x,y available. Index: src/proc_ob_access.int =================================================================== --- src/proc_ob_access.int (revision 374) +++ src/proc_ob_access.int (working copy) @@ -5,6 +5,9 @@ total_dups , map_projection , & get_value , get_x_location , get_y_location , get_longitude , & get_array_index , get_over_water , get_id , get_qc_info , & +!BPR BEGIN +get_elevation , get_fg_3d_t , get_fg_3d_h , & +!BPR END put_value , put_array_index , put_qc_info ) USE observation @@ -32,6 +35,11 @@ get_qc_info LOGICAL , OPTIONAL , DIMENSION (:) :: get_over_water CHARACTER ( LEN = 8 ) , OPTIONAL , DIMENSION (:) :: get_id +!BPR BEGIN + REAL , OPTIONAL , DIMENSION (:) :: get_elevation + REAL , INTENT(IN) , OPTIONAL , DIMENSION(:,:,:) :: get_fg_3d_t, & + get_fg_3d_h +!BPR END REAL , OPTIONAL , DIMENSION (:) :: put_value INTEGER , OPTIONAL , DIMENSION (:) :: put_array_index , & put_qc_info Index: src/proc_obs_sort.F90 =================================================================== --- src/proc_obs_sort.F90 (revision 374) +++ src/proc_obs_sort.F90 (working copy) @@ -3,7 +3,10 @@ SUBROUTINE proc_obs_sort ( obs_filename , unit , & obs , number_of_obs , total_number_of_obs , fatal_if_exceed_max_obs , total_dups , index1 , & print_out_obs_found , print_out_files , & -levels , pressure , & +!BPR BEGIN +!levels , pressure , & +levels , pressure , slp_x , temperature , use_p_tolerance_one_lev , & +!BPR END max_p_extend_t , max_p_extend_w , & height , iew , jns , map_projection , date , time , fdda_loop ) @@ -28,6 +31,11 @@ REAL :: max_p_extend_t , & max_p_extend_w INTEGER :: iew , jns , map_projection + !BPR BEGIN + REAL , INTENT ( IN ) , DIMENSION ( jns, iew ) :: slp_x + REAL , INTENT ( IN ) , DIMENSION ( jns , iew , levels ) :: temperature + LOGICAL , INTENT ( IN ) :: use_p_tolerance_one_lev + !BPR END REAL , DIMENSION ( jns , iew , levels ) :: height INTEGER , INTENT(IN) :: date , time , fdda_loop CHARACTER (LEN=24) :: date_time_char @@ -48,7 +56,10 @@ CALL read_observations ( obs_filename , unit , obs , number_of_obs , & total_number_of_obs , fatal_if_exceed_max_obs , print_out_obs_found , & - height , pressure , iew , jns , levels , map_projection ) +!BPR BEGIN +! height , pressure , iew , jns , levels , map_projection ) + height , pressure , slp_x , temperature , iew , jns , levels , map_projection ) +!BPR END ! There should be at least a single observation to make the rest of ! this routine worth anyone's time. @@ -71,10 +82,26 @@ ! any missing levels that the analysis would like to have available ! for the QC and OA. The requested levels are the analysis levels. ! The only requirements are that there are at least 2 levels of data, - ! and that the observation has not already been discarded. While not - ! really an interpolation, the data that comes in at a specific level - ! (SATOB and AIREP), we should put them at the nearest analysis - ! pressure surface. + ! BPR BEGIN + !! and that the observation has not already been discarded. While not + !! really an interpolation, the data that comes in at a specific level + !! (SATOB and AIREP), we should put them at the nearest analysis + !! pressure surface. + ! and that the observation has not already been discarded. + ! For single-level above surface observations: + ! 1) If the user has not chosen to use a pressure tolerance for using + ! these obs (i.e., the user has set use_p_tolerance_one_lev = .FALSE.) + ! then FM-97 and FM-88 observations will be "adjusted" to the + ! neareast analysis pressure surface. Other single-level above-surface + ! observations will not be used unless they happen to fall exactly on an + ! analysis pressure level. + ! 2) If the user has chosen to use a pressure tolerance for using these + ! obs (i.e., the user has set use_p_tolerance_one_lev = .TRUE. ) + ! then all such obs will be treated the same. Namely, none will be + ! "adjusted" to the nearest analysis pressure surface. Rather, as long + ! as they are close enough to the nearest analysis pressure surface, they + ! will be directly compared against / used for that surface + ! BPR END DO i = 1 , number_of_obs IF ( ( ASSOCIATED ( obs(i)%surface ) ) .AND. & @@ -89,12 +116,23 @@ IF ( ASSOCIATED ( obs(i)%surface%next ) ) THEN CALL interp_all ( pressure , levels , obs(i)%surface , obs(i)%location%longitude ) END IF - - IF ( ( obs(i)%info%platform(1:11) .EQ. 'FM-97 AIREP' ) .OR. & - ( obs(i)%info%platform(1:11) .EQ. 'FM-88 SATOB' ) ) THEN - CALL extend_to_closest ( obs(i)%surface , pressure , levels , & - max_p_extend_t , max_p_extend_w ) + + + !BPR BEGIN + !If the user has chosen to use a pressure tolerance for comparing + !single-level above surface obs to the analysis pressure levels then + !we do not "adjust" the FM-97 and FM-88 data to the nearest analysis + !pressure surface + IF ( .NOT. use_p_tolerance_one_lev ) THEN + !BPR END + IF ( ( obs(i)%info%platform(1:11) .EQ. 'FM-97 AIREP' ) .OR. & + ( obs(i)%info%platform(1:11) .EQ. 'FM-88 SATOB' ) ) THEN + CALL extend_to_closest ( obs(i)%surface , pressure , levels , & + max_p_extend_t , max_p_extend_w ) + END IF + !BPR BEGIN END IF + !BPR END END IF END DO Index: src/proc_obs_sort.int =================================================================== --- src/proc_obs_sort.int (revision 374) +++ src/proc_obs_sort.int (working copy) @@ -1,7 +1,10 @@ SUBROUTINE proc_obs_sort ( obs_filename , unit , & obs , number_of_obs , total_number_of_obs , fatal_if_exceed_max_obs , total_dups , index , & print_out_obs_found , print_out_files , & -levels , pressure , & +!BPR BEGIN +!levels , pressure , & +levels , pressure , slp_x , temperature , use_p_tolerance_one_lev , & +!BPR END max_p_extend_t , max_p_extend_w , & height , iew , jns , map_projection , date , time , fdda_loop) @@ -18,6 +21,11 @@ print_out_files INTEGER , INTENT ( IN ) :: levels REAL , INTENT ( IN ) , DIMENSION ( levels ) :: pressure +!BPR BEGIN + REAL , INTENT ( IN ) , DIMENSION ( jns, iew ) :: slp_x + REAL , INTENT ( IN ) , DIMENSION ( jns , iew , levels ) :: temperature + LOGICAL , INTENT ( IN ) :: use_p_tolerance_one_lev +!BPR END REAL :: max_p_extend_t , & max_p_extend_w INTEGER :: iew , jns , map_projection Index: src/proc_qc.F90 =================================================================== --- src/proc_qc.F90 (revision 374) +++ src/proc_qc.F90 (working copy) @@ -3,15 +3,32 @@ total_dups , map_projection , & qc_test_error_max , qc_test_buddy , & qc_test_vert_consistency , qc_test_convective_adj , & +!BPR BEGIN + qc_psfc , & +!BPR END max_error_t , max_error_uv , & max_error_z , max_error_rh , & +!BPR BEGIN + max_error_dewpoint , & +!BPR END max_error_p , print_error_max , & max_buddy_t , max_buddy_uv , & max_buddy_z , max_buddy_rh , & - max_buddy_p , print_buddy , & +!BPR BEGIN + max_buddy_dewpoint , & +!BPR END +!BPR BEGIN +! max_buddy_p , print_buddy , & + max_buddy_p , & + use_p_tolerance_one_lev , max_p_tolerance_one_lev_qc , & + print_buddy , & +!BPR END print_found_obs , & print_vert , print_dry , & - pressure , date , time , dx , buddy_weight , & +!BPR BEGIN +! pressure , date , time , dx , buddy_weight , & + pressure , pres, date , time , dx , buddy_weight , & +!BPR END obs , index , max_number_of_obs , & t , u , v , h , rh , slp_x , sst , tobbox , odis ) @@ -58,17 +75,31 @@ LOGICAL , INTENT ( IN ) :: qc_test_error_max , & qc_test_buddy , & qc_test_vert_consistency , & - qc_test_convective_adj +!BPR BEGIN +! qc_test_convective_adj + qc_test_convective_adj , & + qc_psfc +!BPR END REAL , INTENT ( IN ) :: max_error_p, & max_error_t, & max_error_uv, & max_error_z, & max_error_rh, & +!BPR BEGIN + max_error_dewpoint, & +!BPR END max_buddy_p, & max_buddy_t, & max_buddy_uv, & max_buddy_z, & - max_buddy_rh +!BPR BEGIN +! max_buddy_rh + max_buddy_rh, & + max_buddy_dewpoint + + LOGICAL , INTENT ( IN ) :: use_p_tolerance_one_lev + INTEGER , INTENT ( IN ) :: max_p_tolerance_one_lev_qc +!BPR END LOGICAL , INTENT ( IN ) :: print_error_max, & print_buddy , & print_found_obs , & @@ -91,6 +122,9 @@ ! xob,yob: x, y locations of station obs ! lonob: longitude of the observation, for local time computation ! station_id: station identifier + !BPR BEGIN + ! station_elevation: Elevation of station in m MSL + !BPR END REAL , DIMENSION (max_number_of_obs) :: obs_value INTEGER , DIMENSION (max_number_of_obs) :: index_criteria @@ -99,6 +133,9 @@ CHARACTER ( LEN = 8 ) , DIMENSION (max_number_of_obs) :: station_id INTEGER , DIMENSION (max_number_of_obs) :: qc_flag LOGICAL , DIMENSION (max_number_of_obs) :: ship_report + !BPR BEGIN + REAL , DIMENSION (max_number_of_obs) :: station_elevation + !BPR END ! Data from the routine that supplies the background field for ! QC checking @@ -106,6 +143,10 @@ ! gridded: background/first field as input, final analysis as output REAL , DIMENSION (iew_alloc,jns_alloc) :: gridded + !BPR BEGIN + !Hold T/RH for calculating dewpoint + REAL , DIMENSION (iew_alloc,jns_alloc) :: gridded_t, gridded_rh + !BPR END ! Internally computed value. @@ -112,7 +153,11 @@ ! error_difference: initial maximum difference between observations and ! the interpolated value of the analysis at the ob ! location. - ! ivar: 1=U, 2=V, 3=T, 4=RH, 5=SLP + !BPR BEGIN + !! ivar: 1=U, 2=V, 3=T, 4=RH, 5=SLP + ! ivar: 1=U, 2=V, 3=T, 4=RH, 5=SLP, 6=PSFC, 7=DEWPOINT + !BPR END + ! kp: vertical index ! local_time: integer time (hhmm), crudely corrected to the local ! time based upon the longitude @@ -124,10 +169,32 @@ INTEGER , DIMENSION (max_number_of_obs) :: local_time CHARACTER ( LEN = 8 ) :: name + !BPR BEGIN + !Maximum difference between pressure of ob and pressure of background + !data it is compared against for QC (Pa) + !Default is 1 + INTEGER :: request_p_diff + !BPR END + !BPR BEGIN + INCLUDE 'constants.inc' + !BPR END + INTERFACE INCLUDE 'proc_ob_access.int' END INTERFACE + !BPR BEGIN + !request_p_diff = 700 + !If the user wants to allow a tolerance between an obs pressure and the + !first-guess field used to QC it, then use the user-specified tolerance. + !If not, then use a tolerance of 1 Pa, which is effectively no tolerance. + IF( use_p_tolerance_one_lev ) THEN + request_p_diff = max_p_tolerance_one_lev_qc + ELSE + request_p_diff = 1 + END IF + !BPR END + ! The first quality control (QC) that can be performed is with ! data that is vertically stacked. These reports have the temperature, ! speed and direction compared against reasonable benchmarks. @@ -175,12 +242,25 @@ vertical_level_loop : DO kp = 1, kbu_alloc - variable_loop : DO ivar = 1, 5 +!BPR BEGIN +! variable_loop : DO ivar = 1, 5 + variable_loop : DO ivar = 1, 7 ! do error checks only for ivar = PMSL and for kp = 1 + !IF ( ( kp .GT. 1 ) .AND. ( ivar .EQ. 5 ) ) EXIT variable_loop + + ! If ivar = PMSLPSFC or PMSL and this is not the lowest level (i.e., kp .ne. 1) + ! then go on to the next variable because PMSL and PMSL from PSFC + ! are only available at the lower level + IF ( ( kp .GT. 1 ) .AND. ( ( ivar .EQ. 5 ) .or. ( ivar .EQ. 6 ) ) ) CYCLE variable_loop + + !If user chose not to QC surface pressure then go to the next + !variable + IF ( ( .NOT. qc_psfc ) .AND. ( ivar .EQ. 6 ) )THEN + CYCLE variable_loop + END IF +!BPR END - IF ( ( kp .GT. 1 ) .AND. ( ivar .EQ. 5 ) ) EXIT variable_loop - which_variable : SELECT CASE ( ivar ) CASE ( 1 ) name = 'UU ' @@ -214,6 +294,30 @@ buddy_difference = max_buddy_p CALL yx2xy ( slp_x(1,1 ) , jns_alloc , iew_alloc , gridded ) gridded = gridded * 100 + !BPR BEGIN + CASE ( 6 ) + name = 'PMSLPSFC' + error_difference = max_error_p + buddy_difference = max_buddy_p + !CALL yx2xy ( pres(1,1,kp) , jns_alloc , iew_alloc , gridded ) + CALL yx2xy ( slp_x(1,1 ) , jns_alloc , iew_alloc , gridded ) + !Convert from hPa to Pa + gridded = gridded * 100 + CASE ( 7 ) + name = 'DEWPOINT' + error_difference = max_error_dewpoint + buddy_difference = max_buddy_dewpoint + CALL yx2xy ( t(1,1,kp) , jns_alloc , iew_alloc , gridded_t ) + CALL yx2xy ( rh(1,1,kp) , jns_alloc , iew_alloc , gridded_rh ) + !Calculate dewpoint using formula in obs_sort_module.F90 + !Since log(0) is undefined, if RH is near zero then set it to + !some small value + WHERE( gridded_rh .LT. very_small_rh ) + gridded = 1. / ( 1./gridded_t - Rv_over_L * LOG ( very_small_rh / 100. ) ) + ELSEWHERE + gridded = 1. / ( 1./gridded_t - Rv_over_L * LOG ( gridded_rh / 100. ) ) + END WHERE + !BPR END END SELECT which_variable ! Obtain observations for kp level and for variable ivar for this @@ -220,12 +324,19 @@ ! time period. CALL proc_ob_access ( 'get', name , print_found_obs , & - pressure(kp) , date , time , 1 , 2**20 , number_of_obs , num_obs_found , obs , & + !BPR BEGIN + !pressure(kp) , date , time , 1 , 2**20 , number_of_obs , num_obs_found , obs , & + pressure(kp) , date , time , request_p_diff , 2**20 , number_of_obs , num_obs_found , obs , & + !BPR END iew_alloc , jns_alloc , kbu_alloc , & total_dups , map_projection , & get_value=obs_value , get_x_location=xob , get_y_location=yob , & get_longitude=lonob , get_array_index=index_criteria , & - get_over_water=ship_report , get_id=station_id , get_qc_info=qc_flag ) + !BPR BEGIN + !get_over_water=ship_report , get_id=station_id , get_qc_info=qc_flag ) + get_over_water=ship_report , get_id=station_id , get_qc_info=qc_flag , & + get_elevation=station_elevation , get_fg_3d_t=t, get_fg_3d_h=h ) + !BPR END ! The local time will be used to modify the acceptable differences ! between the observations and the first guess analysis. @@ -238,12 +349,30 @@ IF ( kp .EQ. 1 ) THEN CALL ob_density ( xob , yob , dx , num_obs_found , tobbox , iew_alloc , jns_alloc ) END IF - + + + ! BPR BEGIN + IF( name .eq. 'PMSLPSFC') THEN + PRINT *,' ' + PRINT *,'****************************************************************************' + PRINT *,'Variable PMSLPSFC in the QC checks refers to the surface pressure in the ' + PRINT *,' obs (PSFC) being converted to sea level pressure (PMSL) for comparison to ' + PRINT *,' first-guess PMSL. This is done because the terrain assumed by the first ' + PRINT *,' guess may vary from the actual terrain and so first guess PSFC may not ' + PRINT *,' match actual PSFC. Using PMSL is less sensitive to this terrain mismatch.' + PRINT *,'****************************************************************************' + PRINT *,' ' + END IF + ! BPR END + ! Perform the maximum error difference QC test with the available ! observations and the background field analysis. IF ( qc_test_error_max ) THEN CALL error_max ( station_id , obs_value , ship_report , xob , yob , & + !BPR BEGIN + station_elevation , & + !BPR END index_criteria , qc_flag , num_obs_found , & gridded , iew_alloc , jns_alloc , & name , error_difference , pressure(kp) , local_time , print_error_max ) @@ -255,6 +384,9 @@ IF ( qc_test_buddy ) THEN CALL buddy_check ( station_id, obs_value , num_obs_found , xob, yob, qc_flag, & + !BPR BEGIN + station_elevation , & + !BPR END gridded, iew_alloc , jns_alloc , name, & pressure(kp), local_time, dx , buddy_weight , buddy_difference , print_buddy ) END IF @@ -263,7 +395,10 @@ ! qc_flags information to be stored. CALL proc_ob_access ( 'put', name , print_found_obs , & - pressure(kp) , date , time , 1 , 2**20 , number_of_obs , num_obs_found , obs , & + !BPR BEGIN + !pressure(kp) , date , time , 1 , 2**20 , number_of_obs , num_obs_found , obs , & + pressure(kp) , date , time , request_p_diff , 2**20 , number_of_obs , num_obs_found , obs , & + !BPR END iew_alloc , jns_alloc , kbu_alloc , & total_dups , map_projection , & put_value=obs_value , put_array_index=index_criteria , put_qc_info=qc_flag ) Index: src/proc_qc.int =================================================================== --- src/proc_qc.int (revision 374) +++ src/proc_qc.int (working copy) @@ -2,15 +2,32 @@ total_dups , map_projection , & qc_test_error_max , qc_test_buddy , & qc_test_vert_consistency , qc_test_convective_adj , & +!BPR BEGIN + qc_psfc , & +!BPR END max_error_t , max_error_uv , & max_error_z , max_error_rh , & +!BPR BEGIN + max_error_dewpoint , & +!BPR END max_error_p , print_error_max , & max_buddy_t , max_buddy_uv , & max_buddy_z , max_buddy_rh , & - max_buddy_p , print_buddy , & +!BPR BEGIN + max_buddy_dewpoint , & +!BPR END +!BPR BEGIN +! max_buddy_p , print_buddy , & + max_buddy_p , & + use_p_tolerance_one_lev , max_p_tolerance_one_lev_qc , & + print_buddy , & +!BPR END print_found_obs , & print_vert , print_dry , & - pressure , date , time , dx , buddy_weight , & +!BPR BEGIN +! pressure , date , time , dx , buddy_weight , & + pressure , pres , date , time , dx , buddy_weight , & +!BPR END obs , index , max_number_of_obs , & t , u , v , h , rh , slp_x , sst , tobbox , odis ) USE observation @@ -21,17 +38,30 @@ LOGICAL , INTENT ( IN ) :: qc_test_error_max , & qc_test_buddy , & qc_test_vert_consistency , & - qc_test_convective_adj +!BPR BEGIN +! qc_test_convective_adj + qc_test_convective_adj , & + qc_psfc +!BPR END REAL , INTENT ( IN ) :: max_error_p, & max_error_t, & max_error_uv, & max_error_z, & max_error_rh, & +!BPR BEGIN + max_error_dewpoint, & +!BPR END max_buddy_p, & max_buddy_t, & max_buddy_uv, & max_buddy_z, & - max_buddy_rh +!BPR BEGIN +! max_buddy_rh + max_buddy_rh, & + max_buddy_dewpoint + LOGICAL , INTENT ( IN ) :: use_p_tolerance_one_lev + INTEGER , INTENT ( IN ) :: max_p_tolerance_one_lev_qc +!BPR END LOGICAL , INTENT ( IN ) :: print_error_max, & print_buddy , & print_found_obs , & @@ -50,7 +80,11 @@ u , & v , & h , & - rh +!BPR BEGIN +! rh + rh , & + pres +!BPR END REAL , DIMENSION(jns,iew) :: slp_x , & sst , & Index: src/qc0_module.F90 =================================================================== --- src/qc0_module.F90 (revision 374) +++ src/qc0_module.F90 (working copy) @@ -33,8 +33,51 @@ END SUBROUTINE modify_qc_flag +!BPR BEGIN !------------------------------------------------------------------------------ +!Add qc_flag_to_add to qc_flag unless qc_flag already includes qc_flag_to_add +SUBROUTINE add_to_qc_flag ( qc_flag, qc_flag_to_add ) + + IMPLICIT NONE + + INTEGER , INTENT ( INOUT ) :: qc_flag + INTEGER , INTENT ( IN ) :: qc_flag_to_add + + INTEGER :: qc_small , & + qc_large + + !Since mod(a,p) = a - [ int(a/p) * p ], + ! qc_small = qc_flag - [ int(qc_flag/qc_flag_to_add) * qc_flag_to_add ] + !For simplicity flag = qc_flag and add = qc_flag_to_add, then + ! qc_small = flag - [ int(flag/add) * add ] + qc_small = mod ( qc_flag , qc_flag_to_add ) + !Since qc_flag, qc_flag_to_add, and "2" are all integers this will be integer + !math and so equivalent to: + !qc_large= 2 * int [ flag / ( 2 * add ) ] + qc_large = ( qc_flag / ( qc_flag_to_add * 2 ) ) * 2 + !The following equation expanded out and rearranged: + ! Term# 1 2 3 4 + ! [ flag ] [ flag ] + !flag_new = 2 * add * int [------------] - add * int [------] + flag + add + ! [ 2 * add ] [ add ] + !So the last two terms do the simple adding of "add" to the current flag + !"flag". The first two terms check to see if with integer math the "2" in the + !integer divide cancels the 2 outside the integer divide. This is checking + !to see if add is in flag already. + !If add is NOT already in flag then the 2's effectively cancel each other and + !thus the first two terms cancel each other out. + !If add is already in flag then 2's do not effectively cancel each other and + !the first term is the second term minus "add" so the sum of the first two + !terms in "-add", which cancels the fourth term and we are left with the + !original flag. + qc_flag = qc_large*qc_flag_to_add + qc_small + qc_flag_to_add + +END SUBROUTINE add_to_qc_flag +!BPR END + +!------------------------------------------------------------------------------ + SUBROUTINE ob_density ( xob , yob , grid_dist , numobs , tobbox , iew , jns ) ! Compute a guess at the observation density for the surface FDDA. @@ -87,7 +130,6 @@ ! no longer needed. ! ! Aijun Deng 06/11/2008 at Penn State -! IMPLICIT NONE INTEGER :: jns, iew ! Array x and y dimensions @@ -96,24 +138,78 @@ REAL, DIMENSION( : , : ), INTENT(OUT) :: distance ! Outgoing distance array INTEGER :: i, j, ii, jj - REAL :: dijmin, di, dj, dij + !BPR BEGIN + !REAL :: dijmin, di, dj, dij + REAL :: dijmin, dij + INTEGER :: di, dj + INTEGER :: ALLOCATE_STATUS, DEALLOCATE_STATUS + !BPR END + REAL :: domain_diagonal_m, domain_diagonal_grid_cells + INTEGER :: delta_i, delta_j + REAL, DIMENSION( jns+1, iew+1 ) :: dij_lookup + LOGICAL, DIMENSION( : , : ), ALLOCATABLE :: ob_in_box - distance(:,:) = sqrt( (iew-1)*dx*(iew-1)*dx + (jns-1)*dx*(jns-1)*dx ) ! Initialize - ! the distance array to be the domain diagonal size (something big). + !BPR BEGIN + !Modified subroutine to increase computational efficiency + + !Calculate the distance from non-adjacent corners of the domain + domain_diagonal_m = sqrt( (iew-1)*dx*(iew-1)*dx + (jns-1)*dx*(jns-1)*dx ) !meters + domain_diagonal_grid_cells = domain_diagonal_m / dx !Grid cells + + !distance(:,:) = sqrt( (iew-1)*dx*(iew-1)*dx + (jns-1)*dx*(jns-1)*dx ) ! Initialize + ! ! the distance array to be the domain diagonal size (something big). + !Initialize the distance array to be the domain diagonal size (i.e., !something big) + distance(:,:) = domain_diagonal_m + + !Create a lookup table with the distance (in grid cells) for two points that + !differ by a given number of points in the x/y direction + !Note that you pass as the first element: the difference in x plus one + !Note that you pass as the second element: the difference in y plus one + DO delta_i = 1, iew + 1 + DO delta_j = 1, jns + 1 + dij_lookup(delta_j,delta_i) = SQRT(REAL(delta_j-1)*REAL(delta_j-1)+REAL(delta_i-1)*REAL(delta_i-1)) + ENDDO + ENDDO + + !ob_in_box = + ALLOCATE(ob_in_box(lbound(tobbox,1):ubound(tobbox,1),lbound(tobbox,2):ubound(tobbox,2)), & + STAT = ALLOCATE_STATUS) + IF( ALLOCATE_STATUS .ne. 0 ) THEN + PRINT *,'ERROR: obs_distance: Could not allocate ob_in_box' + STOP 'ERROR: obs_distance: Could not allocate ob_in_box' + END IF + WHERE( tobbox > 0.5) + ob_in_box = .TRUE. + ELSEWHERE + ob_in_box = .FALSE. + ENDWHERE + dijmin = 0.0 DO i = 1, iew DO j = 1, jns - IF( tobbox(j,i) > 0.5 ) THEN + !BPR BEGIN + !IF( tobbox(j,i) > 0.5 ) THEN + IF( ob_in_box(j,i) ) THEN + !BPR END distance(j,i) = 0.0 ELSE - dijmin = MAXVAL( distance/dx ) + !BPR BEGIN + !dijmin = MAXVAL( distance/dx ) + dijmin = domain_diagonal_grid_cells + !BPR END DO ii = 1, iew DO jj = 1, jns - IF( tobbox(jj,ii) > 0.5 ) THEN - di = ABS( FLOAT(ii-i) ) - dj = ABS( FLOAT(jj-j) ) - dij = SQRT( di*di+dj*dj ) + !BPR BEGIN + !IF( tobbox(jj,ii) > 0.5 ) THEN + IF( ob_in_box(jj,ii) ) THEN + !di = ABS( FLOAT(ii-i) ) + !dj = ABS( FLOAT(jj-j) ) + !dij = SQRT( di*di+dj*dj ) + di = ABS( (ii-i) ) + dj = ABS( (jj-j) ) + dij = dij_lookup( dj+1,di+1 ) + !BPR END dijmin = MIN( dijmin,dij ) ENDIF ENDDO @@ -124,6 +220,14 @@ ENDIF ENDDO ENDDO + + !BPR BEGIN + DEALLOCATE(ob_in_box, STAT = DEALLOCATE_STATUS) + IF( DEALLOCATE_STATUS .ne. 0 ) THEN + PRINT *,'ERROR: obs_distance: Could not deallocate ob_in_box' + STOP 'ERROR: obs_distance: Could not deallocate ob_in_box' + END IF + !BPR END END SUBROUTINE obs_distance @@ -277,22 +381,48 @@ qc_t = current%meas%temperature%qc qc_rh = current%meas%rh%qc - IF ( ( qc_t .GE. fails_error_max ) .AND. ( qc_rh .LT. fails_error_max ) ) THEN - current%meas%rh%qc = current%meas%rh%qc + fails_error_max - current%meas%dew_point%qc = current%meas%dew_point%qc + fails_error_max - ELSE IF ( ( qc_t .GE. fails_buddy_check ) .AND. ( qc_rh .LT. fails_buddy_check ) ) THEN - current%meas%rh%qc = current%meas%rh%qc + fails_buddy_check - current%meas%dew_point%qc = current%meas%dew_point%qc + fails_buddy_check - END IF + + !BPR BEGIN + !The following code appears to contain errors + !For example, consider the case qc_t .eq. fails_buddy_check and qc_rh .lt. fails_error_max + !Since currently fails_buddy_check .gt. fails_error_max, this will lead to qc_rh = fails_error_max + !It would seem that we would want qc_rh = fails_buddy_check + !IF ( ( qc_t .GE. fails_error_max ) .AND. ( qc_rh .LT. fails_error_max ) ) THEN + ! current%meas%rh%qc = current%meas%rh%qc + fails_error_max + ! current%meas%dew_point%qc = current%meas%dew_point%qc + fails_error_max + !ELSE IF ( ( qc_t .GE. fails_buddy_check ) .AND. ( qc_rh .LT. fails_buddy_check ) ) THEN + ! current%meas%rh%qc = current%meas%rh%qc + fails_buddy_check + ! current%meas%dew_point%qc = current%meas%dew_point%qc + fails_buddy_check + !END IF + IF ( ( contains_2n ( qc_t, fails_error_max ) ) .AND. & + ( .NOT. contains_2n ( qc_rh, fails_error_max ) ) ) THEN + CALL add_to_qc_flag( current%meas%rh%qc, fails_error_max ) + CALL add_to_qc_flag( current%meas%dew_point%qc, fails_error_max ) + ENDIF + IF ( ( contains_2n ( qc_t, fails_buddy_check ) ) .AND. & + ( .NOT. contains_2n ( qc_rh, fails_buddy_check ) ) ) THEN + CALL add_to_qc_flag( current%meas%rh%qc, fails_buddy_check ) + CALL add_to_qc_flag( current%meas%dew_point%qc, fails_buddy_check ) + ENDIF + ! If td > t, then td and RH have to be bad. IF ( ( current%meas%dew_point%data .GT. current%meas%temperature%data ) .AND. & ( qc_rh .LT. fails_buddy_check ) ) THEN - current%meas%rh%qc = current%meas%rh%qc + fails_buddy_check - current%meas%dew_point%qc = current%meas%dew_point%qc + fails_buddy_check + !BPR BEGIN + !Since the conditional ensures that rh%qc does not contain the fails_buddy_check + !QC flag, for RH the following line is ok, but for dew_point we do not know + !if dewpoint already includes fails_buddy_check. For consistency, and to + !prevent potential future issues if code is modified, modify rh%qc + !even though it is not currently necessary + !current%meas%rh%qc = current%meas%rh%qc + fails_buddy_check + CALL add_to_qc_flag( current%meas%rh%qc , fails_buddy_check ) + !current%meas%dew_point%qc = current%meas%dew_point%qc + fails_buddy_check + CALL add_to_qc_flag( current%meas%dew_point%qc , fails_buddy_check ) + !BPR END END IF - + ! Point to the next level current => current%next Index: src/qc1_module.F90 =================================================================== --- src/qc1_module.F90 (revision 374) +++ src/qc1_module.F90 (working copy) @@ -7,6 +7,9 @@ !------------------------------------------------------------------------------ SUBROUTINE error_max ( station_id , obs , ship_report , xob , yob , & +!BPR BEGIN +station_elevation , & +!BPR END index , qc_flag , numobs , & gridded, iew , jns , & name , max_difference , pressure , local_time , print_error_max ) @@ -33,6 +36,9 @@ INTEGER :: numobs CHARACTER ( LEN = 8 ) , DIMENSION ( : ) :: station_id REAL , DIMENSION ( : ) :: obs, xob, yob + !BPR BEGIN + REAL , DIMENSION ( : ) :: station_elevation + !BPR END LOGICAL , DIMENSION ( : ) :: ship_report INTEGER , DIMENSION ( : ) :: index INTEGER , DIMENSION ( : ) :: qc_flag @@ -55,11 +61,17 @@ factored_difference , & aob REAL :: plevel - CHARACTER ( LEN = 100 ) :: message + !BPR BEGIN + !CHARACTER ( LEN = 100 ) :: message + CHARACTER ( LEN = 125 ) :: message + !BPR END REAL :: scale CHARACTER ( LEN = 8 ) , DIMENSION ( numobs ) :: station_id_small REAL , DIMENSION ( numobs ) :: obs_small , xob_small , yob_small + !BPR BEGIN + REAL , DIMENSION ( numobs ) :: station_elevation_small + !BPR END LOGICAL , DIMENSION ( numobs ) :: ship_report_small INTEGER , DIMENSION ( numobs ) :: index_small INTEGER , DIMENSION ( numobs ) :: qc_flag_small @@ -79,6 +91,9 @@ obs_small(:numobs) = obs(:numobs) xob_small(:numobs) = xob(:numobs) yob_small(:numobs) = yob(:numobs) + !BPR BEGIN + station_elevation_small(:numobs) = station_elevation(:numobs) + !BPR END ship_report_small(:numobs) = ship_report(:numobs) index_small(:numobs) = index(:numobs) qc_flag_small(:numobs) = qc_flag(:numobs) @@ -122,6 +137,36 @@ !!!!!!! WHERE ( shipreport_small ) factor = 0.75 +!BPR BEGIN + CASE ('PMSLPSFC') error_factor_by_field + ! Loop through each of the station locations (numobs). + station_loop_0a : DO num = 1, numobs + !Increase the maximum allowed error in sea level pressure derived from + !surface pressure as station elevation increases. Since assumptions must + !be made in this calculation above regarding the atmosphere below ground, + !the thicker the layer through which these assumptions are made the + !larger the potential error caused by the assumptions + factor(num)=1.0+(station_elevation(num)/2000.0) + END DO station_loop_0a + + CASE ( 'DEWPOINT' ) error_factor_by_field + ! Loop through each of the station locations (numobs). + ! Make the tolerance larger for lower dewpoints since the same dewpoint + ! difference represents less moisture difference at low dewpoints + station_loop_0b : DO num = 1, numobs + IF((obs_small(num).lt.255.0).and.(obs_small(num).ge.235.0)) THEN + factor(num)=1.25 + ELSEIF((obs_small(num).lt.235.0).and.(obs_small(num).ge.215.0)) THEN + factor(num)=1.50 + ELSEIF(obs_small(num).lt.215.0) THEN + factor(num)=1.75 + ELSE + factor(num)=1.00 + END IF + + END DO station_loop_0b + +!BPR END CASE ( 'RH ' , 'PMSL ' ) error_factor_by_field END SELECT error_factor_by_field @@ -182,7 +227,10 @@ factored_difference = max ( factor * max_difference , 0.15 * aob ) - CASE ( 'TT ' , 'RH ' , 'PMSL ' ) +!BPR BEGIN +! CASE ( 'TT ' , 'RH ' , 'PMSL ' ) + CASE ( 'TT ' , 'RH ' , 'PMSL ', 'PMSLPSFC', 'DEWPOINT' ) +!BPR END ! The maximum difference is scaled by the factor for the ! temperature, sea level pressure and the relative @@ -205,7 +253,10 @@ IF ( print_error_max ) THEN station_loop_2 : DO num = 1 , numobs - IF ( name(1:8) .EQ. 'PMSL ' ) THEN + !BPR BEGIN + !IF ( name(1:8) .EQ. 'PMSL ' ) THEN + IF ( ( name(1:8) .EQ. 'PMSL ' ) .OR. ( name(1:8) .EQ. 'PMSLPSFC' ) ) THEN + !BPR END scale = 100. ELSE scale = 1. @@ -216,11 +267,21 @@ & ",NAME=" ,a8, & & ",ERRMX=" ,f5.1,& & ",LOC=(" ,f6.1,",",f6.1,")",& + !BPR BEGIN + & ",PLEVEL=" ,f6.1,& + !BPR END & ",OBS=" ,f6.0,& & ",1ST GUESS=" ,f6.0,& & ",DIFF=" ,f5.0)' ) & station_id(num),name,factored_difference(num)/scale, & - xob(num),yob(num),obs(num)/scale,aob(num)/scale,err(num)/scale + !BPR BEGIN + !Also add pressure level being examined (plevel) -- note that this is not + !necessarily the pressure of the ob, but it is the pressure level for which we + !were doing QC when we found this ob so it must be at least close to this + !pressure + !xob(num),yob(num),obs(num)/scale,aob(num)/scale,err(num)/scale + xob(num),yob(num),plevel,obs(num)/scale,aob(num)/scale,err(num)/scale + !BPR END error_number = 00361001 error_message(1:31) = 'error_max ' error_message(32:) = message Index: src/qc2_module.F90 =================================================================== --- src/qc2_module.F90 (revision 374) +++ src/qc2_module.F90 (working copy) @@ -7,6 +7,9 @@ !------------------------------------------------------------------------------ SUBROUTINE buddy_check ( station_id, obs, numobs, xob, yob, qc_flag, & +!BPR BEGIN +station_elevation , & +!BPR END gridded, iew, jns, name, & pressure, local_time, dx , buddy_weight , buddy_difference , print_buddy ) @@ -37,6 +40,9 @@ CHARACTER ( LEN = 8 ) , DIMENSION ( : ) :: station_id REAL , DIMENSION ( : ) :: obs, xob, yob INTEGER , DIMENSION ( : ) :: qc_flag + !BPR BEGIN + REAL , DIMENSION ( : ) :: station_elevation + !BPR END REAL , DIMENSION( : , : ) :: gridded INTEGER :: iew, jns CHARACTER ( LEN = 8 ) :: name @@ -70,12 +76,20 @@ aob INTEGER , DIMENSION ( numobs ) :: buddy_num, & buddy_num_final - CHARACTER ( LEN = 100 ) :: message + !BPR BEGIN + INTEGER , DIMENSION ( numobs ) :: buddy_num_index + REAL , DIMENSION ( numobs ) :: factor + !CHARACTER ( LEN = 100 ) :: message + CHARACTER ( LEN = 125 ) :: message + !BPR END REAL :: scale CHARACTER ( LEN = 8 ) , DIMENSION ( numobs ) :: station_id_small REAL , DIMENSION ( numobs ) :: obs_small , xob_small , yob_small INTEGER , DIMENSION ( numobs ) :: qc_flag_small + !BPR BEGIN + REAL , DIMENSION ( numobs ) :: station_elevation_small + !BPR END INTEGER , DIMENSION ( numobs ) :: local_time_small INCLUDE 'error.inc' @@ -89,10 +103,19 @@ xob_small(:numobs) = xob(:numobs) yob_small(:numobs) = yob(:numobs) qc_flag_small(:numobs) = qc_flag(:numobs) + !BPR BEGIN + station_elevation_small(:numobs) = station_elevation(:numobs) + !BPR END local_time_small(:numobs) = local_time(:numobs) plevel = pressure + !BPR BEGIN + !Previously this value could be uninitialized when a conditional checks its + !value + buddy_num_final(:) = 0 + !BPR END + ! Define distance within which buddy check is to find neighboring stations: IF ( plevel .LE. 201.0 ) THEN @@ -118,6 +141,18 @@ difmax = buddy_difference END IF +!BPR BEGIN + CASE ( 'PMSLPSFC' ) + ! Loop through each of the station locations (numobs). + station_loop_0a : DO num = 1, numobs + !Increase the maximum allowed error in sea level pressure derived from + !surface pressure as station elevation increases. Since assumptions + !must be made in this calculation above the atmosphere below ground, + !the thicker the layer through which these assumptions are made the + !larger the potential error caused by the assumptions + difmax(num) = buddy_difference * ( 1.0+(station_elevation(num)/2000.0) ) + END DO station_loop_0a +!BPR END CASE ( 'PMSL ' ) difmax = buddy_difference @@ -131,8 +166,28 @@ END IF CASE ( 'RH ' ) + difmax = buddy_difference +!BPR BEGIN + CASE ( 'DEWPOINT' ) + ! Loop through each of the station locations (numobs). + ! Make the tolerance larger for lower dewpoints since the same dewpoint + ! difference represents less moisture difference at low dewpoints + station_loop_0b : DO num = 1, numobs + IF((obs_small(num).lt.255.0).and.(obs_small(num).ge.235.0)) THEN + factor(num)=1.25 + ELSEIF((obs_small(num).lt.235.0).and.(obs_small(num).ge.215.0)) THEN + factor(num)=1.50 + ELSEIF(obs_small(num).lt.215.0) THEN + factor(num)=1.75 + ELSE + factor(num)=1.00 + END IF + difmax(num) = buddy_difference * factor(num) + END DO station_loop_0b +!BPR END + END SELECT error_allowance_by_field difmax = difmax * buddy_weight @@ -181,6 +236,10 @@ IF ( distance .LE. range .AND. distance .GE. 0.1 ) THEN buddy_num(num) = buddy_num(num) + 1 diff(buddy_num(num)) = obs2(num3) - aob(num3) +!BPR BEGIN + !Keep track of which ob the current buddy is + buddy_num_index(buddy_num(num)) = num3 +!BPR END END IF END IF @@ -199,13 +258,25 @@ ! Check if there is any bad obs among the obs located within the ! the radius surrounding the test ob. - diff_check_1 = difmax (1) * 1.25 - diff_check_2 = difmax (1) + !BPR BEGIN + !This depended on difmax of the first ob, no matter which ob we are + !examining. This does not make sense now the difmax is actually ob + !dependent. + !diff_check_1 = difmax (1) * 1.25 + !diff_check_2 = difmax (1) + !BPR END check_bad_ob : IF ( buddy_num(num) .GE. 2 ) THEN kob = buddy_num(num) remove_bad_ob : DO numj = 1, buddy_num(num) + !BPR BEGIN + !Make the maximum allowed differences dependent on difmax for the + !current buddy rather than for whatever ob is the first in the + !overall ob array + diff_check_1 = difmax (buddy_num_index(numj)) * 1.25 + diff_check_2 = difmax (buddy_num_index(numj)) + !BPR END diff_numj = ABS ( diff ( numj ) - average ) IF ( diff ( numj ) .GT. diff_check_1 & .AND. diff_numj .GT. diff_check_2 ) THEN @@ -214,6 +285,13 @@ END IF END DO remove_bad_ob buddy_num_final(num) = kob + !BPR BEGIN + !Since buddies may have been removed buddy_num_index may no longer + !accurately reflect the indices of the buddies. Since we do not need + !this variable after this point, fill the array with clearly bad + !indices so this variable is not accidentally used. + buddy_num_index = -999 + !BPR END ! We may have removed too many observations. @@ -222,13 +300,27 @@ err(num) = diff_at_obs(num) - average ELSE err(num) = 0. - qc_flag_small(num) = qc_flag_small(num) + no_buddies + !BPR BEGIN + !If qc_flag_small(num) already includes the no_buddies flag this + !will add it again + !qc_flag_small(num) = qc_flag_small(num) + no_buddies + !Instead, the following call will only add the no_buddies flag if it is + !not already in qc_flag_small(num) + CALL add_to_qc_flag( qc_flag_small(num), no_buddies ) + !BPR END END IF ELSE check_bad_ob err(num) = 0. - qc_flag_small(num) = qc_flag_small(num) + no_buddies + !BPR BEGIN + !If qc_flag_small(num) already includes the no_buddies flag this + !will add it again + !qc_flag_small(num) = qc_flag_small(num) + no_buddies + !Instead, the following call will only add the no_buddies flag if it is + !not already in qc_flag_small(num) + CALL add_to_qc_flag( qc_flag_small(num), no_buddies ) + !BPR END END IF check_bad_ob @@ -249,7 +341,10 @@ IF ( print_buddy ) THEN station_loop_4 : DO num = 1 , numobs - IF ( name(1:8) .EQ. 'PMSL ' ) THEN +!BPR BEGIN +! IF ( name(1:8) .EQ. 'PMSL ' ) THEN + IF ( ( name(1:8) .EQ. 'PMSL ' ) .OR. ( name(1:8) .EQ. 'PMSLPSFC' ) ) THEN +!BPR END scale = 100. ELSE scale = 1. @@ -260,11 +355,21 @@ & ",NAME=" ,a8, & & ",BUDDY=" ,f5.1,& & ",LOC=(" ,f6.1,",",f6.1,")",& + !BPR BEGIN + & ",PLEVEL=" ,f6.1,& + !BPR END & ",OBS=" ,f6.0,& & ",1ST GUESS=" ,f6.0,& & ",DIFF=" ,f5.0)' ) & station_id(num),name,difmax(num)/scale, & - xob(num),yob(num),obs(num)/scale,aob(num)/scale,err(num) +!BPR BEGIN +!BUG FIX: Standard code fails to scale error in this diagnostic printout +!Also add pressure level being examined (plevel) -- note that this is not +!necessarily the pressure of the ob, but it is the pressure level for which we +!were doing QC when we found this ob so it must be at least close to this pressure +! xob(num),yob(num),obs(num)/scale,aob(num)/scale,err(num) + xob(num),yob(num),plevel,obs(num)/scale,aob(num)/scale,err(num)/scale +!BPR END error_number = 00364001 error_message(1:31) = 'buddy_check ' error_message(32:) = message Index: src/query_ob.int ===================================================================