diff --git a/Registry/registry.var b/Registry/registry.var index 5524d7446f..e63a3dcace 100644 --- a/Registry/registry.var +++ b/Registry/registry.var @@ -192,6 +192,11 @@ rconfig logical use_atmsobs namelist,wrfvar4 1 .false. - "use rconfig logical use_mwtsobs namelist,wrfvar4 1 .false. - "use_mwtsobs" "" "" rconfig logical use_mwhsobs namelist,wrfvar4 1 .false. - "use_mwhsobs" "" "" rconfig logical use_mwhs2obs namelist,wrfvar4 1 .false. - "use_mwhs2obs" "" "" +rconfig logical use_varbc_tamdar namelist,wrfvar4 1 .false. - "use_varbc_tamdar" "" "" +rconfig integer varbc_tamdar_bm namelist,wrfvar4 1 1 - "For VarBC: 1-ADR; 2-WMP" "" "" +rconfig integer varbc_tamdar_nbgerr namelist,wrfvar4 1 1000 - "Scaling B-beta" "" "" +rconfig integer varbc_tamdar_nobsmin namelist,wrfvar4 1 6 - "OB count req to run VarBC""" "" +rconfig real varbc_tamdar_pred0 namelist,wrfvar4 1 1.0 - "Predictor(1)" "" "" rconfig logical check_max_iv namelist,wrfvar5 1 .true. - "check_max_iv" "" "" rconfig real max_error_t namelist,wrfvar5 1 5.0 - "max_error_t" "" "" rconfig real max_error_uv namelist,wrfvar5 1 5.0 - "max_error_uv" "" "" diff --git a/var/build/da.make b/var/build/da.make index 07880209ee..70eebd2c71 100644 --- a/var/build/da.make +++ b/var/build/da.make @@ -22,6 +22,7 @@ WRFVAR_OBJS = \ da_sound.o \ da_mtgirs.o \ da_tamdar.o \ + da_varbc_tamdar.o \ da_bogus.o \ da_airep.o \ da_pilot.o \ diff --git a/var/build/depend.txt b/var/build/depend.txt index 438b62951e..5b9ee0f61e 100644 --- a/var/build/depend.txt +++ b/var/build/depend.txt @@ -125,7 +125,7 @@ da_interpolation.o : da_interpolation.f90 da_interp_msk_avg_2d_partial.inc da_in da_lapack.o : da_lapack.f90 dlamch.inc dlarf.inc dlarfg.inc dlarft.inc dlarfb.inc dorg2l.inc dorg2r.inc dsytd2.inc dlatrd.inc dorgqr.inc dorgql.inc dlassq.inc dlapy2.inc dlartg.inc dlasrt.inc dlansy.inc dsytrd.inc dsterf.inc dorgtr.inc dlae2.inc dlasr.inc dlaev2.inc dlascl.inc dlanst.inc dlaset.inc iparmq.inc ieeeck.inc ilaenv.inc dsteqr.inc dsyev.inc da_blas.o da_mat_cv3.o : da_mat_cv3.f90 da_metar.o : da_metar.f90 da_calculate_grady_metar.inc da_get_innov_vector_metar.inc da_check_max_iv_metar.inc da_transform_xtoy_metar_adj.inc da_transform_xtoy_metar.inc da_print_stats_metar.inc da_oi_stats_metar.inc da_residual_metar.inc da_jo_and_grady_metar.inc da_ao_stats_metar.inc da_tracing.o da_tools.o da_statistics.o da_physics.o da_grid_definitions.o da_par_util1.o da_par_util.o da_interpolation.o da_define_structures.o da_control.o module_domain.o -da_minimisation.o : da_minimisation.f90 da_read_basicstates.inc da_swap_xtraj.inc da_lanczos_io.inc da_kmat_mul.inc da_amat_mul.inc da_sensitivity.inc da_adjoint_sensitivity.inc da_transform_vtoy_adj.inc da_transform_vtoy.inc da_calculate_grady.inc da_minimise_lz.inc da_minimise_cg.inc da_write_diagnostics.inc da_dot_cv.inc da_dot.inc da_get_innov_vector.inc da_get_var_diagnostics.inc da_calculate_residual.inc da_jo_and_grady.inc da_calculate_gradj.inc da_calculate_j.inc da_transform_vtod_wpec.inc da_transform_vtod_wpec_adj.inc module_io_wrf.o da_4dvar.o da_lapack.o module_symbols_util.o da_wrf_interfaces.o da_vtox_transforms.o da_varbc.o da_transfer_model.o da_tracing.o da_tools_serial.o da_statistics.o da_synop.o da_ssmi.o da_sound.o da_ships.o da_satem.o da_reporting.o da_rain.o da_radar.o da_radiance1.o da_radiance.o da_tamdar.o da_mtgirs.o da_qscat.o da_pseudo.o da_profiler.o da_polaramv.o da_par_util1.o da_par_util.o da_pilot.o da_metar.o da_obs_io.o da_gpsref.o da_gpspw.o da_geoamv.o da_obs.o da_define_structures.o da_control.o da_buoy.o da_bogus.o da_airsr.o da_airep.o module_state_description.o module_domain.o module_dm.o module_configure.o da_join_iv_for_multi_inc.o da_wrfvar_io.o da_gpseph.o +da_minimisation.o : da_minimisation.f90 da_read_basicstates.inc da_swap_xtraj.inc da_lanczos_io.inc da_kmat_mul.inc da_amat_mul.inc da_sensitivity.inc da_adjoint_sensitivity.inc da_transform_vtoy_adj.inc da_transform_vtoy.inc da_calculate_grady.inc da_minimise_lz.inc da_minimise_cg.inc da_write_diagnostics.inc da_dot_cv.inc da_dot.inc da_get_innov_vector.inc da_get_var_diagnostics.inc da_calculate_residual.inc da_jo_and_grady.inc da_calculate_gradj.inc da_calculate_j.inc da_transform_vtod_wpec.inc da_transform_vtod_wpec_adj.inc module_io_wrf.o da_4dvar.o da_lapack.o module_symbols_util.o da_wrf_interfaces.o da_vtox_transforms.o da_varbc.o da_transfer_model.o da_tracing.o da_tools_serial.o da_statistics.o da_synop.o da_ssmi.o da_sound.o da_ships.o da_satem.o da_reporting.o da_rain.o da_radar.o da_radiance1.o da_radiance.o da_tamdar.o da_mtgirs.o da_qscat.o da_pseudo.o da_profiler.o da_polaramv.o da_par_util1.o da_par_util.o da_pilot.o da_metar.o da_obs_io.o da_gpsref.o da_gpspw.o da_geoamv.o da_obs.o da_define_structures.o da_control.o da_buoy.o da_bogus.o da_airsr.o da_airep.o module_state_description.o module_domain.o module_dm.o module_configure.o da_join_iv_for_multi_inc.o da_wrfvar_io.o da_gpseph.o da_varbc_tamdar.o da_module_convert_tool.o : da_module_convert_tool.f90 da_convertor_v_interp.inc da_module_couple_uv.o : da_module_couple_uv.f90 da_couple.inc da_calc_mu_uv.inc da_couple_uv.inc da_module_couple_uv_ad.o : da_module_couple_uv_ad.f90 da_couple_ad.inc da_calc_mu_uv_ad.inc da_couple_uv_ad.inc da_module_couple_uv.o @@ -160,7 +160,8 @@ da_spectral.o : da_spectral.f90 da_apply_power.inc da_legtra_inv_adj.inc da_vtov da_ssmi.o : da_ssmi.f90 da_sigma_v_tl.inc da_epsalt_tl.inc da_effang_tl.inc da_spemiss_tl.inc da_roughem_tl.inc da_effht_tl.inc da_tbatmos_tl.inc da_tb_tl.inc da_tbatmos_adj.inc da_spemiss_adj.inc da_roughem_adj.inc da_epsalt_adj.inc da_effht_adj.inc da_effang_adj.inc da_sigma_v_adj.inc da_tb_adj.inc da_calculate_grady_ssmt2.inc da_calculate_grady_ssmt1.inc da_calculate_grady_ssmi_rv.inc da_calculate_grady_ssmi_tb.inc da_transform_xtoy_ssmt2_adj.inc da_transform_xtoy_ssmt2.inc da_transform_xtoy_ssmt1_adj.inc da_transform_xtoy_ssmt1.inc da_print_stats_ssmt2.inc da_print_stats_ssmt1.inc da_oi_stats_ssmt2.inc da_oi_stats_ssmt1.inc da_ao_stats_ssmt2.inc da_ao_stats_ssmt1.inc da_get_innov_vector_ssmt2.inc da_get_innov_vector_ssmt1.inc da_get_innov_vector_ssmi_tb.inc da_get_innov_vector_ssmi_rv.inc da_check_max_iv_ssmt2.inc da_check_max_iv_ssmt1.inc da_check_max_iv_ssmi_tb.inc da_check_max_iv_ssmi_rv.inc da_residual_ssmt2.inc da_residual_ssmt1.inc da_jo_and_grady_ssmt2.inc da_jo_and_grady_ssmt1.inc da_transform_xtozrhoq_adj.inc da_transform_xtozrhoq_lin.inc da_transform_xtozrhoq.inc da_transform_xtoy_ssmi_tb_adj.inc da_transform_xtoy_ssmi_tb.inc da_transform_xtoy_ssmi_rv_adj.inc da_transform_xtoy_ssmi_rv.inc da_transform_xtotb_adj.inc da_transform_xtotb_lin.inc da_transform_xtotb.inc da_transform_xtoseasfcwind_adj.inc da_transform_xtoseasfcwind_lin.inc da_transform_xtoseasfcwind.inc da_transform_xtospeed_adj.inc da_transform_xtospeed_lin.inc da_transform_xtospeed.inc da_oi_stats_ssmi_tb.inc da_oi_stats_ssmi_rv.inc da_residual_ssmi_tb.inc da_residual_ssmi_rv.inc da_jo_and_grady_ssmi_tb.inc da_jo_and_grady_ssmi_rv.inc da_scan_obs_ssmi.inc da_read_obs_ssmi.inc da_ao_stats_ssmi_tb.inc da_ao_stats_ssmi_rv.inc da_tracing.o da_tools_serial.o da_tools.o da_statistics.o da_reporting.o da_par_util1.o da_par_util.o da_interpolation.o da_define_structures.o da_control.o module_ssmi.o module_domain.o da_statistics.o : da_statistics.f90 da_print_qcstat.inc da_stats_calculate.inc da_data_distribution.inc da_correlation_coeff2d.inc da_correlation_coeff1d.inc da_maxmin_in_field.inc da_analysis_stats.inc da_reporting.o da_tools_serial.o da_tracing.o da_par_util.o da_par_util1.o da_define_structures.o da_control.o module_domain.o da_synop.o : da_synop.f90 da_check_buddy_synop.inc da_calculate_grady_synop.inc da_check_max_iv_synop.inc da_get_innov_vector_synop.inc da_transform_xtoy_synop_adj.inc da_transform_xtoy_synop.inc da_print_stats_synop.inc da_oi_stats_synop.inc da_residual_synop.inc da_jo_synop_uvtq.inc da_jo_and_grady_synop.inc da_ao_stats_synop.inc da_tracing.o da_tools.o da_statistics.o da_physics.o da_grid_definitions.o da_par_util.o da_par_util1.o da_interpolation.o da_define_structures.o da_control.o module_domain.o -da_tamdar.o : da_tamdar.f90 da_calculate_grady_tamdar_sfc.inc da_check_max_iv_tamdar_sfc.inc da_get_innov_vector_tamdar_sfc.inc da_transform_xtoy_tamdar_sfc_adj.inc da_transform_xtoy_tamdar_sfc.inc da_print_stats_tamdar_sfc.inc da_oi_stats_tamdar_sfc.inc da_residual_tamdar_sfc.inc da_jo_tamdar_sfc_uvtq.inc da_jo_and_grady_tamdar_sfc.inc da_ao_stats_tamdar_sfc.inc da_calculate_grady_tamdar.inc da_get_innov_vector_tamdar.inc da_check_max_iv_tamdar.inc da_transform_xtoy_tamdar_adj.inc da_transform_xtoy_tamdar.inc da_print_stats_tamdar.inc da_oi_stats_tamdar.inc da_residual_tamdar.inc da_jo_tamdar_uvtq.inc da_jo_and_grady_tamdar.inc da_ao_stats_tamdar.inc da_tracing.o da_physics.o da_grid_definitions.o da_par_util1.o da_par_util.o da_tools.o da_statistics.o da_interpolation.o module_domain.o da_define_structures.o da_control.o +da_tamdar.o : da_tamdar.f90 da_calculate_grady_tamdar_sfc.inc da_check_max_iv_tamdar_sfc.inc da_get_innov_vector_tamdar_sfc.inc da_transform_xtoy_tamdar_sfc_adj.inc da_transform_xtoy_tamdar_sfc.inc da_print_stats_tamdar_sfc.inc da_oi_stats_tamdar_sfc.inc da_residual_tamdar_sfc.inc da_jo_tamdar_sfc_uvtq.inc da_jo_and_grady_tamdar_sfc.inc da_ao_stats_tamdar_sfc.inc da_calculate_grady_tamdar.inc da_get_innov_vector_tamdar.inc da_check_max_iv_tamdar.inc da_transform_xtoy_tamdar_adj.inc da_transform_xtoy_tamdar.inc da_print_stats_tamdar.inc da_oi_stats_tamdar.inc da_residual_tamdar.inc da_jo_tamdar_uvtq.inc da_jo_and_grady_tamdar.inc da_ao_stats_tamdar.inc da_tracing.o da_physics.o da_grid_definitions.o da_par_util1.o da_par_util.o da_tools.o da_statistics.o da_interpolation.o module_domain.o da_define_structures.o da_control.o da_varbc_tamdar.o +da_varbc_tamdar.o : da_varbc_tamdar.f90 da_varbc_tamdar_init.inc da_varbc_tamdar_pred.inc da_varbc_tamdar_precond.inc da_varbc_tamdar_direct.inc da_varbc_tamdar_adj.inc da_varbc_tamdar_tl.inc da_varbc_tamdar_update.inc da_tracing.o da_tools_serial.o da_tools.o da_reporting.o da_define_structures.o da_control.o module_dm.o da_test.o : da_test.f90 da_test_vxtransform.inc da_check_gradient.inc da_get_y_lhs_value.inc da_check_vtoy_adjoint.inc da_set_tst_trnsf_fld.inc da_check_psfc.inc da_check_sfc_assi.inc da_setup_testfield.inc da_check_xtoy_adjoint_buoy.inc da_check_xtoy_adjoint_profiler.inc da_check_xtoy_adjoint_ssmt2.inc da_check_xtoy_adjoint_ssmt1.inc da_check_xtoy_adjoint_qscat.inc da_check_xtoy_adjoint_pseudo.inc da_dot_cv.inc da_dot.inc da_check.inc da_check_gradient.inc da_transform_xtovp.inc da_check_xtoy_adjoint_rad.inc da_check_xtoy_adjoint_synop.inc da_check_xtoy_adjoint_tamdar_sfc.inc da_check_xtoy_adjoint_tamdar.inc da_check_xtoy_adjoint_mtgirs.inc da_check_xtoy_adjoint_sonde_sfc.inc da_check_xtoy_adjoint_sound.inc da_check_xtoy_adjoint_bogus.inc da_check_xtoy_adjoint_rain.inc da_check_xtoy_adjoint_radar.inc da_check_xtoy_adjoint_ships.inc da_check_xtoy_adjoint_polaramv.inc da_check_xtoy_adjoint_geoamv.inc da_check_xtoy_adjoint_satem.inc da_check_xtoy_adjoint_ssmi_tb.inc da_check_xtoy_adjoint_ssmi_rv.inc da_check_xtoy_adjoint_pilot.inc da_check_xtoy_adjoint_metar.inc da_check_xtoy_adjoint_gpsref.inc da_check_xtoy_adjoint_gpspw.inc da_check_xtoy_adjoint_airep.inc da_check_xtoy_adjoint.inc da_check_xtovptox_errors.inc da_check_vvtovp_adjoint.inc da_check_vp_errors.inc da_check_vptox_adjoint.inc da_check_vtox_adjoint.inc da_check_cvtovv_adjoint.inc da_check_balance.inc da_4dvar.o da_vtox_transforms.o da_wrfvar_io.o da_wrf_interfaces.o da_transfer_model.o da_tracing.o da_tools_serial.o da_statistics.o da_ssmi.o da_spectral.o da_reporting.o da_physics.o da_par_util1.o da_par_util.o da_obs.o da_minimisation.o da_ffts.o da_dynamics.o da_define_structures.o module_state_description.o module_domain.o da_control.o module_comm_dm.o module_dm.o module_configure.o da_rain.o da_check_dynamics_adjoint.inc da_check_xtoy_adjoint_gpseph.inc da_tools.o : da_tools.f90 da_geo2msl1.inc da_msl2geo1.inc da_get_time_slots.inc da_get_julian_time.inc da_get_print_lvl.inc da_get_3d_sum.inc da_get_2d_sum.inc da_set_boundary_3d.inc da_set_boundary_xb.inc da_set_boundary_xa.inc da_ludcmp.inc da_lubksb.inc da_eof_decomposition.inc da_eof_decomposition_test.inc da_buddy_qc.inc da_unifva.inc da_togrid.inc da_togrid_new.inc da_smooth_anl.inc da_openfile.inc da_gaus_noise.inc da_set_randomcv.inc da_random_omb.inc da_max_error_qc.inc da_add_noise_new.inc da_add_noise.inc da_residual_new.inc da_residual.inc da_diff_seconds.inc da_mo_correction.inc da_intpsfc_tem.inc da_intpsfc_prs.inc da_sfcprs.inc da_obs_sfc_correction.inc da_1d_eigendecomposition.inc da_convert_zk.inc da_lc_cone.inc da_set_merc.inc da_map_set.inc da_map_init.inc da_set_ps.inc da_set_lc.inc da_xyll_ps.inc da_xyll_merc.inc da_xyll_lc.inc da_xyll_latlon.inc da_xyll_default.inc da_xyll.inc da_llxy_wrf_new.inc da_llxy_wrf.inc da_llxy_ps_new.inc da_llxy_ps.inc da_llxy_merc_new.inc da_llxy_merc.inc da_llxy_lc_new.inc da_llxy_lc.inc da_llxy_latlon_new.inc da_llxy_latlon.inc da_llxy_rotated_latlon.inc da_llxy_global_new.inc da_llxy_global.inc da_llxy_kma_global_new.inc da_llxy_kma_global.inc da_llxy_default_new.inc da_llxy_default.inc da_llxy_new.inc da_llxy.inc da_map_utils_defines.inc da_lapack.o da_reporting.o da_tracing.o da_tools_serial.o da_define_structures.o da_control.o module_domain.o module_dm.o module_bc.o da_sfc_hori_interp_weights.inc da_tools_serial.o : da_tools_serial.f90 da_find_fft_trig_funcs.inc da_find_fft_factors.inc da_advance_time.inc da_advance_cymdh.inc da_array_print.inc da_change_date.inc da_free_unit.inc da_get_unit.inc da_reporting.o da_control.o @@ -191,7 +192,7 @@ da_wrfvar_esmf.o : da_wrfvar_esmf.f90 da_wrfvar_esmf_super.o : da_wrfvar_esmf_super.f90 da_wrfvar_interface.inc da_esmf_finalize.inc da_esmf_run.inc da_esmf_init.inc da_wrfvar_io.o : copyfile.c da_wrfvar_io.f90 da_med_initialdata_output_lbc.inc da_med_initialdata_output.inc da_med_initialdata_input.inc da_update_firstguess.inc da_4dvar.o da_tracing.o da_reporting.o da_control.o module_io_domain.o module_domain.o module_date_time.o module_configure.o da_wrfvar_main.o : da_wrfvar_main.f90 da_4dvar.o da_wrfvar_top.o da_wrf_interfaces.o da_tracing.o da_control.o module_symbols_util.o -da_wrfvar_top.o : da_wrfvar_top.f90 da_solve_init.inc da_solve_dual_res_init.inc da_solve.inc da_wrfvar_finalize.inc da_wrfvar_interface.inc da_wrfvar_run.inc da_wrfvar_init2.inc da_wrfvar_init1.inc da_wrf_interfaces.o da_rain.o da_synop.o da_ssmi.o da_sound.o da_ships.o da_satem.o da_radar.o da_mtgirs.o da_qscat.o da_profiler.o da_polaramv.o da_pilot.o da_metar.o da_gpsref.o da_gpspw.o da_geoamv.o da_buoy.o da_bogus.o da_airsr.o da_airep.o da_crtm.o da_tools.o da_vtox_transforms.o da_transfer_model.o da_tracing.o da_tools_serial.o da_test.o da_setup_structures.o da_reporting.o da_varbc.o da_radiance1.o da_physics.o da_par_util.o da_obs_io.o da_obs.o da_minimisation.o da_define_structures.o da_control.o module_comm_dm.o module_dm.o module_tiles.o module_state_description.o module_radiance.o da_wrfvar_io.o da_4dvar.o module_symbols_util.o module_driver_constants.o module_domain.o module_configure.o module_io_domain.o da_netcdf_interface.o da_gpseph.o +da_wrfvar_top.o : da_wrfvar_top.f90 da_solve_init.inc da_solve_dual_res_init.inc da_solve.inc da_wrfvar_finalize.inc da_wrfvar_interface.inc da_wrfvar_run.inc da_wrfvar_init2.inc da_wrfvar_init1.inc da_wrf_interfaces.o da_rain.o da_synop.o da_ssmi.o da_sound.o da_ships.o da_satem.o da_radar.o da_mtgirs.o da_qscat.o da_profiler.o da_polaramv.o da_pilot.o da_metar.o da_gpsref.o da_gpspw.o da_geoamv.o da_buoy.o da_bogus.o da_airsr.o da_airep.o da_crtm.o da_tools.o da_vtox_transforms.o da_transfer_model.o da_tracing.o da_tools_serial.o da_test.o da_setup_structures.o da_reporting.o da_varbc.o da_radiance1.o da_physics.o da_par_util.o da_obs_io.o da_obs.o da_minimisation.o da_define_structures.o da_control.o module_comm_dm.o module_dm.o module_tiles.o module_state_description.o module_radiance.o da_wrfvar_io.o da_4dvar.o module_symbols_util.o module_driver_constants.o module_domain.o module_configure.o module_io_domain.o da_netcdf_interface.o da_gpseph.o da_varbc_tamdar.o decode_airs.o : decode_airs.f90 module_read_airs.o f_qv_from_rh.o : f_qv_from_rh.f90 gamma1.o : gamma1.f90 da_control.o diff --git a/var/da/da_control/da_control.f90 b/var/da/da_control/da_control.f90 index 6f31cb9d0b..9d9c1a679d 100644 --- a/var/da/da_control/da_control.f90 +++ b/var/da/da_control/da_control.f90 @@ -278,6 +278,7 @@ module da_control integer :: check_max_iv_unit, check_buddy_unit, rand_unit, omb_unit, & filtered_obs_unit integer :: biasprep_unit, qcstat_conv_unit + integer :: varbc_tamdar_unit integer,parameter :: filename_len = 200 @@ -369,6 +370,7 @@ module da_control integer :: cv_size_domain_jp ! Total jp cv size. integer :: cv_size_domain_js ! Total js cv size. integer :: cv_size_domain_jl ! Total jl cv size. + integer :: cv_size_domain_jt ! Total jt cv size. integer :: cv_size_domain ! Total cv size. ! Hybrid: @@ -419,6 +421,8 @@ module da_control ! missing_r outside_of_domain = -77, & ! Data outside horizontal domain ! or time window, data set to missing_r + fail_varbc_aircraft = -55, & ! Data fail VarBC of aircraft + ! => no action wrong_direction = -15, & ! Wind vector direction <0 or> 360 ! => direction set to missing_r negative_spd = -14, & ! Wind vector norm is negative diff --git a/var/da/da_define_structures/da_define_structures.f90 b/var/da/da_define_structures/da_define_structures.f90 index c67b0f7d06..661221f757 100644 --- a/var/da/da_define_structures/da_define_structures.f90 +++ b/var/da/da_define_structures/da_define_structures.f90 @@ -405,6 +405,25 @@ module da_define_structures type (field_type) , pointer :: q (:) ! q. end type tamdar_type + type varbc_tamdar_type + character(len=40) :: fmt_param ! Format of parameter table + integer :: nmaxpred ! Max. No. of predictors + integer :: nphase ! No. of flight phases + integer :: nair ! No. of aircrafts in table + integer :: npred ! No. of predictors + integer :: nmaxobs ! Max Obs No. + integer , pointer :: nobs (:,:)! Obs No. in proc + integer , pointer :: nobs_sum(:,:)! Total Obs No. + integer , pointer :: tail_id (:)! Tail ID of aircrafts + integer , pointer :: obs_sn(:,:,:)! Serial No. of Obs in proc + integer , pointer :: ifuse (:,:)! run varbc or not + integer , pointer :: index (:,:,:)! Index in CV + real , pointer :: pred (:,:,:)! Predictors + real , pointer :: param (:,:,:)! Parameters + real , pointer :: bgerr (:,:,:)! Bkg err in Hessian + real , pointer :: vtox(:,:,:,:)! Transformation of CV + end type varbc_tamdar_type + type airsr_type real , pointer :: h (:) ! Height in m real , pointer :: p (:) ! pressure. @@ -665,6 +684,8 @@ module da_define_structures type (synop_type) , pointer :: tamdar_sfc(:) type (rain_type) , pointer :: rain(:) + type (varbc_tamdar_type) :: varbc_tamdar + real :: missing real :: ptop end type iv_type @@ -953,6 +974,7 @@ module da_define_structures real :: jl real :: jd real :: jm + real :: jt type (jo_type) :: jo end type j_type @@ -963,6 +985,7 @@ module da_define_structures integer :: size_jp ! Size of CV array for Jp term. integer :: size_js ! Size of CV array for Js term. integer :: size_jl ! Size of CV array for Jl term. + integer :: size_jt ! Size of CV array for Jt term. integer :: size1c ! Complex size of CV array of 1st variable error. integer :: size2c ! Complex size of CV array of 2nd variable error. integer :: size3c ! Complex size of CV array of 3rd variable error. diff --git a/var/da/da_main/da_solve.inc b/var/da/da_main/da_solve.inc index 84479e75f1..2b3526639e 100644 --- a/var/da/da_main/da_solve.inc +++ b/var/da/da_main/da_solve.inc @@ -480,6 +480,7 @@ be % cv % size_jp = 0 be % cv % size_js = 0 be % cv % size_jl = 0 + be % cv % size_jt = 0 !--------------------------------------------------------------------------- ! [5.1] Set up background errors (be): @@ -515,6 +516,13 @@ if (use_rad .and. (use_varbc.or.freeze_varbc)) call da_varbc_init(iv, be) #endif + if (use_tamdarobs .and. use_varbc_tamdar) then + call da_varbc_tamdar_init(iv) + else + use_varbc_tamdar = .false. + end if + if (use_varbc_tamdar) call da_varbc_tamdar_pred(iv, be, ob) + !--------------------------------------------------------------------------- ! [5.3] Set up satellite control variable: !--------------------------------------------------------------------------- @@ -525,7 +533,7 @@ !--------------------------------------------------------------------------- ! [5.4] Total control variable: !--------------------------------------------------------------------------- - be % cv % size = be%cv%size_jb + be%cv%size_je + be%cv%size_jp + be%cv%size_js + be%cv%size_jl + be % cv % size = be%cv%size_jb + be%cv%size_je + be%cv%size_jp + be%cv%size_js + be%cv%size_jl + be%cv%size_jt cv_size = be % cv % size !--------------------------------------------------------------------------- @@ -771,7 +779,7 @@ be,iv,xhat,qhat,cvt,re,y,j,eignvec,eignval,neign) ! Output Cost Function call da_calculate_j(it, 1, cv_size, be%cv%size_jb, be%cv%size_je, be%cv%size_jp, & - be%cv%size_jl, xbx, be, iv, xhat, cvt, re, y, j, grid, config_flags ) + be%cv%size_jl, be%cv%size_jt, xbx, be, iv, xhat, cvt, re, y, j, grid, config_flags ) else call da_minimise_lz(grid, config_flags, it, cv_size, xbx,& @@ -887,6 +895,8 @@ call da_varbc_update(it, cv_size, xhat, iv) #endif + if (use_varbc_tamdar) & + call da_varbc_tamdar_update(cv_size, xhat, iv) !------------------------------------------------------------------------ ! [8.10] Output WRFVAR analysis and analysis increments: @@ -1005,6 +1015,19 @@ if ( anal_type_hybrid_dual_res ) deallocate(aens_locs) + if (use_varbc_tamdar) then + deallocate (iv%varbc_tamdar%nobs) + deallocate (iv%varbc_tamdar%nobs_sum) + deallocate (iv%varbc_tamdar%tail_id) + deallocate (iv%varbc_tamdar%obs_sn) + deallocate (iv%varbc_tamdar%ifuse) + deallocate (iv%varbc_tamdar%index) + deallocate (iv%varbc_tamdar%pred) + deallocate (iv%varbc_tamdar%param) + deallocate (iv%varbc_tamdar%bgerr) + deallocate (iv%varbc_tamdar%vtox) + end if + deallocate ( c1f ) deallocate ( c2f ) deallocate ( c3f ) diff --git a/var/da/da_main/da_wrfvar_finalize.inc b/var/da/da_main/da_wrfvar_finalize.inc index 7a00ee4bf0..dc68c0595b 100644 --- a/var/da/da_main/da_wrfvar_finalize.inc +++ b/var/da/da_main/da_wrfvar_finalize.inc @@ -81,6 +81,10 @@ subroutine da_wrfvar_finalize call da_free_unit (jo_unit) call da_free_unit (check_max_iv_unit) call da_free_unit (check_buddy_unit ) + if (use_varbc_tamdar) then + close (varbc_tamdar_unit) + call da_free_unit (varbc_tamdar_unit) + end if end if #ifdef CRTM diff --git a/var/da/da_main/da_wrfvar_init2.inc b/var/da/da_main/da_wrfvar_init2.inc index 4ae8231725..cba1de11ec 100644 --- a/var/da/da_main/da_wrfvar_init2.inc +++ b/var/da/da_main/da_wrfvar_init2.inc @@ -275,6 +275,10 @@ subroutine da_wrfvar_init2 open(unit=jo_unit,file="jo",status="replace") open(unit=check_max_iv_unit,file="check_max_iv",status="replace") open(unit=check_buddy_unit ,file="buddy_check" ,status="replace") + if (use_varbc_tamdar) then + call da_get_unit (varbc_tamdar_unit) + open(unit=varbc_tamdar_unit,file="varbc_tamdar_detail.log",status="replace") + end if end if if (trace_use) call da_trace_exit("da_wrfvar_init2") diff --git a/var/da/da_main/da_wrfvar_top.f90 b/var/da/da_main/da_wrfvar_top.f90 index 25bc9618a9..b49dfbad13 100644 --- a/var/da/da_main/da_wrfvar_top.f90 +++ b/var/da/da_main/da_wrfvar_top.f90 @@ -110,6 +110,9 @@ module da_wrfvar_top use da_rain, only : da_oi_stats_rain use da_gpseph, only : da_gpseph_final + use da_varbc_tamdar, only : da_varbc_tamdar_init, da_varbc_tamdar_pred, & + da_varbc_tamdar_update + use da_wrf_interfaces use da_netcdf_interface, only : da_get_var_2d_real_cdf diff --git a/var/da/da_minimisation/da_adjoint_sensitivity.inc b/var/da/da_minimisation/da_adjoint_sensitivity.inc index ea239539fb..560eb26502 100644 --- a/var/da/da_minimisation/da_adjoint_sensitivity.inc +++ b/var/da/da_minimisation/da_adjoint_sensitivity.inc @@ -67,8 +67,8 @@ subroutine da_adjoint_sensitivity(grid, config_flags, cv_size, xbx, be, iv, xhat ! Jo for all observations !------------------------ call da_allocate_y (iv, re) - call da_calculate_j(1,1,cv_size,0,0,0,0,xbx,be,iv,xhat,cv,re,y,j_cost,grid,config_flags) - call da_calculate_gradj(1,1,cv_size,0,0,0,0,xbx,be,iv,xhat+cv,y,shat,grid,config_flags,re) + call da_calculate_j(1,1,cv_size,0,0,0,0,0,xbx,be,iv,xhat,cv,re,y,j_cost,grid,config_flags) + call da_calculate_gradj(1,1,cv_size,0,0,0,0,0,xbx,be,iv,xhat+cv,y,shat,grid,config_flags,re) call da_deallocate_y(re) case default; diff --git a/var/da/da_minimisation/da_calculate_gradj.inc b/var/da/da_minimisation/da_calculate_gradj.inc index 401d1d1917..ac77330bd8 100644 --- a/var/da/da_minimisation/da_calculate_gradj.inc +++ b/var/da/da_minimisation/da_calculate_gradj.inc @@ -1,5 +1,5 @@ subroutine da_calculate_gradj(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size_jp, & - cv_size_jl, xbx, be, iv, cv, y, grad, grid, config_flags, re ) + cv_size_jl, cv_size_jt, xbx, be, iv, cv, y, grad, grid, config_flags, re ) !--------------------------------------------------------------------------- ! Purpose: Calculates the gradient of the cost function w/r to cv @@ -15,7 +15,7 @@ subroutine da_calculate_gradj(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size integer, intent(in) :: it ! external iteration #. integer, intent(in) :: iter ! internal iteration #. integer, intent(in) :: cv_size ! Total cv size. - integer, intent(in) :: cv_size_jb, cv_size_je, cv_size_jp, cv_size_jl + integer, intent(in) :: cv_size_jb, cv_size_je, cv_size_jp, cv_size_jl, cv_size_jt type (xbx_type),intent(inout) :: xbx ! For header & non-grid arrays. type (be_type), intent(in) :: be ! background error structure. type (iv_type), intent(inout) :: iv ! innovation vector (o-b). @@ -29,6 +29,7 @@ subroutine da_calculate_gradj(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size integer :: je_start, je_end ! Start/end indices of Je. integer :: jl_start, jl_end ! Start/end indices of Jl. + integer :: jt_start, jt_end, iphase real :: jo_partial ! jo for this processor type (y_type) :: jo_grad_y ! Grad_y(jo) real :: grad_jo(cv_size) @@ -39,7 +40,8 @@ subroutine da_calculate_gradj(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size real :: grad_js(cv_size) real :: grad_jl(cv_size) real :: grad_jm(cv_size) - real :: gnorm_j, gnorm_jo, gnorm_jb, gnorm_je, gnorm_jd, gnorm_jp, gnorm_js, gnorm_jl + real :: grad_jt(cv_size) + real :: gnorm_j, gnorm_jo, gnorm_jb, gnorm_je, gnorm_jd, gnorm_jp, gnorm_js, gnorm_jl, gnorm_jt real :: gnorm_jm logical :: jcdf_flag @@ -62,7 +64,9 @@ subroutine da_calculate_gradj(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size jp_end = cv_size_jb + cv_size_je + cv_size_jp jl_start = cv_size_jb + cv_size_je + cv_size_jp + 1 jl_end = cv_size_jb + cv_size_je + cv_size_jp + cv_size_jl - + jt_start = jl_end + 1 + jt_end = jl_end + cv_size_jt + grad_jo = 0.0 grad_jb = 0.0 grad_je = 0.0 @@ -71,6 +75,7 @@ subroutine da_calculate_gradj(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size grad_js = 0.0 grad_jl = 0.0 grad_jm = 0.0 + grad_jt = 0.0 inc_div = 0.0 @@ -215,17 +220,37 @@ subroutine da_calculate_gradj(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size call da_transform_vtox_adj(grid, cv_size, xbx, be, grid%ep, grid%vp, grid%vv, grad_jm) end if + !------------------------------------------------------------------------- + ! [n.n] calculate grad_t (jt): + !------------------------------------------------------------------------- + + if (cv_size_jt > 0) then + do n = 1, iv%varbc_tamdar%nair + do iphase = 1, iv%varbc_tamdar%nphase + if (iv%varbc_tamdar%nobs(iphase,n) > 0 .and. iv%varbc_tamdar%ifuse(iphase,n) > 0) then + do ipred = 1, iv%varbc_tamdar%npred + id = iv%varbc_tamdar%index(ipred,iphase,n) + bgerr = iv%varbc_tamdar%bgerr(ipred,iphase,n) + if (bgerr > 0.0) & + grad_jt(id) = (1/sqrt(bgerr)) * SUM(cv(id) * & + iv%varbc_tamdar%vtox(ipred,1:iv%varbc_tamdar%npred,iphase,n)) + end do + end if + end do + end do + end if + !-------------------------------------------------------------------------------------------------- ! [7.0] calculate grad_v (j) = grad_v (jb) + grad_v (jo) + grad_v (je) + grad_v (jd) + grad_v (jp) + grad_v (js) + grad_v (jl) !-------------------------------------------------------------------------------------------------- - grad = grad_jo + grad_jb + grad_je + grad_jd + grad_jp + grad_js + grad_jl + grad_jm + grad = grad_jo + grad_jb + grad_je + grad_jd + grad_jp + grad_js + grad_jl + grad_jm + grad_jt !------------------------------------------------------------------------- ! [8.0] write Gradient: !------------------------------------------------------------------------- if (rootproc) then if (it == 1 .and. iter == 0) then - write(unit=grad_unit,fmt='(a)')'Outer EPS Inner G Gb Go Ge Gd Gp Gs Gl Gm' + write(unit=grad_unit,fmt='(a)')'Outer EPS Inner G Gb Go Ge Gd Gp Gs Gl Gm Gt' write(unit=grad_unit,fmt='(a)')'Iter Iter ' end if end if @@ -240,14 +265,13 @@ subroutine da_calculate_gradj(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size gnorm_js = sqrt(da_dot_cv(cv_size, grad_js, grad_js, grid, be%cv_mz, be%ncv_mz)) gnorm_jl = sqrt(da_dot_cv(cv_size, grad_jl, grad_jl, grid, be%cv_mz, be%ncv_mz)) gnorm_jm = sqrt(da_dot_cv(cv_size, grad_jm, grad_jm, grid, be%cv_mz, be%ncv_mz)) + gnorm_jt = sqrt(da_dot_cv(cv_size, grad_jt, grad_jt, grid, be%cv_mz, be%ncv_mz)) if (rootproc) & - write(grad_unit,fmt='(2x,i2,1x,e10.3,2x,i4,9(1x,f10.3))') & - it, eps(it), iter, gnorm_j, gnorm_jb, gnorm_jo, gnorm_je, gnorm_jd, gnorm_jp, gnorm_js, gnorm_jl, gnorm_jm + write(grad_unit,fmt='(2x,i2,1x,e10.3,2x,i4,10(1x,f10.3))') & + it, eps(it), iter, gnorm_j, gnorm_jb, gnorm_jo, gnorm_je, gnorm_jd, gnorm_jp, gnorm_js, gnorm_jl, gnorm_jm, gnorm_jt end if if (trace_use) call da_trace_exit("da_calculate_gradj") end subroutine da_calculate_gradj - - diff --git a/var/da/da_minimisation/da_calculate_j.inc b/var/da/da_minimisation/da_calculate_j.inc index 4f23bffc77..054d2df8f6 100644 --- a/var/da/da_minimisation/da_calculate_j.inc +++ b/var/da/da_minimisation/da_calculate_j.inc @@ -1,5 +1,5 @@ subroutine da_calculate_j(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size_jp, & - cv_size_jl, xbx, be, iv, xhat, cv, & + cv_size_jl, cv_size_jt, xbx, be, iv, xhat, cv, & re, y, j, grid, config_flags ) !--------------------------------------------------------------------------- @@ -15,6 +15,7 @@ subroutine da_calculate_j(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size_jp, integer, intent(in) :: cv_size_je ! Je cv size. integer, intent(in) :: cv_size_jp ! Jp cv size. integer, intent(in) :: cv_size_jl ! Jl cv size. + integer, intent(in) :: cv_size_jt ! Jt cv size. type (xbx_type),intent(inout) :: xbx ! For header & non-grid arrays. type (be_type), intent(in) :: be ! background error structure. type (iv_type), intent(inout) :: iv ! innovation vector (o-b). @@ -29,9 +30,11 @@ subroutine da_calculate_j(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size_jp, integer :: je_start, je_end ! Start/end indices of Je. integer :: jl_start, jl_end ! Start/end indices of Je. + integer :: jt_start, jt_end, iphase real :: jo_partial ! jo for this processor type (y_type) :: jo_grad_y ! Grad_y(jo) - real :: cv_xhat_jb(cv_size_jb), cv_xhat_je(cv_size_je), cv_xhat_jl(cv_size_jl) + real :: cv_xhat_jb(cv_size_jb), cv_xhat_je(cv_size_je), cv_xhat_jl(cv_size_jl), & + cv_xhat_jt(cv_size_jt) integer :: ndynopt real :: dtemp1x integer :: i, jj, k @@ -42,11 +45,13 @@ subroutine da_calculate_j(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size_jp, integer :: jp_start, jp_end ! Start/end indices of Jp. integer :: inst, ichan, npred, ipred, id real :: bgerr, gnorm_jp + real :: gnorm_jt integer :: n, cldtoplevel(1), icld, nclouds, ncv, minlev_cld real :: jd_local real :: js_local real :: jm_local + real :: jt_local real, allocatable :: cc(:) real :: inc_div(ims:ime,jms:jme,kms:kme) @@ -61,7 +66,9 @@ subroutine da_calculate_j(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size_jp, jp_start = cv_size_jb + cv_size_je + 1 jp_end = cv_size_jb + cv_size_je + cv_size_jp jl_start = cv_size_jb + cv_size_je + cv_size_jp + 1 - jl_end = cv_size_jb + cv_size_je + cv_size_jp + cv_size_jl + jl_end = cv_size_jb + cv_size_je + cv_size_jp + cv_size_jl + jt_start = jl_end + 1 + jt_end = jl_end + cv_size_jt inc_div = 0.0 @@ -314,10 +321,34 @@ subroutine da_calculate_j(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size_jp, end if !------------------------------------------------------------------------- - ! [7.0] calculate total cost function j = jo + jb + jc + je + jd + jp + js: + ! [n.n] calculate jt: !------------------------------------------------------------------------- - j % total = j % jb + j % jo % total + j % je + j % jd + j % jp + j % js + j % jm + j % jt = 0.0 + jt_local = 0.0 + if (cv_size_jt > 0) then + do n = 1, iv%varbc_tamdar%nair + do iphase = 1, iv%varbc_tamdar%nphase + if (iv%varbc_tamdar%nobs(iphase,n) > 0 .and. iv%varbc_tamdar%ifuse(iphase,n) > 0) then + do ipred = 1, iv%varbc_tamdar%npred + id = iv%varbc_tamdar%index(ipred,iphase,n) + bgerr = iv%varbc_tamdar%bgerr(ipred,iphase,n) + if (bgerr > 0.0) & + cv_xhat_jt(id-jt_start+1) = (1/sqrt(bgerr)) * & + SUM( (cv(id)+xhat(id)) * iv%varbc_tamdar%vtox(ipred,1:iv%varbc_tamdar%npred,iphase,n) ) + end do + end if + end do + end do + jt_local = 0.5 * da_dot(cv_size_jt, cv_xhat_jt, cv_xhat_jt) + end if + j % jt = wrf_dm_sum_real(jt_local) + + !------------------------------------------------------------------------------- + ! [7.0] calculate total cost function j = jo + jb + jc + je + jd + jp + js + jt: + !------------------------------------------------------------------------------- + + j % total = j % jb + j % jo % total + j % je + j % jd + j % jp + j % js + j % jm + j % jt if (grid%jcdfi_use) j % total = j % total + j % jc if (var4d) j % total = j % total + j % jl @@ -326,12 +357,12 @@ subroutine da_calculate_j(it, iter, cv_size, cv_size_jb, cv_size_je, cv_size_jp, !------------------------------------------------------------------------- if (rootproc) then if (it == 1 .and. iter == 0) then - write(unit=cost_unit,fmt='(a)')'Outer EPS Inner J Jb Jo Jc Je Jd Jp Js Jl Jm' + write(unit=cost_unit,fmt='(a)')'Outer EPS Inner J Jb Jo Jc Je Jd Jp Js Jl Jm Jt' write(unit=cost_unit,fmt='(a)')'Iter Iter ' end if - write(unit=cost_unit,fmt='(2x,i2,1x,e10.3,2x,i4,10(1x,f10.3))') & - it, EPS(it), iter, j % total, j % jb, j % jo % total, j % jc, j % je, j % jd, j % jp, j%js, j%jl, j%jm + write(unit=cost_unit,fmt='(2x,i2,1x,e10.3,2x,i4,11(1x,f10.3))') & + it, EPS(it), iter, j % total, j % jb, j % jo % total, j % jc, j % je, j % jd, j % jp, j%js, j%jl, j%jm, j%jt end if call da_deallocate_y (jo_grad_y) diff --git a/var/da/da_minimisation/da_get_innov_vector.inc b/var/da/da_minimisation/da_get_innov_vector.inc index 7c10528551..92f9e808de 100644 --- a/var/da/da_minimisation/da_get_innov_vector.inc +++ b/var/da/da_minimisation/da_get_innov_vector.inc @@ -196,6 +196,11 @@ subroutine da_get_innov_vector (it, num_qcstat_conv, ob, iv, grid, config_flags) end if #endif + if (.not. anal_type_verify .and. use_varbc_tamdar) then + call da_varbc_tamdar_direct(iv,ob) + if (it == 1) call da_varbc_tamdar_precond(iv) + end if + if ( use_rainobs .and. num_fgat_time > 1 .and. var4d) then deallocate (hr_rainc ) deallocate (hr_rainnc) diff --git a/var/da/da_minimisation/da_get_var_diagnostics.inc b/var/da/da_minimisation/da_get_var_diagnostics.inc index a52f61c443..b15ff07890 100644 --- a/var/da/da_minimisation/da_get_var_diagnostics.inc +++ b/var/da/da_minimisation/da_get_var_diagnostics.inc @@ -218,6 +218,8 @@ subroutine da_get_var_diagnostics(it, iv, j) write(unit=stdout,fmt='(a,f15.5)') ' Final value of Jp = ', j % jp write(unit=stdout,fmt='(a,f15.5)') ' Final value of Jl = ', j % jl write(unit=stdout,fmt='(a,f15.5)') ' Final value of Jm = ', j % jm + write(unit=stdout,fmt='(a,f15.5)') ' Final value of Jt = ', j % jt + if (num_stats_tot > 0) & write(unit=stdout,fmt='(a,f15.5)') ' Final J / total num_obs = ', j % total / & real(num_stats_tot) diff --git a/var/da/da_minimisation/da_kmat_mul.inc b/var/da/da_minimisation/da_kmat_mul.inc index 2cafc6983a..234cdd2715 100644 --- a/var/da/da_minimisation/da_kmat_mul.inc +++ b/var/da/da_minimisation/da_kmat_mul.inc @@ -44,7 +44,7 @@ subroutine da_kmat_mul(grid, config_flags, & call da_calculate_residual(iv,y,re) ! H^T.R^-1.[y^o-H(x_b)] - call da_calculate_gradj(1,1,cv_size,0,0,0,0,xbx,be,iv,cv,y,shat,grid,config_flags,re) + call da_calculate_gradj(1,1,cv_size,0,0,0,0,0,xbx,be,iv,cv,y,shat,grid,config_flags,re) shat = - shat !! Compensate for sign in calculation of grad_v (Jo) ! A.H^T.R^-1.[y^o-H(x_b)] diff --git a/var/da/da_minimisation/da_minimisation.f90 b/var/da/da_minimisation/da_minimisation.f90 index 9364ece06e..4de31e01bb 100644 --- a/var/da/da_minimisation/da_minimisation.f90 +++ b/var/da/da_minimisation/da_minimisation.f90 @@ -42,7 +42,7 @@ module da_minimisation write_filtered_rad,omb_set_rand,use_rad,var_scaling2,var_scaling1, & var_scaling4,var_scaling5,var_scaling3, jo_unit, test_gradient, & print_detail_grad,omb_set_rand,grad_unit,cost_unit, num_pseudo, cv_options, & - cv_size_domain_je,cv_size_domain_jb, cv_size_domain_jp, cv_size_domain_js, cv_size_domain_jl, & + cv_size_domain_je,cv_size_domain_jb, cv_size_domain_jp, cv_size_domain_js, cv_size_domain_jl, cv_size_domain_jt, & sound, mtgirs, sonde_sfc, synop, profiler, gpsref, gpseph, gpspw, polaramv, geoamv, ships, metar, & satem, radar, ssmi_rv, ssmi_tb, ssmt1, ssmt2, airsr, pilot, airep,tamdar, tamdar_sfc, rain, & bogus, buoy, qscat,pseudo, radiance, monitor_on, max_ext_its, use_rttov_kmatrix,& @@ -57,7 +57,8 @@ module da_minimisation write_detail_grad_fn, pseudo_uvtpq, lanczos_ep_filename, use_divc, divc_factor, & cloud_cv_options, use_cv_w, var_scaling6, var_scaling7, var_scaling8, var_scaling9, & var_scaling10, var_scaling11, & - write_gts_omb_oma, write_unpert_obs, write_rej_obs_conv, pseudo_time + write_gts_omb_oma, write_unpert_obs, write_rej_obs_conv, pseudo_time, & + use_varbc_tamdar, varbc_tamdar_nobsmin, varbc_tamdar_unit use da_define_structures, only : iv_type, y_type, j_type, be_type, & xbx_type, jo_type, da_allocate_y,da_zero_x,da_zero_y,da_deallocate_y, & da_zero_vp_type, qhat_type @@ -112,6 +113,8 @@ module da_minimisation da_ao_stats_tamdar_sfc, da_oi_stats_tamdar_sfc,da_oi_stats_tamdar_sfc, & da_get_innov_vector_tamdar_sfc, & da_jo_and_grady_tamdar_sfc, da_residual_tamdar_sfc + use da_varbc_tamdar, only : da_varbc_tamdar_tl,da_varbc_tamdar_adj, & + da_varbc_tamdar_direct,da_varbc_tamdar_precond #if defined(RTTOV) || defined(CRTM) use da_radiance, only : da_calculate_grady_rad, da_write_filtered_rad, & diff --git a/var/da/da_minimisation/da_minimise_cg.inc b/var/da/da_minimisation/da_minimise_cg.inc index 992837cde6..90dae0ab32 100644 --- a/var/da/da_minimisation/da_minimise_cg.inc +++ b/var/da/da_minimisation/da_minimise_cg.inc @@ -83,10 +83,10 @@ subroutine da_minimise_cg(grid, config_flags, & jp_end = be % cv % size_jb + be % cv % size_je + be % cv % size_jp call da_calculate_j(it, 0, cv_size, be%cv%size_jb, be%cv%size_je, be%cv%size_jp, & - be%cv%size_jl, xbx, be, iv, xhat, cv, re, y, j_cost, grid, config_flags) + be%cv%size_jl, be%cv%size_jt, xbx, be, iv, xhat, cv, re, y, j_cost, grid, config_flags) call da_calculate_gradj(it,0,cv_size,be%cv%size_jb, be%cv%size_je,be%cv%size_jp, & - be%cv%size_jl,xbx, be, iv, xhat+cv, y, ghat, grid, config_flags, re) + be%cv%size_jl,be%cv%size_jt,xbx, be, iv, xhat+cv, y, ghat, grid, config_flags, re) j0_total = j_cost%total if (j0_total == 0.0) then @@ -147,7 +147,7 @@ subroutine da_minimise_cg(grid, config_flags, & if (rrmold == 0.0) exit call da_calculate_gradj(it,iter,cv_size,be%cv%size_jb,be%cv%size_je,be%cv%size_jp, & - be%cv%size_jl,xbx,be,iv,phat,y,fhat,grid,config_flags) + be%cv%size_jl,be%cv%size_jt,xbx,be,iv,phat,y,fhat,grid,config_flags) apdotp = da_dot_cv(cv_size, fhat, phat, grid, be%cv_mz, be%ncv_mz, jp_start, jp_end) @@ -185,7 +185,7 @@ subroutine da_minimise_cg(grid, config_flags, & if (calculate_cg_cost_fn) then call da_calculate_j(it, iter, cv_size, be%cv%size_jb, be%cv%size_je, be%cv%size_jp, & - be%cv%size_jl, xbx, be, iv, xhat, cv, re, y, j_cost, grid, config_flags) + be%cv%size_jl, be%cv%size_jt, xbx, be, iv, xhat, cv, re, y, j_cost, grid, config_flags) j_total = j_cost%total else j_total = j0_total + 0.5 * da_dot_cv(cv_size,ghat0,xhat,grid,be%cv_mz,be%ncv_mz,jp_start,jp_end) @@ -193,7 +193,7 @@ subroutine da_minimise_cg(grid, config_flags, & if ( write_detail_grad_fn ) then call da_calculate_gradj(it, iter, cv_size, be%cv%size_jb, be%cv%size_je, be%cv%size_jp, & - be%cv%size_jl, xbx, be, iv, xhat+cv, y, fhat, grid, config_flags, re) + be%cv%size_jl, be%cv%size_jt, xbx, be, iv, xhat+cv, y, fhat, grid, config_flags, re) end if write(unit=stdout,fmt=12) " minimize_cg ", it, iter, j_total, rrmnew_norm, step @@ -225,10 +225,10 @@ subroutine da_minimise_cg(grid, config_flags, & write(unit=stdout,fmt='(A)') " " call da_calculate_j(it, iter, cv_size, be%cv%size_jb, be%cv%size_je, be%cv%size_jp, & - be%cv%size_jl, xbx, be, iv, xhat, cv, re, y, j_cost, grid, config_flags) + be%cv%size_jl, be%cv%size_jt, xbx, be, iv, xhat, cv, re, y, j_cost, grid, config_flags) call da_calculate_gradj(it,iter,cv_size,be%cv%size_jb, be%cv%size_je,be%cv%size_jp, & - be%cv%size_jl,xbx, be, iv, xhat+cv, y, ghat, grid, config_flags, re) + be%cv%size_jl,be%cv%size_jt,xbx, be, iv, xhat+cv, y, ghat, grid, config_flags, re) rrmnew_norm = SQRT(da_dot_cv(cv_size,ghat,ghat,grid,be%cv_mz,be%ncv_mz,jp_start,jp_end)) diff --git a/var/da/da_minimisation/da_minimise_lz.inc b/var/da/da_minimisation/da_minimise_lz.inc index 09d201759c..bbd8efbe6e 100644 --- a/var/da/da_minimisation/da_minimise_lz.inc +++ b/var/da/da_minimisation/da_minimise_lz.inc @@ -68,10 +68,10 @@ subroutine da_minimise_lz(grid, config_flags, & jp_end = be % cv % size_jb + be % cv % size_je + be % cv % size_jp call da_calculate_j(it, 0, cv_size, be%cv%size_jb, be%cv%size_je, be%cv%size_jp, & - be%cv%size_jl, xbx, be, iv, xhat, cv, re, y, j_cost, grid, config_flags) + be%cv%size_jl, be%cv%size_jt, xbx, be, iv, xhat, cv, re, y, j_cost, grid, config_flags) call da_calculate_gradj(it,0,cv_size,be%cv%size_jb,be%cv%size_je,be%cv%size_jp, & - be%cv%size_jl,xbx, be, iv, xhat+cv, y, ghat0, grid, config_flags, re) + be%cv%size_jl, be%cv%size_jt, xbx, be, iv, xhat+cv, y, ghat0, grid, config_flags, re) j0_total = j_cost%total if (j0_total == 0.0) then @@ -98,7 +98,7 @@ subroutine da_minimise_lz(grid, config_flags, & qhat(:,iter) = ghat / beta(iter-1) call da_calculate_gradj(it,iter,cv_size,be%cv%size_jb,be%cv%size_je,be%cv%size_jp, & - be%cv%size_jl, xbx,be,iv,qhat(:,iter),y,fhat,grid,config_flags) + be%cv%size_jl, be%cv%size_jt, xbx,be,iv,qhat(:,iter),y,fhat,grid,config_flags) ! Apply Lanczos recurrence and orthonormalize new gradient (using modified Gramm-Schmidt) !---------------------------------------------------------------------------------------- @@ -161,11 +161,11 @@ subroutine da_minimise_lz(grid, config_flags, & ! ---------------------------------------------- if (calculate_cg_cost_fn) & call da_calculate_j(it, iter, cv_size, be%cv%size_jb, be%cv%size_je, be%cv%size_jp, & - be%cv%size_jl, xbx, be, iv, xhat, cv, re, y, j_cost, grid, config_flags) + be%cv%size_jl, be%cv%size_jt, xbx, be, iv, xhat, cv, re, y, j_cost, grid, config_flags) if ( write_detail_grad_fn ) then call da_calculate_gradj(it, iter, cv_size, be%cv%size_jb, be%cv%size_je, be%cv%size_jp, & - be%cv%size_jl, xbx, be, iv, xhat+cv, y, fhat, grid, config_flags, re) + be%cv%size_jl, be%cv%size_jt, xbx, be, iv, xhat+cv, y, fhat, grid, config_flags, re) end if end do @@ -182,10 +182,10 @@ subroutine da_minimise_lz(grid, config_flags, & write(unit=stdout,fmt='(A)') " " call da_calculate_j(it, iter, cv_size, be%cv%size_jb, be%cv%size_je, be%cv%size_jp, & - be%cv%size_jl, xbx, be, iv, xhat, cv, re, y, j_cost, grid, config_flags) + be%cv%size_jl, be%cv%size_jt, xbx, be, iv, xhat, cv, re, y, j_cost, grid, config_flags) call da_calculate_gradj(it,iter,cv_size,be%cv%size_jb,be%cv%size_je,be%cv%size_jp, & - be%cv%size_jl,xbx, be, iv, xhat+cv, y, ghat, grid, config_flags, re) + be%cv%size_jl, be%cv%size_jt, xbx, be, iv, xhat+cv, y, ghat, grid, config_flags, re) write(unit=stdout,fmt=15) iter, j_cost%total , & SQRT(da_dot_cv(cv_size, ghat, ghat, grid, be%cv_mz, be%ncv_mz, jp_start, jp_end)) diff --git a/var/da/da_minimisation/da_transform_vtoy.inc b/var/da/da_minimisation/da_transform_vtoy.inc index c754bd11f7..5a54a715d4 100644 --- a/var/da/da_minimisation/da_transform_vtoy.inc +++ b/var/da/da_minimisation/da_transform_vtoy.inc @@ -210,6 +210,9 @@ subroutine da_transform_vtoy(cv_size, be, ep, cv, iv, vp, vv, xbx, & #if defined(RTTOV) || defined(CRTM) if (use_varbc) call da_varbc_tl(cv_size, cv, iv, y) #endif + + if (use_varbc_tamdar) call da_varbc_tamdar_tl(cv_size, cv, iv, y) + if (trace_use) call da_trace_exit("da_transform_vtoy") end subroutine da_transform_vtoy diff --git a/var/da/da_minimisation/da_transform_vtoy_adj.inc b/var/da/da_minimisation/da_transform_vtoy_adj.inc index e67b2017e8..3001574498 100644 --- a/var/da/da_minimisation/da_transform_vtoy_adj.inc +++ b/var/da/da_minimisation/da_transform_vtoy_adj.inc @@ -255,6 +255,9 @@ subroutine da_transform_vtoy_adj(cv_size, be, ep, cv, iv, vp, vv, xbx, y, & #if defined(RTTOV) || defined(CRTM) if (use_varbc) call da_varbc_adj(cv_size, cv, iv, y) #endif + + if (use_varbc_tamdar) call da_varbc_tamdar_adj(cv_size, cv, iv, y) + if (trace_use) call da_trace_exit("da_transform_vtoy_adj") end subroutine da_transform_vtoy_adj diff --git a/var/da/da_test/da_check_gradient.inc b/var/da/da_test/da_check_gradient.inc index 399859816f..5b88009fd2 100644 --- a/var/da/da_test/da_check_gradient.inc +++ b/var/da/da_test/da_check_gradient.inc @@ -62,10 +62,10 @@ subroutine da_check_gradient(grid, config_flags, cv_size, xhat, cv, pdx, itertes yhat=xhat call da_calculate_j(1, 1, cv_size, be%cv%size_jb, be%cv%size_je, be%cv%size_jp, & - be%cv%size_jl, xbx, be, iv, yhat, cv, re, y, j_cost, grid, config_flags) + be%cv%size_jl, be%cv%size_jt, xbx, be, iv, yhat, cv, re, y, j_cost, grid, config_flags) call da_calculate_gradj(1,1,cv_size,be%cv%size_jb,be%cv%size_je,be%cv%size_jp, & - be%cv%size_jl,xbx, be, iv, yhat+cv, y, grad, grid, config_flags, re) + be%cv%size_jl, be%cv%size_jt, xbx, be, iv, yhat+cv, y, grad, grid, config_flags, re) zfy = j_cost%total @@ -110,10 +110,10 @@ subroutine da_check_gradient(grid, config_flags, cv_size, xhat, cv, pdx, itertes end do call da_calculate_j(1, 1, cv_size, be%cv%size_jb, be%cv%size_je, be%cv%size_jp, & - be%cv%size_jl, xbx, be, iv, yhat, cv, re, y, j_cost, grid, config_flags) + be%cv%size_jl, be%cv%size_jt, xbx, be, iv, yhat, cv, re, y, j_cost, grid, config_flags) call da_calculate_gradj(1,1,cv_size,be%cv%size_jb,be%cv%size_je,be%cv%size_jp, & - be%cv%size_jl,xbx, be, iv, yhat+cv, y, grad, grid, config_flags, re) + be%cv%size_jl, be%cv%size_jt, xbx, be, iv, yhat+cv, y, grad, grid, config_flags, re) zfy = j_cost%total diff --git a/var/da/da_varbc_tamdar/da_varbc_tamdar.f90 b/var/da/da_varbc_tamdar/da_varbc_tamdar.f90 new file mode 100755 index 0000000000..84014de22e --- /dev/null +++ b/var/da/da_varbc_tamdar/da_varbc_tamdar.f90 @@ -0,0 +1,32 @@ +module da_varbc_tamdar + + !--------------------------------------------------------------------------- + ! Purpose: module for variational bias correction of TAMDAR temperature. + !--------------------------------------------------------------------------- + + use module_dm, only : wrf_dm_sum_real, wrf_dm_sum_integer + use da_control, only : trace_use,missing_r, qc_varbc_bad, stdout, filename_len, & + rootproc, varbc_nbgerr, ierr, comm, max_ext_its, & + obs_qc_pointer, tamdar, tamdar_sfc, missing_r,print_detail_obs, & + use_varbc_tamdar, varbc_tamdar_nbgerr, varbc_tamdar_nobsmin, & + varbc_tamdar_unit, varbc_tamdar_bm, varbc_tamdar_pred0, fail_varbc_aircraft + use da_define_structures, only : iv_type, y_type, be_type, varbc_tamdar_type + use da_reporting, only : da_error, message, da_warning, da_message + use da_tools, only : da_eof_decomposition, da_diff_seconds + use da_tools_serial, only : da_free_unit, da_get_unit + use da_tracing, only : da_trace_entry, da_trace_exit, da_trace, & + da_trace_int_sort + + implicit none + +contains + +#include "da_varbc_tamdar_init.inc" +#include "da_varbc_tamdar_pred.inc" +#include "da_varbc_tamdar_direct.inc" +#include "da_varbc_tamdar_precond.inc" +#include "da_varbc_tamdar_tl.inc" +#include "da_varbc_tamdar_adj.inc" +#include "da_varbc_tamdar_update.inc" + +end module da_varbc_tamdar diff --git a/var/da/da_varbc_tamdar/da_varbc_tamdar_adj.inc b/var/da/da_varbc_tamdar/da_varbc_tamdar_adj.inc new file mode 100755 index 0000000000..0ebf7d5353 --- /dev/null +++ b/var/da/da_varbc_tamdar/da_varbc_tamdar_adj.inc @@ -0,0 +1,53 @@ + subroutine da_varbc_tamdar_adj (cv_size, cv, iv, y) + + implicit none + + integer, intent(in) :: cv_size + real, intent(inout) :: cv(1:cv_size) + type (iv_type), intent(inout) :: iv + type (y_type), intent(in) :: y ! y = h (xa) + + integer :: isn,iflt,iobs,ipred,iphase,npred + real, allocatable :: varbc_param_adj(:) + real :: yh, cv_local + + + if (trace_use) call da_trace_entry("da_varbc_tamdar_adj") + + npred = iv%varbc_tamdar%npred + + allocate( varbc_param_adj(npred) ) + + do iflt = 1, iv%varbc_tamdar%nair + do iphase = 1, iv%varbc_tamdar%nphase + if (iv%varbc_tamdar%nobs_sum(iphase,iflt) >= varbc_tamdar_nobsmin) then + + varbc_param_adj(:) = 0.0 + + if (iv%varbc_tamdar%nobs(iphase,iflt) > 0 .and. iv%varbc_tamdar%ifuse(iphase,iflt) > 0) then + do iobs = 1, iv%varbc_tamdar%nobs(iphase,iflt) + isn = iv%varbc_tamdar%obs_sn(iobs,iphase,iflt) + yh = 0 + if( iv%tamdar(isn)%t(1)%qc >= 0 ) yh = y%tamdar(isn)%t(1) + varbc_param_adj(1:npred) = varbc_param_adj(1:npred) + & + yh * iv%varbc_tamdar%pred(1:npred,iphase,iflt) + end do + end if + + do ipred = 1, npred + cv_local = SUM( varbc_param_adj(1:npred) * iv%varbc_tamdar%vtox(ipred,1:npred,iphase,iflt) ) + cv_local = wrf_dm_sum_real(cv_local) + + if (iv%varbc_tamdar%nobs(iphase,iflt) > 0 .and. iv%varbc_tamdar%ifuse(iphase,iflt) > 0) & + cv(iv%varbc_tamdar%index(ipred,iphase,iflt)) = cv( iv%varbc_tamdar%index(ipred,iphase,iflt) ) + cv_local + end do + + end if + end do + end do + + deallocate(varbc_param_adj) + + if (trace_use) call da_trace_exit("da_varbc_tamdar_adj") + + end subroutine da_varbc_tamdar_adj diff --git a/var/da/da_varbc_tamdar/da_varbc_tamdar_direct.inc b/var/da/da_varbc_tamdar/da_varbc_tamdar_direct.inc new file mode 100755 index 0000000000..ce824f23d6 --- /dev/null +++ b/var/da/da_varbc_tamdar/da_varbc_tamdar_direct.inc @@ -0,0 +1,94 @@ + subroutine da_varbc_tamdar_direct (iv,ob) + + !-----------------------------------------------! + ! Apply bias correction to TAMDAR innovations ! + !-----------------------------------------------! + + implicit none + + type (iv_type), intent(inout) :: iv + type (y_type), intent(inout) :: ob + + real :: bc,bias + real :: contri(5) + integer :: i,isn,iflt,iobs,ipred,iphase,nobs + character(len=3) :: cphz(3) + character(len=30) :: stringn + + + if (trace_use) call da_trace_entry("da_varbc_tamdar_direct") + + if (rootproc) & + write(unit=varbc_tamdar_unit,fmt='(//5X,A/)')'Calculating corrected innovation' + + cphz = (/'des','asc','cru'/) + + do iflt = 1, iv%varbc_tamdar%nair + do iphase = 1, iv%varbc_tamdar%nphase + if (iv%varbc_tamdar%nobs_sum(iphase,iflt) >= varbc_tamdar_nobsmin) then + + bc = SUM( iv%varbc_tamdar%param(1:iv%varbc_tamdar%npred,iphase,iflt) * & + iv%varbc_tamdar%pred(1:iv%varbc_tamdar%npred,iphase,iflt) ) + + nobs = 0 + bias = 0. + + if (iv%varbc_tamdar%nobs(iphase,iflt) > 0 .and. iv%varbc_tamdar%ifuse(iphase,iflt) > 0 ) then + do iobs = 1, iv%varbc_tamdar%nobs(iphase,iflt) + isn = iv%varbc_tamdar%obs_sn(iobs,iphase,iflt) + if (iv%tamdar(isn)%t(1)%qc >= 0) then + bias = bias + iv%tamdar(isn)%t(1)%inv + nobs = nobs + 1 + + iv%tamdar(isn)%t(1)%inv = iv%tamdar(isn)%t(1)%inv - bc + end if + end do + end if + + if (wrf_dm_sum_integer(nobs) > 0) & + bias = wrf_dm_sum_real(bias)/wrf_dm_sum_integer(nobs) + + if (rootproc .and. iv%varbc_tamdar%ifuse(iphase,iflt) > 0) then + write(unit=varbc_tamdar_unit,fmt='(10X,A,I5,A,2X,A,2X,A,2F10.4)') & + 'bias corrected for ',iv%varbc_tamdar%tail_id(iflt), & + ' at', cphz(iphase), ': (BC,OMB)', bc, bias + end if + + end if + end do + end do + + if (rootproc) write(unit=varbc_tamdar_unit,fmt='(/10X,A)') & + " ID Phase Predictor/Parameter/Bias model (npred) & BC" + + ipred = 3*iv%varbc_tamdar%npred+1 + write(stringn,'(I3)') ipred + stringn = '(10X,I4,2X,A,2X,'//trim(ADJUSTL(stringn))//'f8.3)' + stringn = trim(adjustl(stringn)) + + do iflt = 1, iv%varbc_tamdar%nair + do iphase = 1, iv%varbc_tamdar%nphase + if (rootproc .and. iv%varbc_tamdar%ifuse(iphase,iflt) > 0) then + + contri(:) = 0. + bc = 0. + + do ipred = 1, iv%varbc_tamdar%npred + contri(ipred) = & + iv%varbc_tamdar%param(ipred,iphase,iflt) * & + iv%varbc_tamdar%pred(ipred,iphase,iflt) + bc = bc + contri(ipred) + end do + + write(unit=varbc_tamdar_unit,fmt=stringn) & + iv%varbc_tamdar%tail_id(iflt), cphz(iphase), & + (iv%varbc_tamdar%pred(ipred,iphase,iflt), & + iv%varbc_tamdar%param(ipred,iphase,iflt), & + contri(ipred), ipred=1,iv%varbc_tamdar%npred),bc + end if + end do + end do + + if (trace_use) call da_trace_exit("da_varbc_tamdar_direct") + + end subroutine da_varbc_tamdar_direct diff --git a/var/da/da_varbc_tamdar/da_varbc_tamdar_init.inc b/var/da/da_varbc_tamdar/da_varbc_tamdar_init.inc new file mode 100755 index 0000000000..733c468440 --- /dev/null +++ b/var/da/da_varbc_tamdar/da_varbc_tamdar_init.inc @@ -0,0 +1,137 @@ + subroutine da_varbc_tamdar_init (iv) + + !-----------------------------------------------------! + ! Initialize Varbc for TAMDAR from VARBC_TAMDAR.tbl ! + !-----------------------------------------------------! + + implicit none + + type (iv_type), intent(inout) :: iv + + integer, parameter :: nmaxpred = 5 ! 1,w,(optional: Mach,dT/dt,T) + integer, parameter :: nphase = 3 ! descent/ascent/cruise + + integer :: npred,nair + integer :: iunit,i,id,iflt,iobs,iphase + + logical :: ifexist_table + character(len=120) :: filename + character(len=5), allocatable :: tail_id(:) + + + if (trace_use) call da_trace_entry("da_varbc_tamdar_init") + + if (rootproc) then + write(unit=varbc_tamdar_unit,fmt='(//5X,A/)') 'VARBC_TAMDAR namelist options' + write(unit=varbc_tamdar_unit,fmt='(10X,A,7X,L)') 'use_varbc_tamdar : ', use_varbc_tamdar + write(unit=varbc_tamdar_unit,fmt='(10X,A,7X,I1)') 'varbc_tamdar_bm : ', varbc_tamdar_bm ! bias models + write(unit=varbc_tamdar_unit,fmt='(10X,A,F8.3)') 'varbc_tamdar_pred0 : ', varbc_tamdar_pred0 ! predictor(1) + write(unit=varbc_tamdar_unit,fmt='(10X,A,I8)') 'varbc_tamdar_nbgerr : ', varbc_tamdar_nbgerr + write(unit=varbc_tamdar_unit,fmt='(10X,A,I8)') 'varbc_tamdar_nobsmin : ', varbc_tamdar_nobsmin + end if + + filename = 'VARBC_TAMDAR.tbl' + inquire(file=trim(adjustl(filename)), exist=ifexist_table) + + if (ifexist_table) then + + if (rootproc) write(unit=varbc_tamdar_unit,fmt='(//5X,A/)') & + 'Reading in VARBC_TAMDAR.tbl and initializating' + + iv%varbc_tamdar%nmaxpred = nmaxpred + iv%varbc_tamdar%nphase = nphase + + call da_get_unit(iunit) + open(unit=iunit,file=filename, status='old') + do i = 1, 11; read(iunit,*); end do + read(iunit,*) nair + read(iunit,*) + + iv%varbc_tamdar%nair = nair + + npred = 2 + if (varbc_tamdar_bm == 2) npred = 5 + + iv%varbc_tamdar%npred = npred + + allocate ( iv%varbc_tamdar%tail_id (nair) ) + allocate ( iv%varbc_tamdar%ifuse (nphase,nair) ) + allocate ( iv%varbc_tamdar%nobs (nphase+1,nair) ) + allocate ( iv%varbc_tamdar%nobs_sum (nphase+1,nair) ) + allocate ( iv%varbc_tamdar%param (nmaxpred,nphase,nair) ) + allocate ( iv%varbc_tamdar%pred (nmaxpred,nphase,nair) ) + allocate ( iv%varbc_tamdar%bgerr (npred,nphase,nair) ) + allocate ( iv%varbc_tamdar%index (npred,nphase,nair) ) + allocate ( iv%varbc_tamdar%vtox (npred,npred,nphase,nair) ) + + iv%varbc_tamdar%tail_id(:) = 0 + iv%varbc_tamdar%ifuse(:,:) = 0 + iv%varbc_tamdar%nobs(:,:) = 0 + iv%varbc_tamdar%nobs_sum(:,:) = 0 + iv%varbc_tamdar%param(:,:,:) = 0. + iv%varbc_tamdar%pred(:,:,:) = 0. + iv%varbc_tamdar%pred(1,:,:) = varbc_tamdar_pred0 + iv%varbc_tamdar%bgerr(:,:,:) = 0. + iv%varbc_tamdar%index(:,:,:) = 0 + iv%varbc_tamdar%vtox(:,:,:,:) = 0. + + + write(iv%varbc_tamdar%fmt_param,*) '(1X,A5,1X,3(1X,I1),15F7.3))' + + allocate( tail_id (nair) ) + tail_id = '' + do iflt = 1, nair + read(iunit, fmt=trim(adjustl(iv%varbc_tamdar%fmt_param))) & + tail_id(iflt), iv%varbc_tamdar%ifuse(1:nphase,iflt), & + (iv%varbc_tamdar%param(1:nmaxpred,iphase,iflt), iphase=1,nphase) + + read(tail_id(iflt),'(I5)') iv%varbc_tamdar%tail_id(iflt) + + do iobs = 1, iv%info(tamdar)%nlocal + read(iv%info(tamdar)%id(iobs),'(I5)') id + if (iv%varbc_tamdar%tail_id(iflt) .eq. id) & + iv%varbc_tamdar%nobs(nphase+1,iflt) = iv%varbc_tamdar%nobs(nphase+1,iflt) + 1 + end do + + iv%varbc_tamdar%nobs_sum(nphase+1,iflt) = wrf_dm_sum_integer(iv%varbc_tamdar%nobs(nphase+1,iflt)) + end do + + iv%varbc_tamdar%nmaxobs = MAXVAL( iv%varbc_tamdar%nobs_sum(nphase+1,:) ) + + allocate ( iv%varbc_tamdar%obs_sn(iv%varbc_tamdar%nmaxobs,nphase+1,nair) ) + iv%varbc_tamdar%obs_sn(:,:,:) = 0 + + do iflt = 1, nair + i = 0 + do iobs = 1, iv%info(tamdar)%nlocal + read(iv%info(tamdar)%id(iobs),'(I5)') id + if (iv%varbc_tamdar%tail_id(iflt) .eq. id) then + i = i + 1 + iv%varbc_tamdar%obs_sn(i,nphase+1,iflt) = iobs + end if + end do + end do + + if (rootproc) then + write(unit=varbc_tamdar_unit,fmt='(10X,A,I4)') & + 'Number of aircrafts in table : ', nair + write(unit=varbc_tamdar_unit,fmt='(10X,A,I4)') & + 'Max. Obs. number in aircraft : ', iv%varbc_tamdar%nmaxobs + end if + + close(iunit) + call da_free_unit(iunit) + + deallocate (tail_id) + + else + + use_varbc_tamdar = .false. + if (rootproc) write(unit=varbc_tamdar_unit,fmt='(/5X,A/)') & + 'Could not find VARBC_TAMDAR.tbl file. VARBC_TAMDAR is switched off.' + + end if + + if (trace_use) call da_trace_exit("da_varbc_tamdar_init") + + end subroutine da_varbc_tamdar_init diff --git a/var/da/da_varbc_tamdar/da_varbc_tamdar_precond.inc b/var/da/da_varbc_tamdar/da_varbc_tamdar_precond.inc new file mode 100755 index 0000000000..0ca71d30c4 --- /dev/null +++ b/var/da/da_varbc_tamdar/da_varbc_tamdar_precond.inc @@ -0,0 +1,113 @@ + subroutine da_varbc_tamdar_precond (iv) + + !----------------------------------------------------! + ! Estimate Hessian for Preconditioning ! + ! ! + ! Preconditioning = inverse square root of Hessian ! + ! ! + ! Hermite real matrix: ! + ! A = Q*D*Q^T ! + ! A^(1/k) = Q*D^(1/K)*Q^T ! + !----------------------------------------------------! + + implicit none + + type (iv_type), intent (inout) :: iv + + integer :: i,j,iflt,iobs,isn,ipred,iphase,npred + real :: lhessian,lbgerr,predi,predj,verr + + real, allocatable :: hessian(:,:) + real*8, allocatable :: eignvec(:,:), eignval(:) + + + if (trace_use) call da_trace_entry("da_varbc_tamdar_precond") + + if (rootproc) & + write(unit=varbc_tamdar_unit,fmt='(//5X,A)') 'Estimating hessian for preconditioning' + + npred = iv%varbc_tamdar%npred + + allocate ( hessian(npred, npred) ) + allocate ( eignvec(npred, npred) ) + allocate ( eignval(npred) ) + + do iflt = 1, iv%varbc_tamdar%nair + do iphase = 1, iv%varbc_tamdar%nphase + if (iv%varbc_tamdar%nobs_sum(iphase,iflt) >= varbc_tamdar_nobsmin) then + + hessian(:,:) = 0. + eignvec(:,:) = 0. + eignval(:) = 0. + + lbgerr = 0.0 + if (iv%varbc_tamdar%nobs(iphase,iflt) > 0 .and. iv%varbc_tamdar%ifuse(iphase,iflt) > 0) then + do iobs = 1, iv%varbc_tamdar%nobs(iphase,iflt) + isn = iv%varbc_tamdar%obs_sn(iobs,iphase,iflt) + verr = 0.0 + if (iv%tamdar(isn)%t(1)%qc >= 0) verr = iv%tamdar(isn)%t(1)%error + lbgerr = lbgerr + verr**2.0/varbc_tamdar_nbgerr + end do + end if + lbgerr = wrf_dm_sum_real(lbgerr) + + do i = 1, npred + do j = i, npred + lhessian = 0.0 + + if (iv%varbc_tamdar%nobs(iphase,iflt) > 0 .and. iv%varbc_tamdar%ifuse(iphase,iflt) > 0) then + do iobs = 1, iv%varbc_tamdar%nobs(iphase,iflt) + isn = iv%varbc_tamdar%obs_sn(iobs,iphase,iflt) + predi = iv%varbc_tamdar%pred(i,iphase,iflt) + predj = iv%varbc_tamdar%pred(j,iphase,iflt) + + verr = 0.0 + if (iv%tamdar(isn)%t(1)%qc >= 0) verr = iv%tamdar(isn)%t(1)%error + if (verr > 0.) lhessian = lhessian + predi*predj/verr**2.0 + end do + end if + hessian(i,j) = wrf_dm_sum_real(lhessian) + hessian(j,i) = hessian(i,j) + end do + + iv%varbc_tamdar%bgerr(i,iphase,iflt) = lbgerr/iv%varbc_tamdar%nobs_sum(iphase,iflt) + if (iv%varbc_tamdar%bgerr(i,iphase,iflt) > 0.) & + hessian(i,i)=hessian(i,i)+1.0/iv%varbc_tamdar%bgerr(i,iphase,iflt) + end do + + if (iv%varbc_tamdar%nobs(iphase,iflt) > 0 .and. iv%varbc_tamdar%ifuse(iphase,iflt) > 0) then + + call da_eof_decomposition(npred, hessian(1:npred,1:npred), & + eignvec(1:npred,1:npred),eignval(1:npred)) + + if (ANY( eignval(1:npred) <= 0 )) then + write(unit=stdout,fmt='(A,I8,A,2F12.5)') & + 'VARBC_TAMDAR: non-positive Hessian for tail_id ', iv%varbc_tamdar%tail_id(iflt), & + '. Eigenvalues =',eignval(1:npred) + + do i = 1, npred + if (hessian(i,i) > 0) iv%varbc_tamdar%vtox(i,i,iphase,iflt) = 1.0/sqrt(hessian(i,i)) + end do + else + do i = 1, npred + do j = i, npred + iv%varbc_tamdar%vtox(i,j,iphase,iflt) = SUM( eignvec(i,1:npred)* & + sqrt(1.0/eignval(1:npred)) * & + eignvec(j,1:npred) ) + iv%varbc_tamdar%vtox(j,i,iphase,iflt) = iv%varbc_tamdar%vtox(i,j,iphase,iflt) + end do + end do + end if + end if + end if + end do + end do + + deallocate(hessian, eignvec, eignval) + + if (rootproc) & + write(unit=varbc_tamdar_unit,fmt='(/5X,A)') 'End estimating hessian for preconditioning' + + if (trace_use) call da_trace_exit("da_varbc_tamdar_precond") + + end subroutine da_varbc_tamdar_precond diff --git a/var/da/da_varbc_tamdar/da_varbc_tamdar_pred.inc b/var/da/da_varbc_tamdar/da_varbc_tamdar_pred.inc new file mode 100755 index 0000000000..720c682fc6 --- /dev/null +++ b/var/da/da_varbc_tamdar/da_varbc_tamdar_pred.inc @@ -0,0 +1,231 @@ + subroutine da_varbc_tamdar_pred(iv, be, ob) + + !--------------------------------------! + ! Calculate Predictors in Bias Model ! + !--------------------------------------! + + implicit none + + type (iv_type), intent(inout) :: iv + type (be_type), intent(inout) :: be + type (y_type), intent(inout) :: ob + + integer :: jt_start,size_jt + integer :: n,nphase,iflt,iobs,ipred,isn,iqc,iphase + real :: dhdt,dtdt,p1,p2,p3,p4 + + integer, allocatable :: tdiff(:),osn(:) + real , allocatable :: op(:),ot(:),oh(:),mach(:) + + character(len=24) :: otime,stime + character(len=40) :: stringn + + + if (trace_use) call da_trace_entry("da_varbc_tamdar_pred") + + if (rootproc) write(unit=varbc_tamdar_unit,fmt='(//5X,A/)') & + 'Calculating predictors' + + size_jt = 0 + jt_start = be%cv%size_jb + be%cv%size_je + be%cv%size_jp + be%cv%size_js + be%cv%size_jl + + nphase = iv%varbc_tamdar%nphase + + allocate( op (iv%varbc_tamdar%nmaxobs) ) + allocate( ot (iv%varbc_tamdar%nmaxobs) ) + allocate( oh (iv%varbc_tamdar%nmaxobs) ) + allocate( osn (iv%varbc_tamdar%nmaxobs) ) + allocate( mach (iv%varbc_tamdar%nmaxobs) ) + allocate( tdiff (iv%varbc_tamdar%nmaxobs) ) + + do iflt = 1, iv%varbc_tamdar%nair + if (iv%varbc_tamdar%nobs_sum(nphase+1,iflt) >= varbc_tamdar_nobsmin) then + + n = 0 + op(:) = 0. + ot(:) = 0. + oh(:) = 0. + osn(:) = 0 + mach(:) = 0. + tdiff(:) = 0 + + if (iv%varbc_tamdar%nobs(nphase+1,iflt) > 1) then + do iobs = 1, iv%varbc_tamdar%nobs(nphase+1,iflt) + + isn = iv%varbc_tamdar%obs_sn(iobs,nphase+1,iflt) + iqc = iv%tamdar(isn)%t(1)%qc + + if (iv%tamdar(isn)%p(1) > missing_r .and. iqc >= 0) then + n = n + 1 + osn(n) = isn + op(n) = iv%tamdar(isn)%p(1) + ot(n) = ob%tamdar(isn)%t(1) + oh(n) = (1-(op(n)/101325.0)**0.190284)*145366.45*0.3048 ! pressure altitude + + if (n == 1) stime = iv%info(tamdar)%date_char(isn) + otime = iv%info(tamdar)%date_char(isn) + tdiff(n) = int(da_diff_seconds(otime, stime)) + + ! May assign Mach number on the position of 'pw' in ob.ascii + if (op(n) >= 95000.0) mach(n) = 0.208 + if (op(n) < 95000.0 .and. op(n) >= 90000.0) mach(n) = 0.265 + if (op(n) < 90000.0 .and. op(n) >= 85000.0) mach(n) = 0.288 + if (op(n) < 85000.0 .and. op(n) >= 80000.0) mach(n) = 0.319 + if (op(n) < 80000.0 .and. op(n) >= 75000.0) mach(n) = 0.351 + if (op(n) < 75000.0 .and. op(n) >= 70000.0) mach(n) = 0.387 + if (op(n) < 70000.0 .and. op(n) >= 65000.0) mach(n) = 0.419 + if (op(n) < 65000.0 .and. op(n) >= 60000.0) mach(n) = 0.446 + if (op(n) < 60000.0 .and. op(n) >= 55000.0) mach(n) = 0.474 + if (op(n) < 55000.0 .and. op(n) >= 50000.0) mach(n) = 0.492 + if (op(n) < 50000.0 .and. op(n) >= 45000.0) mach(n) = 0.505 + if (op(n) < 45000.0 .and. op(n) >= 40000.0) mach(n) = 0.529 + if (op(n) < 40000.0 .and. op(n) >= 30000.0) mach(n) = 0.552 + if (op(n) < 30000.0) mach(n) = 0.580 + end if + end do + end if + + if (n > 1) then + do iobs = 2, n + + if (tdiff(iobs) == tdiff(iobs-1) .or. abs(tdiff(iobs) - tdiff(iobs-1)) > 360) cycle + + dhdt = (oh(iobs)-oh(iobs-1))/10.0/(tdiff(iobs)-tdiff(iobs-1)) + dtdt = -1.0*(ot(iobs)-ot(iobs-1))/(tdiff(iobs)-tdiff(iobs-1)) + + p1 = (dhdt + 2.0) / 5.0 + p2 = 1.0 / (10.0 * mach(iobs)) + p3 = dtdt * 5.0 + p4 = ot(iobs) / 2000.0 + + if (dhdt < -0.2) then + + iv%varbc_tamdar%nobs(1,iflt) = iv%varbc_tamdar%nobs(1,iflt) + 1 + iv%varbc_tamdar%obs_sn(iv%varbc_tamdar%nobs(1,iflt),1,iflt) = osn(iobs) + + if (varbc_tamdar_bm == 1) then + iv%varbc_tamdar%pred(2,1,iflt) = iv%varbc_tamdar%pred(2,1,iflt) + dhdt + elseif (varbc_tamdar_bm == 2) then + iv%varbc_tamdar%pred(2,1,iflt) = iv%varbc_tamdar%pred(2,1,iflt) + p1 + iv%varbc_tamdar%pred(3,1,iflt) = iv%varbc_tamdar%pred(3,1,iflt) + p2 + iv%varbc_tamdar%pred(4,1,iflt) = iv%varbc_tamdar%pred(4,1,iflt) + p3 + iv%varbc_tamdar%pred(5,1,iflt) = iv%varbc_tamdar%pred(5,1,iflt) + p4 + end if + + elseif (dhdt > 0.3) then + + iv%varbc_tamdar%nobs(2,iflt) = iv%varbc_tamdar%nobs(2,iflt) + 1 + iv%varbc_tamdar%obs_sn(iv%varbc_tamdar%nobs(2,iflt),2,iflt) = osn(iobs) + + if (varbc_tamdar_bm == 1) then + iv%varbc_tamdar%pred(2,2,iflt) = iv%varbc_tamdar%pred(2,2,iflt) + dhdt + elseif (varbc_tamdar_bm == 2) then + iv%varbc_tamdar%pred(2,2,iflt) = iv%varbc_tamdar%pred(2,2,iflt) + p1 + iv%varbc_tamdar%pred(3,2,iflt) = iv%varbc_tamdar%pred(3,2,iflt) + p2 + iv%varbc_tamdar%pred(4,2,iflt) = iv%varbc_tamdar%pred(4,2,iflt) + p3 + iv%varbc_tamdar%pred(5,2,iflt) = iv%varbc_tamdar%pred(5,2,iflt) + p4 + end if + + else + + iv%varbc_tamdar%nobs(3,iflt) = iv%varbc_tamdar%nobs(3,iflt) + 1 + iv%varbc_tamdar%obs_sn(iv%varbc_tamdar%nobs(3,iflt),3,iflt) = osn(iobs) + + if (varbc_tamdar_bm == 1) then + iv%varbc_tamdar%pred(2,3,iflt) = iv%varbc_tamdar%pred(2,3,iflt) + dhdt + elseif (varbc_tamdar_bm == 2) then + iv%varbc_tamdar%pred(2,3,iflt) = iv%varbc_tamdar%pred(2,3,iflt) + p1 + iv%varbc_tamdar%pred(3,3,iflt) = iv%varbc_tamdar%pred(3,3,iflt) + p2 + iv%varbc_tamdar%pred(4,3,iflt) = iv%varbc_tamdar%pred(4,3,iflt) + p3 + iv%varbc_tamdar%pred(5,3,iflt) = iv%varbc_tamdar%pred(5,3,iflt) + p4 + end if + + end if + + end do + end if + + do iphase = 1, nphase + + iv%varbc_tamdar%nobs_sum(iphase,iflt) = wrf_dm_sum_integer(iv%varbc_tamdar%nobs(iphase,iflt)) + + if (iv%varbc_tamdar%nobs_sum(iphase,iflt) >= varbc_tamdar_nobsmin) then + + do ipred = 1, iv%varbc_tamdar%npred + iv%varbc_tamdar%pred(ipred,iphase,iflt) = & + wrf_dm_sum_real(iv%varbc_tamdar%pred(ipred,iphase,iflt)) / & + iv%varbc_tamdar%nobs_sum(iphase,iflt) + if (ipred == 1) iv%varbc_tamdar%pred(ipred,iphase,iflt) = varbc_tamdar_pred0 + end do + + if (iv%varbc_tamdar%nobs(iphase,iflt) > 0 .and. iv%varbc_tamdar%ifuse(iphase,iflt) > 0) then + do ipred = 1, iv%varbc_tamdar%npred + size_jt = size_jt + 1 + iv%varbc_tamdar%index(ipred,iphase,iflt) = jt_start + size_jt + iv%varbc_tamdar%vtox(ipred,ipred,iphase,iflt) = 1.0 + end do + end if + + ! if (iv%varbc_tamdar%nobs(iphase,iflt) > 0 .and. iv%varbc_tamdar%ifuse(iphase,iflt) <= 0) then + ! do iobs = 1, iv%varbc_tamdar%nobs(iphase,iflt) + ! isn = iv%varbc_tamdar%obs_sn(iobs,iphase,iflt) + ! if (iv%tamdar(isn)%t(1)%qc >= 0) iv%tamdar(isn)%t(1)%qc = fail_varbc_aircraft + ! end do + ! end if + + else + + iv%varbc_tamdar%ifuse(iphase,iflt) = 0 + + if (iv%varbc_tamdar%nobs(iphase,iflt) > 0) then + do iobs = 1, iv%varbc_tamdar%nobs(iphase,iflt) + isn = iv%varbc_tamdar%obs_sn(iobs,iphase,iflt) + if (iv%tamdar(isn)%t(1)%qc >= 0) iv%tamdar(isn)%t(1)%qc = fail_varbc_aircraft + end do + end if + + end if + end do + + else + + iv%varbc_tamdar%ifuse(:,iflt) = 0 + + if (iv%varbc_tamdar%nobs(nphase+1,iflt) > 0) then + do iobs = 1, iv%varbc_tamdar%nobs(nphase+1,iflt) + isn = iv%varbc_tamdar%obs_sn(iobs,nphase+1,iflt) + if (iv%tamdar(isn)%t(1)%qc >= 0) iv%tamdar(isn)%t(1)%qc = fail_varbc_aircraft + end do + end if + + if (rootproc) & + write(unit=varbc_tamdar_unit,fmt='(10X,A,I4,A,I4)') & + 'Measurement from ', iv%varbc_tamdar%tail_id(iflt),' is insufficient: ', & + iv%varbc_tamdar%nobs_sum(nphase+1,iflt) + + end if + end do + + be % cv % size_jt = size_jt + + if (rootproc .and. ANY(iv%varbc_tamdar%nobs_sum(nphase+1,:) >= varbc_tamdar_nobsmin)) & + write(unit=varbc_tamdar_unit,fmt='(/10X,A)') " ID Predictors(npred,nphase) & Obs Count(nphase)" + + write(stringn,'(I3)') iv%varbc_tamdar%npred + stringn = '(10X,I4,3('//trim(ADJUSTL(stringn))//'f8.3),3I5)' + stringn = trim(adjustl(stringn)) + + do iflt = 1, iv%varbc_tamdar%nair + if (rootproc .and. iv%varbc_tamdar%nobs_sum(nphase+1,iflt) >= varbc_tamdar_nobsmin) then + write(unit=varbc_tamdar_unit,fmt=stringn) & + iv%varbc_tamdar%tail_id(iflt), & + iv%varbc_tamdar%pred(1:iv%varbc_tamdar%npred,1:nphase,iflt), & + iv%varbc_tamdar%nobs_sum(1:nphase,iflt) + end if + end do + + deallocate(op,ot,osn,mach,tdiff) + + if (trace_use) call da_trace_exit("da_varbc_tamdar_pred") + + end subroutine da_varbc_tamdar_pred diff --git a/var/da/da_varbc_tamdar/da_varbc_tamdar_tl.inc b/var/da/da_varbc_tamdar/da_varbc_tamdar_tl.inc new file mode 100755 index 0000000000..ba39a75b48 --- /dev/null +++ b/var/da/da_varbc_tamdar/da_varbc_tamdar_tl.inc @@ -0,0 +1,56 @@ + subroutine da_varbc_tamdar_tl (cv_size, cv, iv, y) + + !---------------------------------------------------! + ! METHOD: delta_d = y - delta_bc ! + ! y = H (delta_x) ! + ! delta_bc = SUM( delta_beta_i Pred_i ) ! + !---------------------------------------------------! + + implicit none + + integer, intent(in) :: cv_size + real, intent(in) :: cv(1:cv_size) + type (iv_type), intent(in) :: iv + type (y_type), intent(inout) :: y ! y = h (xa) + + integer :: isn,ivar,iflt,ipred,iobs,iphase + real :: delta_bc + + real, allocatable :: varbc_param_tl(:) + + + if (trace_use) call da_trace_entry("da_varbc_tamdar_tl") + + allocate( varbc_param_tl(iv%varbc_tamdar%npred) ) + + do iflt = 1, iv%varbc_tamdar%nair + do iphase = 1, iv%varbc_tamdar%nphase + + if (iv%varbc_tamdar%nobs(iphase,iflt) > 0 .and. iv%varbc_tamdar%ifuse(iphase,iflt) > 0) then + + varbc_param_tl(:) = 0.0 + + do ipred = 1, iv%varbc_tamdar%npred + varbc_param_tl(ipred) = & + SUM( cv(iv%varbc_tamdar%index(ipred,iphase,iflt)) * & + iv%varbc_tamdar%vtox(ipred,1:iv%varbc_tamdar%npred,iphase,iflt) ) + end do + + delta_bc = SUM( varbc_param_tl(1:iv%varbc_tamdar%npred)* & + iv%varbc_tamdar%pred(1:iv%varbc_tamdar%npred,iphase,iflt) ) + + do iobs = 1, iv%varbc_tamdar%nobs(iphase,iflt) + isn = iv%varbc_tamdar%obs_sn(iobs,iphase,iflt) + if (iv%tamdar(isn)%t(1)%qc >= 0) & + y%tamdar(isn)%t(1) = y%tamdar(isn)%t(1) + delta_bc + end do + + end if + end do + end do + + deallocate(varbc_param_tl) + + if (trace_use) call da_trace_exit("da_varbc_tamdar_tl") + + end subroutine da_varbc_tamdar_tl diff --git a/var/da/da_varbc_tamdar/da_varbc_tamdar_update.inc b/var/da/da_varbc_tamdar/da_varbc_tamdar_update.inc new file mode 100755 index 0000000000..2dadb36dc7 --- /dev/null +++ b/var/da/da_varbc_tamdar/da_varbc_tamdar_update.inc @@ -0,0 +1,143 @@ + subroutine da_varbc_tamdar_update (cv_size, cv, iv) + + !-----------------------------------! + ! Update and write out parameters ! + !-----------------------------------! + + implicit none + + integer, intent(in) :: cv_size + real, intent(in) :: cv(1:cv_size) + type (iv_type), intent(inout) :: iv + + integer :: i,ipred,iflt,iphase,npred + integer :: iunit,iost + character(len=filename_len) :: filename + character(len=99) :: fmt_param + character(len=3) :: cphz(3) + character(len=30) :: stringn + + character(len=5), allocatable :: tail_id(:) + real, allocatable :: sum_tl(:),varbc_param_tl(:,:,:),param(:,:,:) + integer, allocatable :: ifuse(:,:) + + + if (trace_use) call da_trace_entry("da_varbc_tamdar_update") + + if (rootproc) write(unit=varbc_tamdar_unit,fmt='(//5X,A)') 'Updating parameters' + + npred = iv%varbc_tamdar%npred + cphz = (/'des','asc','cru'/) + + write(stringn,'(I3)') iv%varbc_tamdar%npred + stringn = '(9X,I5,2X,A,2X,A,'//trim(ADJUSTL(stringn))//'f9.3)' + stringn = trim(adjustl(stringn)) + + allocate( sum_tl(npred) ) + allocate( varbc_param_tl(npred,iv%varbc_tamdar%nphase,iv%varbc_tamdar%nair) ) + + write(unit=varbc_tamdar_unit,fmt='(/10X,A)') & + " ID Phase Parameters updated (npred)" + + varbc_param_tl(:,:,:) = 0.0 + + do iflt = 1, iv%varbc_tamdar%nair + do iphase = 1, iv%varbc_tamdar%nphase + if (iv%varbc_tamdar%nobs_sum(iphase,iflt) >= varbc_tamdar_nobsmin) then + + if (iv%varbc_tamdar%nobs(iphase,iflt) > 0 .and. iv%varbc_tamdar%ifuse(iphase,iflt) > 0) then + do ipred = 1, npred + varbc_param_tl(ipred,iphase,iflt) = & + SUM( cv(iv%varbc_tamdar%index(ipred,iphase,iflt)) * & + iv%varbc_tamdar%vtox(ipred,1:npred,iphase,iflt) ) + end do + end if + + do ipred = 1, npred + sum_tl(ipred) = wrf_dm_sum_real( varbc_param_tl(ipred,iphase,iflt) ) + iv%varbc_tamdar%param(ipred,iphase,iflt) = iv%varbc_tamdar%param(ipred,iphase,iflt) + sum_tl(ipred) + end do + + if (rootproc .and. iv%varbc_tamdar%ifuse(iphase,iflt) > 0) then + write(unit=varbc_tamdar_unit,fmt=stringn) & + iv%varbc_tamdar%tail_id(iflt),cphz(iphase),' parameters updated : ',sum_tl(:) + end if + end if + end do + end do + + if (.not. rootproc) then + deallocate(varbc_param_tl) + if (trace_use) call da_trace_exit("da_varbc_tamdar_update") + return + end if + + !-----------------------------------! + ! Write VARBC_TAMDAR.tbl.out file ! + !-----------------------------------! + + allocate( tail_id (iv%varbc_tamdar%nair) ) + allocate( ifuse (iv%varbc_tamdar%nphase,iv%varbc_tamdar%nair) ) + allocate( param (iv%varbc_tamdar%nmaxpred,iv%varbc_tamdar%nphase,iv%varbc_tamdar%nair) ) + + call da_get_unit(iunit) + open(unit=iunit,file='VARBC_TAMDAR.tbl', form='formatted', status='old') + + do i = 1, 13 + read(iunit, *) + end do + do iflt = 1, iv%varbc_tamdar%nair + read(iunit, *) tail_id(iflt), & + ifuse(1:iv%varbc_tamdar%nphase,iflt), & + ( param(1:iv%varbc_tamdar%nmaxpred,iphase,iflt), iphase=1,iv%varbc_tamdar%nphase ) + end do + + close(iunit) + call da_free_unit(iunit) + + + ! write updated parameters + + call da_get_unit(iunit) + filename = 'VARBC_TAMDAR.tbl.out' + + open(unit=iunit,file=trim(adjustl(filename)),iostat=iost,status='replace') + if (iost /= 0) then + message(1)="Cannot open TAMDAR bias correction file "//adjustl(filename) + call da_error(__FILE__,__LINE__,message(1:1)) + end if + + write(iunit, '(A)')'*' + write(iunit, '(A)')' PARAMETER TABLE FOR TAMDAR VARBC' + write(iunit, '(A)')' ' + write(iunit, '(A)')' Table format:' + write(iunit, '(A)')'- ID (1X,A5)' + write(iunit, '(A)')'- Ifuse(ascent/desent/cruise) 1X,3(1X,I1)' + write(iunit, '(A)')'- Parameters 3(5F7.3)' + write(iunit, '(A)')' ' + write(iunit, '(A)')' Preditors(1:5): 1.0, w, (optional: Mach, dT/dt, T)' + write(iunit, '(A)')' ' + write(iunit, '(A)')' Number of aircrafts currently involved in this table' + write(iunit, '(I5)') iv%varbc_tamdar%nair + write(iunit, '(A)')'*' + + do iflt = 1, iv%varbc_tamdar%nair + do iphase = 1, iv%varbc_tamdar%nphase + param(1:npred,iphase,iflt) = iv%varbc_tamdar%param(1:npred,iphase,iflt) + end do + + write(iunit,fmt=trim(adjustl(iv%varbc_tamdar%fmt_param))) & + tail_id(iflt), ifuse(iflt,1:iv%varbc_tamdar%nphase), & + ( param(1:iv%varbc_tamdar%nmaxpred,iphase,iflt), iphase=1,iv%varbc_tamdar%nphase ) + end do + + close(iunit) + call da_free_unit(iunit) + + deallocate(sum_tl, tail_id, ifuse, param, varbc_param_tl) + + if (rootproc) write(unit=varbc_tamdar_unit,fmt='(//5X,A/)') 'VARBC_TAMDAR is done' + + if (trace_use) call da_trace_exit("da_varbc_tamdar_update") + + end subroutine da_varbc_tamdar_update diff --git a/var/run/VARBC_TAMDAR.tbl b/var/run/VARBC_TAMDAR.tbl new file mode 100644 index 0000000000..f08ceba5dd --- /dev/null +++ b/var/run/VARBC_TAMDAR.tbl @@ -0,0 +1,33 @@ +* + PARAMETER TABLE FOR TAMDAR VARBC + + Table format: +- ID (1X,A5) +- Ifuse(ascent/desent/cruise) 1X,3(1X,I1) +- Parameters 3(5F7.3) + + Preditors(1:5): 1.0, w, (optional: Mach, dT/dt, T) + + Number of aircrafts currently involved in this table + 20 +* + 00751 1 1 1 0.041 0.149 0.137 0.239 0.176 0.338 0.273 0.261 0.240 0.223 0.436 0.304 0.270 0.208 0.236 + 00862 1 1 1 0.158 0.191 0.186 0.228 0.191 0.200 0.202 0.204 0.226 0.197 0.014 0.130 0.155 0.201 0.174 + 00871 1 1 1 0.408 0.267 0.287 0.166 0.232 0.523 0.382 0.340 0.336 0.248 0.523 0.337 0.296 0.198 0.246 + 00873 1 1 1 -0.147 0.090 0.052 0.241 0.152 0.093 0.144 0.153 0.209 0.186 0.099 0.166 0.168 0.204 0.186 + 00875 1 1 1 0.384 0.263 0.282 0.121 0.221 0.436 0.331 0.304 0.396 0.241 0.562 0.351 0.304 0.209 0.243 + 00877 1 1 1 0.156 0.175 0.204 0.110 0.199 0.351 0.283 0.279 0.244 0.219 0.139 0.182 0.191 0.210 0.193 + 00878 1 1 1 0.307 0.231 0.245 0.171 0.217 0.299 0.254 0.241 0.230 0.218 0.325 0.250 0.238 0.198 0.223 + 00881 1 1 1 0.083 0.158 0.149 0.191 0.192 0.176 0.185 0.187 0.217 0.195 0.092 0.154 0.168 0.205 0.185 + 00888 1 1 1 -0.019 0.126 0.108 0.216 0.173 -0.039 0.075 0.099 0.190 0.167 0.195 0.205 0.198 0.203 0.194 + 00935 1 1 1 -0.095 0.105 0.061 0.274 0.154 0.206 0.202 0.208 0.121 0.203 0.262 0.223 0.224 0.200 0.209 + 01004 1 1 1 0.362 0.251 0.264 0.191 0.221 0.415 0.318 0.304 0.253 0.224 0.821 0.452 0.392 0.209 0.281 + 01045 1 1 1 0.136 0.185 0.175 0.227 0.194 0.394 0.306 0.292 0.246 0.225 0.350 0.264 0.249 0.199 0.222 + 01071 1 1 1 0.071 0.161 0.147 0.232 0.180 0.238 0.217 0.216 0.224 0.208 0.338 0.254 0.246 0.199 0.221 + 05000 1 1 1 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 + 05001 1 1 1 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 + 05002 1 1 1 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 + 05003 1 1 1 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 + 05004 1 1 1 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 + 05005 1 1 1 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 + 05020 1 1 1 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200