diff --git a/Registry/registry.var b/Registry/registry.var index e3c6c9cfa3..00a39434a7 100644 --- a/Registry/registry.var +++ b/Registry/registry.var @@ -170,6 +170,7 @@ rconfig logical use_radar_rv namelist,wrfvar4 1 .false. - "rad rconfig logical use_radar_rf namelist,wrfvar4 1 .false. - "reflectivity" "" "" rconfig logical use_radar_rqv namelist,wrfvar4 1 .false. - "retrieved water vapor" "" "" rconfig logical use_radar_rhv namelist,wrfvar4 1 .false. - "retr. hydrometeor var" "" "" +rconfig integer radar_rhv_opt namelist,wrfvar4 1 1 - "hydrometeor retrieval option" "2 is for background-dependent scheme" "" rconfig integer radar_rf_opt namelist,wrfvar4 1 1 - "reflectivity DA option" "" "" rconfig real rf_qthres namelist,wrfvar4 1 1e-12 - "mixing ratio threshold" "" "" rconfig real rfmin namelist,wrfvar4 1 0.0 - "min rf for no-rain echo" "" "" diff --git a/var/da/da_radar/da_get_innov_vector_radar.inc b/var/da/da_radar/da_get_innov_vector_radar.inc index 44fcc0c83f..bddcda85c3 100644 --- a/var/da/da_radar/da_get_innov_vector_radar.inc +++ b/var/da/da_radar/da_get_innov_vector_radar.inc @@ -50,7 +50,7 @@ subroutine da_get_innov_vector_radar (it, grid, ob, iv) integer :: irv, irvf integer :: irf, irff - real :: alog_10, czr,czs,czg, zrr,zds,zws,zg,rze + real :: alog_10, czrn, czds, czws, czgr, zrn, zds, zws, zgr, rze real :: ob_radar_rf, bg_rze, bg_rf real :: cwr, cws ! weighting coefficient for mixing ratio @@ -63,10 +63,30 @@ subroutine da_get_innov_vector_radar (it, grid, ob, iv) logical :: echo_non_precip, echo_rf_good + !-------------------------------------------------------- + ! for background-dependent hydrmeteor retrieval scheme + !-------------------------------------------------------- + character(len=filename_len) :: hydro_weight_file + integer :: hydro_weight_unit + integer :: tot_h_index, tot_z_index + integer :: ii, jj, kk, nk + integer :: h_index, z_index + logical :: file_exist, qg_exist + real :: zern_ratio, zews_ratio, zeds_ratio, zegr_ratio ! contributions of each hydrometeor to total reflectivity + real, allocatable :: num_sample(:,:) ! number of samples from the background + real, allocatable :: avg_zern(:,:) ! ze contributed by bin-averaged rainwater + real, allocatable :: avg_zeds(:,:) ! ze contributed by bin-averaged dry snow + real, allocatable :: avg_zews(:,:) ! ze contributed by bin-averaged wet snow + real, allocatable :: avg_zegr(:,:) ! ze contributed by bin-averaged graupel + real, allocatable :: avg_qrn(:,:) ! bin-averaged rainwater + real, allocatable :: avg_qds(:,:) ! bin-averaged dry snow + real, allocatable :: avg_qws(:,:) ! bin-averaged wet snow + real, allocatable :: avg_qgr(:,:) ! bin-averaged graupel + real, allocatable :: ave_rho(:,:) ! bin-averaged air density + !------------------------ ! for jung et al 2008 !------------------------ - real :: qvp,qra,qsn,qgr ! mixing ratio real :: dqra,dqsn,dqgr,dtmk,dqvp real :: dqnr,dqns,dqng @@ -83,10 +103,10 @@ subroutine da_get_innov_vector_radar (it, grid, ob, iv) ! Ze=zv*(ro*v)**1.75 ! Zdb=10*log10(Ze) - zrr = 3.63*1.00e+9 ! rainwater + zrn = 3.63*1.00e+9 ! rainwater zds = 9.80*1.00e+8 ! dry snow zws = 4.26*1.00e+11 ! wet snow - zg = 4.33*1.00e+10 ! grauple + zgr = 4.33*1.00e+10 ! grauple !------------------------ ! for jung et al 2008 @@ -240,6 +260,128 @@ END IF end do end if ! lcl for use_radar_rqv + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! background-dependent hydrometer retrieval scheme ! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + if (use_radar_rhv .and. radar_rhv_opt == 2 ) then + !! allocate variables + tot_h_index = 40 ! from 500 m to 20 km, at an interval of 500 m + tot_z_index = 7 ! from -5 dBZ to 65 dBZ, at an interval of 10 dBZ + allocate (num_sample(tot_h_index,tot_z_index)) + allocate (avg_zern(tot_h_index,tot_z_index)) + allocate (avg_zeds(tot_h_index,tot_z_index)) + allocate (avg_zews(tot_h_index,tot_z_index)) + allocate (avg_zegr(tot_h_index,tot_z_index)) + allocate (avg_qrn(tot_h_index,tot_z_index)) + allocate (avg_qds(tot_h_index,tot_z_index)) + allocate (avg_qws(tot_h_index,tot_z_index)) + allocate (avg_qws(tot_h_index,tot_z_index)) + allocate (avg_qgr(tot_h_index,tot_z_index)) + allocate (ave_rho(tot_h_index,tot_z_index)) + + !! variable initialization + num_sample = 0. + avg_qrn = 0. + avg_qws = 0. + avg_qds = 0. + avg_qgr = 0. + ave_rho = 0. + avg_zern = 0. + avg_zews = 0. + avg_zeds = 0. + avg_zegr = 0. + + !! read historical statistics from hydro_mean.dat if available + hydro_weight_file='hydro_mean.dat' + inquire(file=trim(hydro_weight_file), exist=file_exist) + if (file_exist) then + call da_get_unit(hydro_weight_unit) + open(unit=hydro_weight_unit, file=trim(hydro_weight_file), form='FORMATTED') + read(unit=hydro_weight_unit, fmt='(A)') + do z_index=1,tot_z_index + do h_index=1,tot_h_index + read(hydro_weight_unit, fmt='(2(10x), 4(f19.9,2x))') & + avg_zern(h_index,z_index), avg_zews(h_index,z_index), avg_zeds(h_index,z_index), avg_zegr(h_index,z_index) + end do + end do + close(hydro_weight_unit) + call da_free_unit(hydro_weight_unit) + end if + + !! calculate sum of background states in the current processor + do kk=kds, kde + do jj=jps, jpe + do ii=ips, ipe + !! calculate background reflectivity + call da_radar_rf(grid%xb%qrn(ii,jj,kk), grid%xb%qsn(ii,jj,kk), grid%xb%qgr(ii,jj,kk), & + grid%xb%t(ii,jj,kk)-273.15, grid%xb%rho(ii,jj,kk), bg_rze) + bg_rf = 10.*log10(bg_rze) + !! get the index of reflectvity + z_index = nint(bg_rf/10.)+1 + z_index = max(z_index, 0) ! set to non-precip if below -5 dBZ + z_index = min(z_index, tot_z_index) ! set to 65 dBZ if above + !! get the height index + h_index = nint(grid%xb%h(ii,jj,kk)/500.) + h_index = max(h_index, 1) ! set to 500 m if below + h_index = min(h_index,tot_h_index) ! set to 20 km if above + + !! Sum of the model states of different model levels and reflectivity thresholds + if (z_index .ne. 0 ) then + avg_qrn(h_index,z_index) = avg_qrn(h_index,z_index) + grid%xb%qrn(ii,jj,kk) + if ( grid%xb%t(ii,jj,kk) > 273.15 ) then + avg_qws(h_index,z_index) = avg_qws(h_index,z_index) + grid%xb%qsn(ii,jj,kk) + else + avg_qds(h_index,z_index) = avg_qds(h_index,z_index) + grid%xb%qsn(ii,jj,kk) + end if + avg_qgr(h_index,z_index) = avg_qgr(h_index,z_index) + grid%xb%qgr(ii,jj,kk) + ave_rho(h_index,z_index) = ave_rho(h_index,z_index) + grid%xb%rho(ii,jj,kk) + num_sample(h_index,z_index) = num_sample(h_index,z_index) + 1. + end if + end do ! west-east + end do ! south-north + end do ! bottom-top + + !! sum of all processors and get the averaged background states + do z_index=1,tot_z_index + do h_index=1,tot_h_index + num_sample(h_index,z_index) = wrf_dm_sum_real(num_sample(h_index,z_index)) + if (num_sample(h_index,z_index) .gt. 0) then + ave_rho(h_index,z_index) = wrf_dm_sum_real(ave_rho(h_index,z_index)) / num_sample(h_index,z_index) + avg_qrn(h_index,z_index) = wrf_dm_sum_real(avg_qrn(h_index,z_index)) / num_sample(h_index,z_index) + avg_qws(h_index,z_index) = wrf_dm_sum_real(avg_qws(h_index,z_index)) / num_sample(h_index,z_index) + avg_qds(h_index,z_index) = wrf_dm_sum_real(avg_qds(h_index,z_index)) / num_sample(h_index,z_index) + avg_qgr(h_index,z_index) = wrf_dm_sum_real(avg_qgr(h_index,z_index)) / num_sample(h_index,z_index) + end if + end do + end do + + !! calculate the contributions of each hydrometeor to total reflectivity and save them to hydro_mean.dat.update + hydro_weight_file='hydro_mean.dat.update' + if (rootproc) call da_get_unit(hydro_weight_unit) + if (rootproc) open(unit=hydro_weight_unit, file=trim(hydro_weight_file), form='FORMATTED') + if (rootproc) write(unit=hydro_weight_unit, fmt='(2(a8,2x), 4(a19,2x))') & + "z_index:", "h_index:", "===Rainwater===", "===Wet snow===", "===Dry snow===", "===Graupel===" + do z_index=1,tot_z_index + do h_index=1,tot_h_index + if (num_sample(h_index,z_index) .gt. 10.) then + if (avg_qrn(h_index,z_index) > 0.) & !! rain water + avg_zern(h_index,z_index) = zrn*(ave_rho(h_index,z_index)*avg_qrn(h_index,z_index))**1.75 + if (avg_qws(h_index,z_index) > 0.) & !! wet snow + avg_zews(h_index,z_index) = zws*(ave_rho(h_index,z_index)*avg_qws(h_index,z_index))**1.75 + if (avg_qds(h_index,z_index) > 0.) & !! dry snow + avg_zeds(h_index,z_index) = zds*(ave_rho(h_index,z_index)*avg_qds(h_index,z_index))**1.75 + if (avg_qgr(h_index,z_index) > 0.) & !! graupel + avg_zegr(h_index,z_index) = zgr*(ave_rho(h_index,z_index)*avg_qgr(h_index,z_index))**1.75 + end if + if (rootproc) & + write(unit=hydro_weight_unit, fmt='(2(i8, 2x), 4(f19.9,2x))') z_index, h_index, & + avg_zern(h_index,z_index),avg_zews(h_index,z_index), avg_zeds(h_index,z_index), avg_zegr(h_index,z_index) + end do + end do !bottom-top + if (rootproc) close(hydro_weight_unit) + if (rootproc) call da_get_unit(hydro_weight_unit) + end if !! use_radar_rhv .and. radar_rhv_opt == 2 + do n=iv%info(radar)%n1,iv%info(radar)%n2 if ( use_radar_rf .and. radar_rf_opt==1) then @@ -257,7 +399,6 @@ END IF dxm = iv%info(radar)%dxm(1,n) dym = iv%info(radar)%dym(1,n) - model_ps(n) = dxm *(dym * grid%xb % psac(i, j) + dy * grid%xb%psac(i+1, j)) + & dx *(dym * grid%xb % psac(i,j+1) + dy * grid%xb%psac(i+1,j+1)) + & grid%xb % ptop @@ -393,7 +534,7 @@ END IF end if end if - ! calculate background/model reflectivity + ! Calculate background/model reflectivity if (use_radar_rhv .or. use_radar_rqv) then if ( echo_rf_good ) then call da_radar_rf (model_qrn(k,n),model_qsn(k,n),model_qgr(k,n),model_tc(k,n),model_rho(k,n),bg_rze) @@ -403,8 +544,8 @@ END IF end if end if - ! calculate retrieved hydrometeorological variables - ! Jidong Gao JAS 2013 + ! Calculate retrieved hydrometeorological variables + ! Background-dependent retrieval scheme (Chen et al. 2020 AR; Chen et al. 2021 QJRMS) if (use_radar_rhv) then if ( echo_rf_good ) then @@ -431,48 +572,100 @@ END IF end if !if echo_non_precip end if - ob_radar_rf = min(ob_radar_rf, 55.0) ! if dBZ>55.0, set to 55.0 + ! The original WRFDA hydrometeor retrieval scheme + if (model_tc(k,n) .ge. 5.0) then + czrn = 1.0 + czws = 0.0 + czds = 0.0 + czgr = 0.0 + else if (model_tc(k,n) .ge. 0.0) then + czrn = (model_tc(k,n)+5.0)/10.0 + czws = (1.0-czrn)*zws/(zws+zgr) + czds = 0.0 + czgr = (1.0-czrn)*zgr/(zws+zgr) + else if (model_tc(k,n) .ge. -5.0) then + czrn = (model_tc(k,n)+5.0)/10.0 + czws = 0.0 + czds = (1.0-czrn)*zds/(zds+zgr) + czgr = (1.0-czrn)*zgr/(zds+zgr) + else if (model_tc(k,n) .lt. -5.0) then + czrn = 0.0 + czws = 0.0 + czds = zds/(zds+zgr) + czgr = zgr/(zds+zgr) + end if + + if (radar_rhv_opt == 2) then + ! backgound-dependent reflectivity retrival scheme (Chen et al. 2020, AR; Chen et al. 2021, QJRMS) + !! get the index of reflectvity + z_index = nint(ob_radar_rf/10.+1) + z_index = max(z_index, 0) + z_index = min(z_index, tot_z_index) + !! get the height index + h_index = nint(iv%radar(n)%height(k)/500.) + h_index = max(h_index, 1) + h_index = min(h_index, tot_h_index) + + if (z_index > 0) then + zern_ratio = avg_zern(h_index, z_index) + zews_ratio = avg_zews(h_index, z_index) + zeds_ratio = avg_zeds(h_index, z_index) + zegr_ratio = avg_zegr(h_index, z_index) + ! detect whether rain/snow/graupel exists in certain temperatures. + qg_exist = .true. + ! when T < 273.15K + if (model_tc(k,n) .lt. -5.0) zern_ratio = 0. + if (model_tc(k,n) .lt. 0.0) zews_ratio = 0. + ! when T >= 273.15K + if (model_tc(k,n) .ge. 0.0) then + zeds_ratio = 0. + qg_exist = .false. + do nk = k, iv%info(radar)%levels(n) + if (model_tc(nk,n) .lt. -5.0 .and. ob % radar(n) % rf(nk) .ge. 40.) qg_exist = .true. + end do + end if + if (model_tc(k,n) .ge. 5.0) zews_ratio = 0. + if (.not. qg_exist .or. model_tc(k,n) .ge. 10.0) zegr_ratio = 0. + + ! determine the contributions of each hydrometeor to reflectivity + if ((zern_ratio+zews_ratio+zeds_ratio+zegr_ratio) .gt. 0.) then + czrn = zern_ratio/(zern_ratio+zews_ratio+zeds_ratio+zegr_ratio) + czws = zews_ratio/(zern_ratio+zews_ratio+zeds_ratio+zegr_ratio) + czds = zeds_ratio/(zern_ratio+zews_ratio+zeds_ratio+zegr_ratio) + czgr = zegr_ratio/(zern_ratio+zews_ratio+zeds_ratio+zegr_ratio) + end if + else + ob_radar_rf = -15.0 !! Assign reflectivity below -5.0 dBZ to -15.0 dbZ for suppression + !! No need to tune the weights because of very small impacts + end if + end if + + ! convert dBZ to Z + ob_radar_rf = min(ob_radar_rf, 65.0) ! if dBZ>65.0, set to 65.0 rze = 10.0**(ob_radar_rf*0.1) ! dBZ to Z - if (model_tc(k,n).ge.5.0) then - ! contribution from rain only - ! Z_Qr = 3.63*1.0e9*(rho*Qr)**1.75 - iv % radar(n) % rrno(k) = exp ( log(rze/zrr)/1.75 )/model_rho(k,n) + ! Rainwater mixing ratio + if (czrn .gt. 0.) then + iv % radar(n) % rrno(k) = exp ( log(czrn*rze/zrn)/1.75 )/model_rho(k,n) iv % radar(n) % rrn(k) % qc = 0 + end if - ! rrn and rrno were assigned missing values in read_obs_radar_ascii.inc - ! maximum value check, use the data under threshold 15g/kg - iv % radar(n) % rrno(k) = min(iv%radar(n)%rrno(k), 0.015) - - else if (model_tc(k,n).lt.5.0 .and. model_tc(k,n).gt.-5.0 ) then - ! contribution from rain, snow and graupel - ! Ze = c * Z_Qr + (1-c) * (Z_Qs+Z_Qg) - ! the factor c varies linearly between 0 at t=-5C and 1 at t=5C - czr=(model_tc(k,n)+5)/10.0 - if (model_tc(k,n).le.0.0) then - czs = (1.0-czr)*zds/(zds+zg) ! dry snow - czg = (1.0-czr)*zg/(zds+zg) - iv % radar(n) % rsno(k) = exp ( log(czs*rze/zds)/1.75 )/model_rho(k,n) + ! Snow mixing ratio + if ((czws+czds) .gt. 0.) then + if (model_tc(k,n) .gt. 0.) then + iv % radar(n) % rsno(k) = exp ( log(czws*rze/zws)/1.75 )/model_rho(k,n) + iv % radar(n) % rsn(k) % qc = 0 else - czs = (1.0-czr)*zws/(zws+zg) ! wet snow - czg = (1.0-czr)*zg/(zws+zg) - iv % radar(n) % rsno(k) = exp ( log(czs*rze/zws)/1.75 )/model_rho(k,n) + iv % radar(n) % rsno(k) = exp ( log(czds*rze/zds)/1.75 )/model_rho(k,n) + iv % radar(n) % rsn(k) % qc = 0 end if - iv % radar(n) % rrno(k) = exp ( log(czr*rze/zrr)/1.75 )/model_rho(k,n) - iv % radar(n) % rgro(k) = exp ( log(czg*rze/zg )/1.75 )/model_rho(k,n) - iv % radar(n) % rrn(k) % qc = 0 - iv % radar(n) % rsn(k) % qc = 0 - iv % radar(n) % rgr(k) % qc = 0 + end if - else if (model_tc(k,n).le.-5.0) then - ! contribution from snow and graupel - czs = zds/(zds+zg) - czg = 1.0 - czs - iv % radar(n) % rsno(k) = exp ( log(czs*rze/zds)/1.75 )/model_rho(k,n) - iv % radar(n) % rgro(k) = exp ( log(czg*rze/zg )/1.75 )/model_rho(k,n) - iv % radar(n) % rsn(k) % qc = 0 + ! Graupel mixing ratio + if (czgr .gt. 0.) then + iv % radar(n) % rgro(k) = exp ( log(czgr*rze/zgr)/1.75 )/model_rho(k,n) iv % radar(n) % rgr(k) % qc = 0 - end if ! temp + end if if ( radar_rhv_err_opt == 1 ) then ! rainwater error @@ -643,6 +836,17 @@ END IF deallocate (model_qsn) deallocate (model_qgr) + if ( allocated(num_sample) ) deallocate (num_sample) + if ( allocated(avg_zern) ) deallocate (avg_zern) + if ( allocated(avg_zeds) ) deallocate (avg_zeds) + if ( allocated(avg_zews) ) deallocate (avg_zews) + if ( allocated(avg_zegr) ) deallocate (avg_zegr) + if ( allocated(avg_qrn) ) deallocate (avg_qrn) + if ( allocated(avg_qds) ) deallocate (avg_qds) + if ( allocated(avg_qws) ) deallocate (avg_qws) + if ( allocated(avg_qgr) ) deallocate (avg_qgr) + if ( allocated(ave_rho) ) deallocate (ave_rho) + if ( use_radar_rqv ) then deallocate (model_lcl) deallocate (model_qs_ice) diff --git a/var/da/da_radar/da_radar.f90 b/var/da/da_radar/da_radar.f90 index d971f6f604..507a5df0b8 100644 --- a/var/da/da_radar/da_radar.f90 +++ b/var/da/da_radar/da_radar.f90 @@ -1,20 +1,21 @@ module da_radar use module_domain, only : domain - + use module_dm, only : wrf_dm_sum_real use da_control, only : obs_qc_pointer,max_ob_levels,missing_r, & v_interp_p, v_interp_h, check_max_iv_print, trace_use, & missing, max_error_uv, max_error_t, rootproc, & max_error_p,max_error_q, check_max_iv_unit,check_max_iv, & max_stheight_diff,missing_data,max_error_bq,max_error_slp, & max_error_bt, max_error_buv, radar,fails_error_max, & - use_radar_rv, use_radar_rf,radar_rf_opt,radar_rf_rscl,radar_rv_rscl,rf_noice,rfmin, rf_qthres, use_radar_rhv, use_radar_rqv, & + use_radar_rv, use_radar_rf,radar_rf_opt,radar_rf_rscl,radar_rv_rscl,rf_noice,rfmin, rf_qthres, & + use_radar_rhv, use_radar_rqv, radar_rhv_opt,& below_model_surface,mkz,above_model_lid,& fg_format,fg_format_wrf_arw_regional,fg_format_wrf_nmm_regional,fg_format_wrf_arw_global,& fg_format_kma_global,max_error_rv,max_error_rf, & far_below_model_surface,kms,kme,kts,kte, trace_use_dull,filename_len,& myproc, analysis_date, num_procs , ierr, comm, es_beta, es_gamma, a_ew - use da_control, only : its, ite, jts, jte, ids, ide, jds, jde, ims, ime, jms, jme + use da_control, only : its, ite, jts, jte, ids, ide, jds, jde, ims, ime, jms, jme, ips, ipe, jps, jpe, kds, kde use da_control, only : cloudbase_calc_opt, & radar_non_precip_rf, radar_non_precip_opt, radar_rqv_thresh1, radar_rqv_thresh2, & radar_rqv_rh1, radar_rqv_rh2, radar_non_precip_rh_w, radar_non_precip_rh_i, & diff --git a/var/run/hydro_mean.dat b/var/run/hydro_mean.dat new file mode 100644 index 0000000000..8eb4898a78 --- /dev/null +++ b/var/run/hydro_mean.dat @@ -0,0 +1,281 @@ +z_index: h_index: ===Rainwater=== ===Wet snow=== ===Dry snow=== ===Graupel=== + 1 1 1.110364621 0.000000000 0.000000000 0.000000000 + 1 2 1.191557061 0.000000000 0.000000000 0.000000000 + 1 3 1.161192223 0.000000422 0.000000000 0.000000012 + 1 4 1.184153769 0.000000491 0.000000000 0.000000861 + 1 5 1.215461736 0.000001245 0.000000000 0.000032498 + 1 6 1.249331021 0.000003465 0.000000000 0.000434329 + 1 7 1.188016575 0.004044376 0.000000000 0.006002718 + 1 8 0.658994592 0.258770821 0.000000000 0.023394588 + 1 9 0.076126130 0.493034557 0.041413365 0.089647445 + 1 10 0.018769813 0.001490633 0.602470884 0.246930429 + 1 11 0.012606547 0.000000000 0.682771950 0.218481389 + 1 12 0.006803642 0.000000000 0.774319410 0.166622734 + 1 13 0.002194630 0.000000000 0.973009700 0.103653619 + 1 14 0.000539458 0.000000000 1.057978759 0.062300466 + 1 15 0.000082431 0.000000000 1.124476300 0.031163022 + 1 16 0.000006093 0.000000000 1.181231493 0.016328655 + 1 17 0.000000633 0.000000000 1.186018000 0.008443903 + 1 18 0.000000026 0.000000000 1.184349250 0.005027429 + 1 19 0.000000000 0.000000000 1.204235782 0.002984170 + 1 20 0.000000000 0.000000000 1.247787931 0.001683112 + 1 21 0.000000000 0.000000000 1.182761010 0.001269201 + 1 22 0.000000000 0.000000000 1.147730129 0.001016334 + 1 23 0.000000000 0.000000000 1.081208081 0.001700866 + 1 24 0.000000000 0.000000000 1.058804569 0.002569211 + 1 25 0.000000000 0.000000000 1.069680539 0.004297983 + 1 26 0.000000000 0.000000000 1.105473007 0.007288858 + 1 27 0.000000000 0.000000000 1.123479694 0.012151141 + 1 28 0.000000000 0.000000000 1.077309190 0.021671793 + 1 29 0.000000000 0.000000000 1.052393202 0.048595214 + 1 30 0.000000000 0.000000000 1.075984323 0.070006682 + 1 31 0.000000000 0.000000000 0.931566007 0.133371018 + 1 32 0.000000000 0.000000000 0.747305077 0.325729422 + 1 33 0.000000000 0.000000000 0.644699148 0.487301524 + 1 34 0.000000000 0.000000000 0.450743857 0.445289492 + 1 35 0.000000000 0.000000000 0.000000000 0.000000000 + 1 36 0.000000000 0.000000000 0.000000000 0.000000000 + 1 37 0.000000000 0.000000000 0.000000000 0.000000000 + 1 38 0.000000000 0.000000000 0.000000000 0.000000000 + 1 39 0.000000000 0.000000000 0.000000000 0.000000000 + 1 40 0.000000000 0.000000000 0.000000000 0.000000000 + 2 1 11.485129451 0.000000000 0.000000000 0.000003007 + 2 2 11.396482615 0.000000000 0.000000000 0.000002900 + 2 3 11.185311921 0.000000000 0.000000000 0.000009683 + 2 4 11.235017166 0.000000000 0.000000000 0.000194856 + 2 5 11.906867881 0.000000000 0.000000000 0.002666470 + 2 6 11.398847925 0.000004541 0.000000000 0.022052947 + 2 7 11.083775847 0.033225508 0.000000000 0.146095935 + 2 8 5.601100129 3.143440880 0.000000000 0.453400773 + 2 9 0.385149278 4.427440268 0.842305357 1.350949073 + 2 10 0.074431311 0.005787899 7.263531568 2.738786067 + 2 11 0.061621934 0.000000000 8.502123869 2.235457196 + 2 12 0.025156071 0.000000000 9.505179215 1.626797976 + 2 13 0.010992198 0.000000000 10.101719191 1.078340551 + 2 14 0.003244363 0.000000000 10.231194394 0.741196853 + 2 15 0.000662645 0.000000000 10.552275893 0.472961157 + 2 16 0.000098238 0.000000000 11.064867118 0.251980352 + 2 17 0.000012384 0.000000000 11.666488175 0.141560635 + 2 18 0.000000634 0.000000000 12.120665829 0.083122906 + 2 19 0.000000010 0.000000000 12.226585112 0.057658141 + 2 20 0.000000001 0.000000000 11.835602963 0.052477660 + 2 21 0.000000000 0.000000000 10.691880050 0.053115267 + 2 22 0.000000000 0.000000000 10.217975933 0.045675562 + 2 23 0.000000000 0.000000000 10.638193393 0.080667913 + 2 24 0.000000000 0.000000000 10.270044193 0.097538724 + 2 25 0.000000000 0.000000000 10.790678433 0.129937144 + 2 26 0.000000000 0.000000000 10.138613304 0.195182422 + 2 27 0.000000000 0.000000000 10.388859246 0.390972109 + 2 28 0.000000000 0.000000000 10.026214473 0.655879823 + 2 29 0.000000000 0.000000000 10.448562136 1.164410662 + 2 30 0.000000000 0.000000000 9.255934098 1.542579014 + 2 31 0.000000000 0.000000000 6.739038416 3.212334382 + 2 32 0.000000000 0.000000000 4.253683315 5.936382508 + 2 33 0.000000000 0.000000000 2.620775698 6.022844513 + 2 34 0.000000000 0.000000000 1.516023585 4.217000919 + 2 35 0.000000000 0.000000000 0.000000000 0.000000000 + 2 36 0.000000000 0.000000000 0.000000000 0.000000000 + 2 37 0.000000000 0.000000000 0.000000000 0.000000000 + 2 38 0.000000000 0.000000000 0.000000000 0.000000000 + 2 39 0.000000000 0.000000000 0.000000000 0.000000000 + 2 40 0.000000000 0.000000000 0.000000000 0.000000000 + 3 1 109.932737059 0.000000000 0.000000000 0.007300182 + 3 2 122.092308392 0.000000000 0.000000000 0.004750457 + 3 3 108.299617642 0.000000000 0.000000000 0.012055900 + 3 4 101.792502801 0.000000000 0.000000000 0.030280001 + 3 5 100.160428959 0.000000000 0.000000000 0.116610347 + 3 6 103.856724007 0.000012841 0.000000000 0.481703396 + 3 7 98.090295840 0.378111314 0.000000000 3.373334736 + 3 8 37.170010844 43.881337887 0.000000000 7.231537025 + 3 9 1.986512900 59.245030215 3.788333765 14.542833224 + 3 10 0.257588682 0.079139965 59.148690698 33.963358072 + 3 11 0.142639374 0.000000000 74.379298950 24.212025312 + 3 12 0.062449161 0.000000000 85.695554670 16.296630658 + 3 13 0.044433424 0.000000000 90.834574411 11.697786901 + 3 14 0.013513088 0.000000000 100.614643042 7.783517598 + 3 15 0.002632258 0.000000000 105.313449943 4.882578074 + 3 16 0.000608011 0.000000000 104.875563707 3.515404060 + 3 17 0.000127344 0.000000000 100.332384837 2.697175298 + 3 18 0.000011229 0.000000000 94.444999278 2.283756425 + 3 19 0.000000849 0.000000000 91.343376648 2.281776948 + 3 20 0.000000114 0.000000000 90.858339122 2.950994033 + 3 21 0.000000018 0.000000000 95.499549882 3.737628459 + 3 22 0.000000000 0.000000000 95.559053301 3.206116511 + 3 23 0.000000000 0.000000000 93.723659784 4.629463735 + 3 24 0.000000000 0.000000000 93.300850497 5.858860033 + 3 25 0.000000000 0.000000000 95.178401387 7.444054514 + 3 26 0.000000000 0.000000000 97.461499396 8.145390824 + 3 27 0.000000000 0.000000000 87.133350811 11.663855996 + 3 28 0.000000000 0.000000000 80.834832519 13.460996914 + 3 29 0.000000000 0.000000000 84.814338244 15.834966383 + 3 30 0.000000000 0.000000000 69.391744019 33.927804794 + 3 31 0.000000000 0.000000000 30.083735714 58.070284810 + 3 32 0.000000000 0.000000000 12.257326573 55.287309841 + 3 33 0.000000000 0.000000000 0.000000000 0.000000000 + 3 34 0.000000000 0.000000000 0.000000000 0.000000000 + 3 35 0.000000000 0.000000000 0.000000000 0.000000000 + 3 36 0.000000000 0.000000000 0.000000000 0.000000000 + 3 37 0.000000000 0.000000000 0.000000000 0.000000000 + 3 38 0.000000000 0.000000000 0.000000000 0.000000000 + 3 39 0.000000000 0.000000000 0.000000000 0.000000000 + 3 40 0.000000000 0.000000000 0.000000000 0.000000000 + 4 1 582.593809462 0.000000000 0.000000000 0.275695684 + 4 2 862.795555369 0.000000000 0.000000000 0.445421273 + 4 3 891.383213768 0.000000000 0.000000000 1.101761842 + 4 4 847.098276169 0.000000000 0.000000000 3.643758990 + 4 5 788.052922707 0.000000001 0.000000000 9.016362689 + 4 6 775.932257572 0.000351128 0.000000000 27.869731729 + 4 7 684.552482610 6.937820624 0.000000000 96.504729633 + 4 8 165.141711916 602.918935837 0.000000000 77.354169482 + 4 9 11.110636789 841.094650399 6.179351624 73.141773422 + 4 10 2.592036028 2.737345087 258.937787057 546.620437475 + 4 11 0.844255872 0.000000000 393.879602971 428.287945544 + 4 12 0.317489641 0.000000000 499.466047843 305.905580600 + 4 13 0.191266101 0.000000000 549.092489790 245.526019251 + 4 14 0.082059535 0.000000000 610.522719671 187.512000878 + 4 15 0.053360187 0.000000000 631.848753741 158.688350577 + 4 16 0.034991483 0.000000000 658.186633724 154.982817570 + 4 17 0.017887056 0.000000000 682.103320599 165.294918242 + 4 18 0.003462044 0.000000000 686.448995011 171.743749840 + 4 19 0.000334653 0.000000000 671.701302583 174.346198287 + 4 20 0.000028286 0.000000000 660.096023344 185.927444415 + 4 21 0.000004650 0.000000000 701.677255663 200.345430022 + 4 22 0.000000131 0.000000000 721.448748841 200.200056715 + 4 23 0.000000000 0.000000000 692.662575998 173.064245844 + 4 24 0.000000000 0.000000000 737.784120807 182.548421609 + 4 25 0.000000000 0.000000000 770.023711050 158.227222727 + 4 26 0.000000000 0.000000000 793.757906639 155.100836077 + 4 27 0.000000000 0.000000000 741.519916987 201.064163937 + 4 28 0.000000000 0.000000000 672.930450614 237.650181985 + 4 29 0.000000000 0.000000000 443.019863425 306.677426487 + 4 30 0.000000000 0.000000000 190.427776503 401.203266619 + 4 31 0.000000000 0.000000000 73.085580942 376.938006133 + 4 32 0.000000000 0.000000000 0.000000000 0.000000000 + 4 33 0.000000000 0.000000000 0.000000000 0.000000000 + 4 34 0.000000000 0.000000000 0.000000000 0.000000000 + 4 35 0.000000000 0.000000000 0.000000000 0.000000000 + 4 36 0.000000000 0.000000000 0.000000000 0.000000000 + 4 37 0.000000000 0.000000000 0.000000000 0.000000000 + 4 38 0.000000000 0.000000000 0.000000000 0.000000000 + 4 39 0.000000000 0.000000000 0.000000000 0.000000000 + 4 40 0.000000000 0.000000000 0.000000000 0.000000000 + 5 1 0.000000000 0.000000000 0.000000000 0.000000000 + 5 2 5578.212864738 0.000000000 0.000000000 5.008952656 + 5 3 5372.103162093 0.000000000 0.000000000 15.041257061 + 5 4 5409.303983283 0.000000000 0.000000000 53.293005337 + 5 5 5618.973616023 0.000000000 0.000000000 143.608558988 + 5 6 4879.392758682 0.000059368 0.000000000 577.809533511 + 5 7 4349.351164642 2.596033085 0.000000000 1801.203238199 + 5 8 854.397245908 4577.243654527 0.000000000 921.156400737 + 5 9 75.919974478 9679.349939645 0.367917266 405.301875485 + 5 10 30.335561303 461.116686680 412.412251661 5431.263523222 + 5 11 21.568149623 0.000000000 943.750748633 7496.776413933 + 5 12 6.205921446 0.000000000 1263.739261662 7039.395051868 + 5 13 3.792888765 0.000000000 1407.049350326 6724.745434684 + 5 14 3.038712738 0.000000000 1790.658381642 5919.476026227 + 5 15 2.853361147 0.000000000 2068.818886995 5413.289648478 + 5 16 2.403568418 0.000000000 2216.382731337 5005.192492141 + 5 17 1.565186269 0.000000000 2320.988829920 4802.670353956 + 5 18 0.376286682 0.000000000 2400.454833272 4344.940778358 + 5 19 0.036770184 0.000000000 2588.226752918 3615.891464519 + 5 20 0.002523353 0.000000000 2787.011057487 3022.057443383 + 5 21 0.000221743 0.000000000 2508.866032337 3016.475286259 + 5 22 0.000017548 0.000000000 2586.553484818 3026.145356612 + 5 23 0.000000002 0.000000000 3070.760438617 1818.735698531 + 5 24 0.000000000 0.000000000 2372.313361303 2504.319989388 + 5 25 0.000000000 0.000000000 2670.223119035 1809.783614435 + 5 26 0.000000000 0.000000000 2097.739600125 2207.333859109 + 5 27 0.000000000 0.000000000 1778.791933499 2413.184645816 + 5 28 0.000000000 0.000000000 1151.672803739 2910.940558316 + 5 29 0.000000000 0.000000000 712.314792483 2915.130593453 + 5 30 0.000000000 0.000000000 0.000000000 0.000000000 + 5 31 0.000000000 0.000000000 0.000000000 0.000000000 + 5 32 0.000000000 0.000000000 0.000000000 0.000000000 + 5 33 0.000000000 0.000000000 0.000000000 0.000000000 + 5 34 0.000000000 0.000000000 0.000000000 0.000000000 + 5 35 0.000000000 0.000000000 0.000000000 0.000000000 + 5 36 0.000000000 0.000000000 0.000000000 0.000000000 + 5 37 0.000000000 0.000000000 0.000000000 0.000000000 + 5 38 0.000000000 0.000000000 0.000000000 0.000000000 + 5 39 0.000000000 0.000000000 0.000000000 0.000000000 + 5 40 0.000000000 0.000000000 0.000000000 0.000000000 + 6 1 0.000000000 0.000000000 0.000000000 0.000000000 + 6 2 0.000000000 0.000000000 0.000000000 0.000000000 + 6 3 0.000000000 0.000000000 0.000000000 0.000000000 + 6 4 0.000000000 0.000000000 0.000000000 0.000000000 + 6 5 0.000000000 0.000000000 0.000000000 0.000000000 + 6 6 0.000000000 0.000000000 0.000000000 0.000000000 + 6 7 0.000000000 0.000000000 0.000000000 0.000000000 + 6 8 585.441438425 54783.345060723 0.000000000 435.573114800 + 6 9 274.433475945 58482.456351137 0.000000000 1681.712939361 + 6 10 336.916776730 72160.296616061 7.104634902 7034.338489053 + 6 11 525.273471409 0.000000000 471.183532054 50751.765613477 + 6 12 100.010286359 0.000000000 617.046488536 51050.961892201 + 6 13 34.699959319 0.000000000 870.849970379 46579.609004041 + 6 14 32.881914560 0.000000000 1047.803060470 42096.679451858 + 6 15 34.787300955 0.000000000 998.529025322 39468.443642318 + 6 16 34.937999378 0.000000000 912.658347862 37910.400787624 + 6 17 29.476320004 0.000000000 946.200593796 37480.851970615 + 6 18 13.600879797 0.000000000 830.800168234 36555.339790068 + 6 19 0.000000000 0.000000000 0.000000000 0.000000000 + 6 20 0.000000000 0.000000000 0.000000000 0.000000000 + 6 21 0.000000000 0.000000000 0.000000000 0.000000000 + 6 22 0.000000000 0.000000000 0.000000000 0.000000000 + 6 23 0.000000000 0.000000000 0.000000000 0.000000000 + 6 24 0.000000000 0.000000000 0.000000000 0.000000000 + 6 25 0.000000000 0.000000000 0.000000000 0.000000000 + 6 26 0.000000000 0.000000000 0.000000000 0.000000000 + 6 27 0.000000000 0.000000000 0.000000000 0.000000000 + 6 28 0.000000000 0.000000000 0.000000000 0.000000000 + 6 29 0.000000000 0.000000000 0.000000000 0.000000000 + 6 30 0.000000000 0.000000000 0.000000000 0.000000000 + 6 31 0.000000000 0.000000000 0.000000000 0.000000000 + 6 32 0.000000000 0.000000000 0.000000000 0.000000000 + 6 33 0.000000000 0.000000000 0.000000000 0.000000000 + 6 34 0.000000000 0.000000000 0.000000000 0.000000000 + 6 35 0.000000000 0.000000000 0.000000000 0.000000000 + 6 36 0.000000000 0.000000000 0.000000000 0.000000000 + 6 37 0.000000000 0.000000000 0.000000000 0.000000000 + 6 38 0.000000000 0.000000000 0.000000000 0.000000000 + 6 39 0.000000000 0.000000000 0.000000000 0.000000000 + 6 40 0.000000000 0.000000000 0.000000000 0.000000000 + 7 1 0.000000000 0.000000000 0.000000000 0.000000000 + 7 2 0.000000000 0.000000000 0.000000000 0.000000000 + 7 3 0.000000000 0.000000000 0.000000000 0.000000000 + 7 4 0.000000000 0.000000000 0.000000000 0.000000000 + 7 5 0.000000000 0.000000000 0.000000000 0.000000000 + 7 6 0.000000000 0.000000000 0.000000000 0.000000000 + 7 7 0.000000000 0.000000000 0.000000000 0.000000000 + 7 8 0.000000000 0.000000000 0.000000000 0.000000000 + 7 9 2.583919647 368559.529798930 0.000000000 293.474202984 + 7 10 75.463310355 380388.338450024 0.000000000 3151.421917646 + 7 11 0.000000000 0.000000000 0.000000000 0.000000000 + 7 12 0.000000000 0.000000000 0.000000000 0.000000000 + 7 13 0.000000000 0.000000000 0.000000000 0.000000000 + 7 14 0.000000000 0.000000000 0.000000000 0.000000000 + 7 15 0.000000000 0.000000000 0.000000000 0.000000000 + 7 16 0.000000000 0.000000000 0.000000000 0.000000000 + 7 17 0.000000000 0.000000000 0.000000000 0.000000000 + 7 18 0.000000000 0.000000000 0.000000000 0.000000000 + 7 19 0.000000000 0.000000000 0.000000000 0.000000000 + 7 20 0.000000000 0.000000000 0.000000000 0.000000000 + 7 21 0.000000000 0.000000000 0.000000000 0.000000000 + 7 22 0.000000000 0.000000000 0.000000000 0.000000000 + 7 23 0.000000000 0.000000000 0.000000000 0.000000000 + 7 24 0.000000000 0.000000000 0.000000000 0.000000000 + 7 25 0.000000000 0.000000000 0.000000000 0.000000000 + 7 26 0.000000000 0.000000000 0.000000000 0.000000000 + 7 27 0.000000000 0.000000000 0.000000000 0.000000000 + 7 28 0.000000000 0.000000000 0.000000000 0.000000000 + 7 29 0.000000000 0.000000000 0.000000000 0.000000000 + 7 30 0.000000000 0.000000000 0.000000000 0.000000000 + 7 31 0.000000000 0.000000000 0.000000000 0.000000000 + 7 32 0.000000000 0.000000000 0.000000000 0.000000000 + 7 33 0.000000000 0.000000000 0.000000000 0.000000000 + 7 34 0.000000000 0.000000000 0.000000000 0.000000000 + 7 35 0.000000000 0.000000000 0.000000000 0.000000000 + 7 36 0.000000000 0.000000000 0.000000000 0.000000000 + 7 37 0.000000000 0.000000000 0.000000000 0.000000000 + 7 38 0.000000000 0.000000000 0.000000000 0.000000000 + 7 39 0.000000000 0.000000000 0.000000000 0.000000000 + 7 40 0.000000000 0.000000000 0.000000000 0.000000000