SCALE-RM
Data Types | Functions/Subroutines | Variables
scale_letkf Module Reference

module LETKF for Data-Assimilation More...

Data Types

type  obs_da_value
 
type  obs_grid_type
 
type  obs_info
 

Functions/Subroutines

subroutine, public letkf_setup (OBS_IN_NUM, ensemble_comm, ensemble_nprocs, ensemble_myrank, local_comm, local_nprocs, local_myrank, PRC_NUM_X, PRC_NUM_Y, KA, KS, KE, IA, IS, IE, JA, JS, JE, KMAX, IMAX, JMAX, KHALO, IHALO, JHALO, delta_x, delta_y, Zsfc)
 Setup. More...
 
subroutine, public letkf_finalize ()
 
subroutine, public letkf_obs_readfile (OBS_IN_NUM, OBS_IN_FORMAT, OBS_IN_BASENAME, OBS_IN_MASKFILE)
 
subroutine, public letkf_obs_clear (OBS_IN_NUM)
 
subroutine, public letkf_obs_operator (OBS_IN_NUM, OBS_IN_FORMAT, U, V, W, TEMP, PRES, QV, QC, QR, QI, QS, QG, RH, HGT, TOPO, PS, RAIN, U10M, V10M, T2M, Q2M, nobs_extern)
 
subroutine, public letkf_obs_initialize (OBS_IN_NUM, nobs_extern)
 
subroutine, public letkf_system (OBS_IN_NUM, OBS_IN_FORMAT, U, V, W, TEMP, PRES, QV, QC, QR, QI, QS, QG)
 
subroutine, public letkf_param_estimation_system (PEST_PMAX, PEST_VAR0)
 
subroutine, public letkf_add_inflation_setup (addi3d, addi2d)
 
subroutine radar_superobing (na, nr, ne, radlon, radlat, radz, ze, vr, qcflag, attenuation, nlon, nlat, nlev, lon, lat, z, dlon, dlat, dz, missing, input_is_dbz, lon0, lat0, nobs_sp, grid_index, grid_ref, grid_lon_ref, grid_lat_ref, grid_z_ref, grid_count_ref, grid_vr, grid_lon_vr, grid_lat_vr, grid_z_vr, grid_count_vr)
 

Variables

integer, parameter, public nv3d = 11
 
integer, parameter, public nv2d = 0
 
integer, parameter, public nid_obs = 16
 
integer, parameter, public nobtype = 24
 
integer, parameter, public max_obs_info_meta = 3
 
integer, parameter, public n_qc_steps = 2
 
integer, parameter, public i_before_qc = 1
 
integer, parameter, public i_after_qc = 2
 
integer, parameter n_search_incr = 8
 

Detailed Description

module LETKF for Data-Assimilation

Description
Author
Team SCALE imported from the DA system compiled by Data Assimilation Team
NAMELIST
  • PARAM_LETKF
    nametypedefault valuecomment
    LETKF_DEBUG_LOG logical .false.
    LETKF_DETERMINISTIC_RUN logical
    LETKF_MEM_NODES integer 0 Number of nodes used for one member (0: automatically determined)
    SLOT_START integer 1
    SLOT_END integer 1
    SLOT_BASE integer 1
    SLOT_TINTERVAL real(RP) 3600.0_RP unit: [sec]
    DEPARTURE_STAT_RADAR logical .false.
    INFL_MUL real(RP) 1.0_RP > 0: globally constant covariance inflation
    INFL_MUL_MIN real(RP) -1.0_RP minimum inlfation factor (<= 0: not used)
    INFL_MUL_ADAPTIVE logical .false. if true, outout adaptively estimated 3D inlfation field to 'INFL_MUL_OUT_BASENAME' file
    INFL_ADD real(RP) 0.0_RP additive inflation
    INFL_ADD_SHUFFLE logical .false. shuffle the additive inflation members?
    INFL_ADD_Q_RATIO logical .false.
    INFL_ADD_REF_ONLY logical .false.
    RELAX_ALPHA real(RP) 0.0_RP RTPP relaxation parameter
    RELAX_ALPHA_SPREAD real(RP) 0.0_RP RTPS relaxation parameter
    RELAX_TO_INFLATED_PRIOR logical .false. .true. : relaxation to multiplicatively inflated prior
    GROSS_ERROR real(RP) 5.0_RP
    GROSS_ERROR_RAIN real(RP) -1.0_RP 0: same as GROSS_ERROR
    GROSS_ERROR_RADAR_REF real(RP) -1.0_RP 0: same as GROSS_ERROR
    GROSS_ERROR_RADAR_VR real(RP) -1.0_RP 0: same as GROSS_ERROR
    GROSS_ERROR_RADAR_PRH real(RP) -1.0_RP 0: same as GROSS_ERROR
    GROSS_ERROR_TCX real(RP) -1.0_RP debug ! < 0: same as GROSS_ERROR
    GROSS_ERROR_TCY real(RP) -1.0_RP debug ! < 0: same as GROSS_ERROR
    GROSS_ERROR_TCP real(RP) -1.0_RP debug ! < 0: same as GROSS_ERROR
    Q_UPDATE_TOP real(RP) 0.0_RP water vapor and hydrometeors are updated only below this pressure level (Pa)
    Q_SPRD_MAX real(RP) -1.0_RP maximum q (ensemble spread)/(ensemble mean) (only effective when > 0)
    BOUNDARY_BUFFER_WIDTH real(RP) 0.0_RP
    HORI_LOCAL real(RP), dimension(NOBTYPE) (/ 0.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP /) >0: localization length scale (m)
    0: no localization XXX not implemented yet XXX
    0: same as HORI_LOCAL(1)
    VERT_LOCAL real(RP), dimension(NOBTYPE) (/ 0.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP /) >0: localization length scale [ln(p) or m depends on obstype]
    0: no localization
    0: same as VERT_LOCAL(1)
    TIME_LOCAL real(RP), dimension(NOBTYPE) (/ 0.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP /) >0: localization length scale (sec) XXX not implemented yet XXX
    0: no localization
    0: same as TIME_LOCAL(1)
    HORI_LOCAL_RADAR_OBSNOREF real(RP) -1.0_RP 0: same as HORI_LOCAL(22=PHARAD)
    HORI_LOCAL_RADAR_VR real(RP) -1.0_RP 0: same as HORI_LOCAL(22=PHARAD)
    VERT_LOCAL_RADAR_VR real(RP) -1.0_RP 0: same as VERT_LOCAL(22=PHARAD)
    MAX_NOBS_PER_GRID integer, dimension(NOBTYPE) (/ 0, -1, -1, -1, -1, -1, -1, -1, -1, -1,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,-1, -1, -1, -1/) >0: observation number limit
    0: do not limit observation numbers
    0: same as MAX_NOBS_PER_GRID(1)
    MAX_NOBS_PER_GRID_CRITERION integer 1 1: normalized 3D distance (from closest)
    OBS_MIN_SPACING real(RP), dimension(NOBTYPE) (/ 0.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP /) >0: typical minimum spacing of the obsetvation types in the densest observed area (not tuned carefully yet)
    *this is only used for automatically determine OBS_SORT_GRID_SPACING. if using pre-set OBS_SORT_GRID_SPACING, this has no effect.
    =0: same as OBS_MIN_SPACING(1)
    OBS_SORT_GRID_SPACING real(RP), dimension(NOBTYPE) (/ 0.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP,-1.0_RP, -1.0_RP, -1.0_RP, -1.0_RP /) >0: optimal grid spacing for bucket sorting of observations
    0: automatically determined based on HORI_LOCAL, MAX_NOBS_PER_GRID, and OBS_MIN_SPACING
    0: same as OBS_SORT_GRID_SPACING(1)
    USE_RADAR_REF logical .true.
    USE_RADAR_VR logical .true.
    METHOD_REF_CALC integer 3
    USE_TERMINAL_VELOCITY logical .false.
    USE_OBSERR_RADAR_REF logical .false.
    USE_OBSERR_RADAR_VR logical .false.
    RADAR_REF_THRES_DBZ real(RP) 15.0_RP Threshold of rain/no rain
    MIN_RADAR_REF_MEMBER_OBSRAIN integer 1 Minimum rainy ensemble members for assimilating rainy radar obs
    MIN_RADAR_REF_MEMBER_OBSNORAIN integer 1 Minimum rainy ensemble members for assimilating clear-sky radar obs
    MIN_RADAR_REF_DBZ real(RP) 0.0_RP Minimum reflectivity
    MIN_RADAR_REF_DBZ_VR real(RP) 5.0_RP Minimum reflectivity (dBZ) for Doppler velocity observation
    LOW_REF_SHIFT real(RP) 0.0_RP
    RADAR_ZMAX real(RP) 99.e+3_RP Height limit of radar data to be used
    RADAR_ZMIN real(RP) -99.e+3_RP Height limit of radar data to be used
    RADAR_SO_SIZE_HORI real(RP) 1000.0_RP
    RADAR_SO_SIZE_VERT real(RP) 1000.0_RP
    RADAR_MAX_ABS_VR real(RP) 100.0_RP
    USE_METHOD3_REF_MELT logical .false. Use radar operator considering melting (Xue et al. 2009QJRMS)
    RADAR_BIAS_COR_RAIN logical .false. Simple bias correction for radar obs (rain)
    RADAR_BIAS_COR_CLR logical .false. Simple bias correction for radar obs (clear sky)
    RADAR_BIAS_RAIN_CONST_DBZ real(RP) 0.0_RP Simply bias correction for radar obs (rain)
    RADAR_BIAS_CLR_CONST_DBZ real(RP) 0.0_RP Simply bias correction for radar obs (clear sky)
    RADAR_THIN_LETKF_METHOD integer 0 Thinning method
    RADAR_THIN_LETKF_HGRID integer 1 Horizontal thinning level in obs_local
    RADAR_THIN_LETKF_VGRID integer 1 Vertical thinning level in obs_local
    RADAR_THIN_LETKF_HNEAR integer 1
    RADAR_THIN_LETKF_VNEAR integer 1
    RADAR_THIN_HORI integer 1 Thinning horizontal interval (# of grids)
    RADAR_THIN_VERT integer 1 Thinning vertical interval (# of grids)
    OBSERR_U real(RP) 1.0_RP
    OBSERR_V real(RP) 1.0_RP
    OBSERR_T real(RP) 1.0_RP
    OBSERR_Q real(RP) 0.001_RP
    OBSERR_RH real(RP) 10.0_RP
    OBSERR_PS real(RP) 100.0_RP (Pa)
    OBSERR_RADAR_REF real(RP) 5.0_RP
    OBSERR_RADAR_VR real(RP) 3.0_RP
    OBSERR_PQ real(RP) 0.001_RP (kg/m3)
    LETKF_ENTIRE_GRID_SEARCH_X logical .false. Gather all obs to analyze global constant parameters via ETKF (as in Kotsuki et al.)
    LETKF_ENTIRE_GRID_SEARCH_Y logical .false.

History Output
No history output

Function/Subroutine Documentation

◆ letkf_setup()

subroutine, public scale_letkf::letkf_setup ( integer, intent(in)  OBS_IN_NUM,
integer, intent(in)  ensemble_comm,
integer, intent(in)  ensemble_nprocs,
integer, intent(in)  ensemble_myrank,
integer, intent(in)  local_comm,
integer, intent(in)  local_nprocs,
integer, intent(in)  local_myrank,
integer, intent(in)  PRC_NUM_X,
integer, intent(in)  PRC_NUM_Y,
integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  IA,
integer, intent(in)  IS,
integer, intent(in)  IE,
integer, intent(in)  JA,
integer, intent(in)  JS,
integer, intent(in)  JE,
integer, intent(in)  KMAX,
integer, intent(in)  IMAX,
integer, intent(in)  JMAX,
integer, intent(in)  KHALO,
integer, intent(in)  IHALO,
integer, intent(in)  JHALO,
real(rp), intent(in)  delta_x,
real(rp), intent(in)  delta_y,
real(rp), dimension(:,:), intent(in)  Zsfc 
)

Setup.

Definition at line 615 of file scale_letkf.F90.

615  use mpi
616  use scale_prc, only: &
617  prc_abort
618  use scale_prc_cartesc, only: &
619  prc_2drank
620  implicit none
621 
622  integer, intent(in) :: OBS_IN_NUM
623  integer, intent(in) :: ensemble_comm
624  integer, intent(in) :: ensemble_nprocs
625  integer, intent(in) :: ensemble_myrank
626  integer, intent(in) :: local_comm
627  integer, intent(in) :: local_nprocs
628  integer, intent(in) :: local_myrank
629  integer, intent(in) :: PRC_NUM_X
630  integer, intent(in) :: PRC_NUM_Y
631  integer, intent(in) :: KA, KS, KE
632  integer, intent(in) :: IA, IS, IE
633  integer, intent(in) :: JA, JS, JE
634  integer, intent(in) :: KMAX
635  integer, intent(in) :: IMAX
636  integer, intent(in) :: JMAX
637  integer, intent(in) :: KHALO
638  integer, intent(in) :: IHALO
639  integer, intent(in) :: JHALO
640 
641  real(RP), intent(in) :: delta_x
642  real(RP), intent(in) :: delta_y
643  real(RP), intent(in) :: Zsfc(:,:)
644 
645  logical :: LETKF_DETERMINISTIC_RUN
646 
647  integer :: LETKF_MEM_NODES = 0 ! Number of nodes used for one member (0: automatically determined)
648 
649  namelist / param_letkf / &
650  letkf_debug_log, &
651  letkf_deterministic_run, &
652  letkf_mem_nodes, &
653  slot_start, &
654  slot_end, &
655  slot_base, &
656  slot_tinterval, &
657  departure_stat_radar, &
658  infl_mul, &
659  infl_mul_min, &
660  infl_mul_adaptive, &
661  infl_add, &
662  infl_add_shuffle, &
663  infl_add_q_ratio, &
664  infl_add_ref_only, &
665  relax_alpha, &
666  relax_alpha_spread, &
667  relax_to_inflated_prior, &
668  gross_error, &
669  gross_error_rain, &
670  gross_error_radar_ref, &
671  gross_error_radar_vr, &
672  gross_error_radar_prh, &
673  gross_error_tcx, &
674  gross_error_tcy, &
675  gross_error_tcp, &
676  q_update_top, &
677  q_sprd_max, &
678  boundary_buffer_width, &
679  hori_local, &
680  vert_local, &
681  time_local, &
682  hori_local_radar_obsnoref, &
683  hori_local_radar_vr, &
684  vert_local_radar_vr, &
685  max_nobs_per_grid, &
686  max_nobs_per_grid_criterion, &
687  obs_min_spacing, &
688  obs_sort_grid_spacing, &
689  use_radar_ref, &
690  use_radar_vr, &
691  method_ref_calc, &
692  use_terminal_velocity, &
693  use_obserr_radar_ref, &
694  use_obserr_radar_vr, &
695  radar_ref_thres_dbz, &
696  min_radar_ref_member_obsrain, &
697  min_radar_ref_member_obsnorain, &
698  min_radar_ref_dbz, &
699  min_radar_ref_dbz_vr, &
700  low_ref_shift, &
701  radar_zmax, &
702  radar_zmin, &
703  radar_so_size_hori, &
704  radar_so_size_vert, &
705  radar_max_abs_vr, &
706  use_method3_ref_melt, &
707  radar_bias_cor_rain, &
708  radar_bias_cor_clr, &
709  radar_bias_rain_const_dbz, &
710  radar_bias_clr_const_dbz, &
711  radar_thin_letkf_method, &
712  radar_thin_letkf_hgrid, &
713  radar_thin_letkf_vgrid, &
714  radar_thin_letkf_hnear, &
715  radar_thin_letkf_vnear, &
716  radar_thin_hori, &
717  radar_thin_vert, &
718  obserr_u, &
719  obserr_v, &
720  obserr_t, &
721  obserr_q, &
722  obserr_rh, &
723  obserr_ps, &
724  obserr_radar_ref, &
725  obserr_radar_vr, &
726  obserr_pq, &
727  letkf_entire_grid_search_x, &
728  letkf_entire_grid_search_y
729 
730  integer :: n_mem
731  integer :: n_mempn
732 
733  integer :: i, j, k, n
734  integer :: ierr
735  !---------------------------------------------------------------------------
736 
737  log_newline
738  log_info("LETKF_setup",*) 'Setup'
739 
740  if ( rp == sp ) then
741  datatype = mpi_real
742  else if( rp == dp ) then
743  datatype = mpi_double_precision
744  else
745  log_error("obsope_tool",*) 'The precision has not been implemented yet:', rp
746  call prc_abort
747  endif
748 
749  letkf_deterministic_run = .false.
750 
751  !--- read namelist
752  rewind(io_fid_conf)
753  read(io_fid_conf,nml=param_letkf,iostat=ierr)
754  if( ierr < 0 ) then !--- missing
755  log_info("LETKF_setup",*) 'Not found namelist. Default used.'
756  elseif( ierr > 0 ) then !--- fatal error
757  log_error("LETKF_setup",*) 'Not appropriate names in namelist PARAM_LETKF. Check!'
758  call prc_abort
759  endif
760  log_nml(param_letkf)
761 
762  nproc_x = prc_num_x
763  nproc_y = prc_num_y
764  nens = ensemble_nprocs
765  nensobs = ensemble_nprocs
766  nmem = ensemble_nprocs
767  nlon = imax
768  nlat = jmax
769  nlev = kmax
770  nlonh = ia
771  nlath = ja
772  nlevh = ka
773  nlong = imax * prc_num_x
774  nlatg = jmax * prc_num_y
775  nlevall = nlev * nv3d + nv2d
776  xhalo = ihalo
777  yhalo = jhalo
778  zhalo = khalo
779  start_x = is
780  end_x = ie
781  start_y = js
782  end_y = je
783  start_z = ks
784  end_z = ke
785  mmean = nens + 1
786 
787  comm_ens = ensemble_comm
788  nprc_ens = ensemble_nprocs
789  rank_ens = ensemble_myrank
790  comm_lcl = local_comm
791  nprc_lcl = local_nprocs
792  rank_lcl = local_myrank
793 
794  dx = delta_x
795  dy = delta_y
796 
797  if( letkf_deterministic_run ) then ! set deterministic run
798  ens_with_mdet = .true.
799  ! deterministic run is set to the last of ensemble member
800  mmdet = nens
801  mmdetin = nens
802  mmdetobs = nens
803  nmem = nens - 1 ! member size except for deterministic run
804  end if
805 
806  if( letkf_mem_nodes == 0 ) then
807  letkf_mem_nodes = 1 ! (nprocs_m-1) / PPN + 1
808  end if
809  if( letkf_mem_nodes > 1 ) then
810  n_mem = nprc_ens / letkf_mem_nodes
811  n_mempn = 1
812  else
813  n_mem = nprc_ens
814  n_mempn = 1 ! PPN / nprocs_m
815  end if
816  nitmax = ( nprc_ens - 1 ) / ( n_mem * n_mempn ) + 1
817 
818  i = mod( nlon*nlat, nprc_ens )
819  nij1max = ( nlon*nlat - i ) / nprc_ens + 1
820  if( rank_ens < i ) then
821  nij1 = nij1max
822  else
823  nij1 = nij1max - 1
824  end if
825 
826  allocate( obs( obs_in_num ) )
827  allocate( nij1node( nprc_ens ) )
828 
829  do n = 1, nprc_ens
830  if( n-1 < i ) then
831  nij1node(n) = nij1max
832  else
833  nij1node(n) = nij1max - 1
834  end if
835  end do
836 
837  allocate( rig1( nij1 ) )
838  allocate( rjg1( nij1 ) )
839  allocate( topo1( nij1 ) )
840  allocate( hgt1( nij1, nlev ) )
841 
842  allocate( v3dg( nlev, nlon, nlat, nv3d ) )
843  allocate( v2dg( nlon, nlat, nv2d ) )
844 
845  allocate( v3d( nij1, nlev, nv3d ) )
846  allocate( v2d( nij1, nv2d ) )
847 
848  do j = 1, jmax
849  do i = 1, imax
850  v3dg(1,i,j,1) = real( i + prc_2drank(local_myrank,1) * imax + ihalo, kind=rp )
851  v3dg(1,i,j,2) = real( j + prc_2drank(local_myrank,2) * jmax + jhalo, kind=rp )
852  v3dg(1,i,j,3) = zsfc(i+ihalo,j+jhalo)
853  end do
854  end do
855 
856  call scatter_grd_mpi( mod( nens, n_mem*n_mempn ), v3dg, v2dg, v3d, v2d )
857 
858  rig1(:) = v3d(:,1,1)
859  rjg1(:) = v3d(:,1,2)
860  topo1(:) = v3d(:,1,3)
861 
862  call calc_z_grd( nij1, topo1, hgt1 )
863 
864  allocate( var_local( nv3d+nv2d, nid_obs_varlocal ) )
865 
866  if( radar_ref_thres_dbz < min_radar_ref_dbz ) then
867  radar_ref_thres_dbz = min_radar_ref_dbz
868  end if
869 
870  min_radar_ref = 10.0_rp ** ( min_radar_ref_dbz / 10.0_rp )
871  min_radar_ref_vr = 10.0_rp ** ( min_radar_ref_dbz_vr / 10.0_rp )
872  radar_ref_thres = 10.0_rp ** ( radar_ref_thres_dbz / 10.0_rp )
873  radar_bias_rain_const = 10.0_rp ** ( radar_bias_rain_const_dbz / 10.0_rp )
874  radar_bias_clr_const = 10.0_rp ** ( radar_bias_clr_const_dbz / 10.0_rp )
875 
876  return

References scale_precision::dp, scale_io::io_fid_conf, nv2d, nv3d, scale_prc_cartesc::prc_2drank, scale_prc::prc_abort(), scale_precision::rp, and scale_precision::sp.

Referenced by mod_da_driver::da_driver_setup().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ letkf_finalize()

subroutine, public scale_letkf::letkf_finalize

Definition at line 881 of file scale_letkf.F90.

881  implicit none
882 
883  deallocate( rig1 )
884  deallocate( rjg1 )
885  deallocate( topo1 )
886  deallocate( hgt1 )
887  deallocate( v3dg )
888  deallocate( v2dg )
889  deallocate( v3d )
890  deallocate( v2d )
891 
892  deallocate( var_local )
893 
894  deallocate( obs )
895  deallocate( nij1node )
896 
897  return

Referenced by mod_da_driver::da_driver_finalize().

Here is the caller graph for this function:

◆ letkf_obs_readfile()

subroutine, public scale_letkf::letkf_obs_readfile ( integer, intent(in)  OBS_IN_NUM,
character(len=h_long), dimension(:), intent(in)  OBS_IN_FORMAT,
character(len=h_long), dimension(:), intent(in)  OBS_IN_BASENAME,
character(len=h_long), intent(in)  OBS_IN_MASKFILE 
)

Definition at line 906 of file scale_letkf.F90.

906  use scale_prc, only: &
907  prc_abort
908  use scale_time, only: &
909  time_gettimelabel
910  implicit none
911 
912  integer, intent(in) :: OBS_IN_NUM
913 
914  character(len=H_LONG), intent(in) :: OBS_IN_FORMAT(:)
915  character(len=H_LONG), intent(in) :: OBS_IN_BASENAME(:)
916  character(len=H_LONG), intent(in) :: OBS_IN_MASKFILE
917 
918  character(len=H_LONG) :: obsfile
919  character(len=H_LONG) :: timelabel_obsfile
920 
921  logical :: err
922 
923  integer :: n
924  integer :: ierr
925  !---------------------------------------------------------------------------
926 
927  log_progress(*) 'data-assimilation / LETKF / obs / readfile'
928 
929  if( rank_ens == 0 .and. rank_lcl == 0 ) then
930  timelabel_obsfile = '_???????????????????.dat'
931  call time_gettimelabel( timelabel_obsfile(2:20) )
932 
933  do n = 1, obs_in_num
934  obsfile = trim(obs_in_basename(n))//trim(timelabel_obsfile)
935 
936  if (obs_in_format(n) /= 'PAWR_TOSHIBA' .and. &
937  obs_in_format(n) /= 'MP_PAWR_TOSHIBA' ) then
938  inquire( file=obsfile, exist=err )
939  if( .not. err ) then
940  log_info("LETKF_obs_readfile",*) 'Warning: File (',trim(obsfile),') is not found. Skip.'
941  ! skip process
942  obs(n)%nobs = 0
943  call obs_info_allocate(obs(n), extended=.true.)
944  cycle
945  end if
946  end if
947 
948  select case( obs_in_format(n) )
949  case ( 'PREPBUFR' )
950  call get_nobs( obsfile, 8, obs(n)%nobs )
951  call obs_info_allocate( obs(n), extended=.true. )
952  call read_obs( obsfile, obs(n) )
953  case ( 'RADAR' )
954  call get_nobs_radar( obsfile, obs(n)%nobs, obs(n)%meta(1), obs(n)%meta(2), obs(n)%meta(3) )
955  call obs_info_allocate( obs(n), extended=.true. )
956  call read_obs_radar( obsfile, obs(n) )
957  case ( 'PAWR_JRC' )
958  log_error("LETKF_obs_readfile",*) 'Error: This system has not been implemented yet. (OBS_IN_FORMAT(:) = PAWR_JRC)'
959  call prc_abort
960  ! call read_obs_radar_jrc( obsfile, obs(n) )
961  case ( 'PAWR_TOSHIBA' )
962  call read_obs_radar_toshiba_pawr( obs(n), obsfile )
963  case ( 'MP_PAWR_TOSHIBA' )
964  call read_obs_radar_toshiba_mp_pawr( obs(n), obsfile, obs_in_maskfile )
965  case default
966  log_error("LETKF_obs_readfile",*) 'Error: Unsupported observation file format.'
967  call prc_abort
968  end select
969  end do
970  end if
971 
972  do n = 1, obs_in_num
973  ! communicate obs. data to ensemble world in each local domain master at first
974  if( rank_lcl == 0 ) then
975  call mpi_bcast( obs(n)%nobs, 1, mpi_integer, 0, comm_ens, ierr )
976 
977  if( rank_ens /= 0 ) then
978  call obs_info_allocate( obs(n), extended=.true. )
979  end if
980 
981  call mpi_bcast( obs(n)%elm, obs(n)%nobs, mpi_integer, 0, comm_ens, ierr )
982  call mpi_bcast( obs(n)%lon, obs(n)%nobs, datatype, 0, comm_ens, ierr )
983  call mpi_bcast( obs(n)%lat, obs(n)%nobs, datatype, 0, comm_ens, ierr )
984  call mpi_bcast( obs(n)%lev, obs(n)%nobs, datatype, 0, comm_ens, ierr )
985  call mpi_bcast( obs(n)%dat, obs(n)%nobs, datatype, 0, comm_ens, ierr )
986  call mpi_bcast( obs(n)%err, obs(n)%nobs, datatype, 0, comm_ens, ierr )
987  call mpi_bcast( obs(n)%typ, obs(n)%nobs, mpi_integer, 0, comm_ens, ierr )
988  call mpi_bcast( obs(n)%dif, obs(n)%nobs, datatype, 0, comm_ens, ierr )
989  call mpi_bcast( obs(n)%meta, max_obs_info_meta, datatype, 0, comm_ens, ierr )
990  end if
991 
992  ! broadcast obs. data to local domain
993  call mpi_bcast( obs(n)%nobs, 1, mpi_integer, 0, comm_lcl, ierr )
994 
995  if( rank_lcl /= 0 ) then
996  call obs_info_allocate( obs(n), extended=.true. )
997  end if
998 
999  call mpi_bcast( obs(n)%elm, obs(n)%nobs, mpi_integer, 0, comm_lcl, ierr )
1000  call mpi_bcast( obs(n)%lon, obs(n)%nobs, datatype, 0, comm_lcl, ierr )
1001  call mpi_bcast( obs(n)%lat, obs(n)%nobs, datatype, 0, comm_lcl, ierr )
1002  call mpi_bcast( obs(n)%lev, obs(n)%nobs, datatype, 0, comm_lcl, ierr )
1003  call mpi_bcast( obs(n)%dat, obs(n)%nobs, datatype, 0, comm_lcl, ierr )
1004  call mpi_bcast( obs(n)%err, obs(n)%nobs, datatype, 0, comm_lcl, ierr )
1005  call mpi_bcast( obs(n)%typ, obs(n)%nobs, mpi_integer, 0, comm_lcl, ierr )
1006  call mpi_bcast( obs(n)%dif, obs(n)%nobs, datatype, 0, comm_lcl, ierr )
1007  call mpi_bcast( obs(n)%meta, max_obs_info_meta, datatype, 0, comm_lcl, ierr )
1008 
1009  log_info('LETKF_obs_readfile',*) 'observations input:', obs(n)%nobs
1010  end do
1011 
1012  return

References max_obs_info_meta, scale_prc::prc_abort(), and scale_time::time_gettimelabel().

Referenced by mod_da_driver::da_driver_update().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ letkf_obs_clear()

subroutine, public scale_letkf::letkf_obs_clear ( integer, intent(in)  OBS_IN_NUM)

Definition at line 1016 of file scale_letkf.F90.

1016  implicit none
1017 
1018  integer, intent(in) :: OBS_IN_NUM
1019 
1020  integer :: n
1021  !---------------------------------------------------------------------------
1022 
1023  log_progress(*) 'data-assimilation / LETKF / cleaning'
1024 
1025  do n = 1, obs_in_num
1026  call obs_info_deallocate( obs(n) )
1027  end do
1028 
1029  call obs_da_value_deallocate( obsda_sort )
1030 
1031  return

Referenced by mod_da_driver::da_driver_update().

Here is the caller graph for this function:

◆ letkf_obs_operator()

subroutine, public scale_letkf::letkf_obs_operator ( integer, intent(in)  OBS_IN_NUM,
character(len=h_long), dimension(:), intent(in)  OBS_IN_FORMAT,
real(rp), dimension (nlevh,nlonh,nlath), intent(in)  U,
real(rp), dimension (nlevh,nlonh,nlath), intent(in)  V,
real(rp), dimension (nlevh,nlonh,nlath), intent(in)  W,
real(rp), dimension(nlevh,nlonh,nlath), intent(in)  TEMP,
real(rp), dimension(nlevh,nlonh,nlath), intent(in)  PRES,
real(rp), dimension (nlevh,nlonh,nlath), intent(in)  QV,
real(rp), dimension (nlevh,nlonh,nlath), intent(in)  QC,
real(rp), dimension (nlevh,nlonh,nlath), intent(in)  QR,
real(rp), dimension (nlevh,nlonh,nlath), intent(in)  QI,
real(rp), dimension (nlevh,nlonh,nlath), intent(in)  QS,
real(rp), dimension (nlevh,nlonh,nlath), intent(in)  QG,
real(rp), dimension (nlevh,nlonh,nlath), intent(in)  RH,
real(rp), dimension (nlevh,nlonh,nlath), intent(in)  HGT,
real(rp), dimension(nlonh,nlath), intent(in)  TOPO,
real(rp), dimension (nlonh,nlath), intent(in)  PS,
real(rp), dimension(nlonh,nlath), intent(in)  RAIN,
real(rp), dimension(nlonh,nlath), intent(in)  U10M,
real(rp), dimension(nlonh,nlath), intent(in)  V10M,
real(rp), dimension (nlonh,nlath), intent(in)  T2M,
real(rp), dimension (nlonh,nlath), intent(in)  Q2M,
integer, intent(in), optional  nobs_extern 
)

Definition at line 1061 of file scale_letkf.F90.

1061  implicit none
1062 
1063  integer, intent(in) :: OBS_IN_NUM
1064  character(len=H_LONG), intent(in) :: OBS_IN_FORMAT(:)
1065 
1066  real(RP), intent(in) :: U (nlevh,nlonh,nlath)
1067  real(RP), intent(in) :: V (nlevh,nlonh,nlath)
1068  real(RP), intent(in) :: W (nlevh,nlonh,nlath)
1069  real(RP), intent(in) :: TEMP(nlevh,nlonh,nlath)
1070  real(RP), intent(in) :: PRES(nlevh,nlonh,nlath)
1071  real(RP), intent(in) :: QV (nlevh,nlonh,nlath)
1072  real(RP), intent(in) :: QC (nlevh,nlonh,nlath)
1073  real(RP), intent(in) :: QR (nlevh,nlonh,nlath)
1074  real(RP), intent(in) :: QI (nlevh,nlonh,nlath)
1075  real(RP), intent(in) :: QS (nlevh,nlonh,nlath)
1076  real(RP), intent(in) :: QG (nlevh,nlonh,nlath)
1077  real(RP), intent(in) :: RH (nlevh,nlonh,nlath)
1078  real(RP), intent(in) :: HGT (nlevh,nlonh,nlath)
1079 
1080  real(RP), intent(in) :: TOPO(nlonh,nlath)
1081  real(RP), intent(in) :: PS (nlonh,nlath)
1082  real(RP), intent(in) :: RAIN(nlonh,nlath)
1083  real(RP), intent(in) :: U10M(nlonh,nlath)
1084  real(RP), intent(in) :: V10M(nlonh,nlath)
1085  real(RP), intent(in) :: T2M (nlonh,nlath)
1086  real(RP), intent(in) :: Q2M (nlonh,nlath)
1087 
1088  integer, optional, intent(in) :: nobs_extern
1089 
1090  type(obs_da_value) :: obsda_tmp
1091 
1092  integer :: i, j, k
1093  integer :: it, im, iof, islot, ierr
1094  integer :: n, nn, nsub, nmod, n1, n2
1095 
1096  integer :: nobs ! observation number processed in this subroutine
1097  integer :: nobs_all
1098  integer :: nobs_max_per_file
1099  integer :: nobs_max_per_file_sub
1100  integer :: slot_nobsg
1101 
1102  integer :: ip, ibufs
1103  integer, allocatable :: cntr(:), dspr(:)
1104  integer, allocatable :: cnts(:), dsps(:)
1105  integer, allocatable :: bsn(:,:), bsna(:,:), bsnext(:,:)
1106  integer :: islot_time_out, islot_domain_out
1107 
1108  integer, allocatable :: obrank_bufs(:)
1109  real(RP), allocatable :: ri_bufs(:)
1110  real(RP), allocatable :: rj_bufs(:)
1111 
1112  integer, allocatable :: obset_bufs(:)
1113  integer, allocatable :: obidx_bufs(:)
1114 
1115  integer :: slot_id(SLOT_START:SLOT_END)
1116  real(RP) :: slot_lb(SLOT_START:SLOT_END)
1117  real(RP) :: slot_ub(SLOT_START:SLOT_END)
1118 
1119  real(RP), allocatable :: v3dg(:,:,:,:)
1120  real(RP), allocatable :: v2dg(:,:,:)
1121 
1122  real(RP) :: ril, rjl, rk, rkz
1123 
1124  character(len=4) :: nstr
1125  !-------------------------------------------------------------------------------
1126 
1127  log_progress(*) 'data-assimilation / LETKF / obs / operator'
1128 
1129  !-------------------------------------------------------------------------------
1130  ! First scan of all observation data: Compute their horizontal location and time
1131  !-------------------------------------------------------------------------------
1132 
1133  nobs_all = 0
1134  nobs_max_per_file = 0
1135  do iof = 1, obs_in_num
1136  if (obs(iof)%nobs > nobs_max_per_file) then
1137  nobs_max_per_file = obs(iof)%nobs
1138  end if
1139  if (obsda_run(iof)) then
1140  nobs_all = nobs_all + obs(iof)%nobs
1141  end if
1142  end do
1143 
1144  nobs_max_per_file_sub = (nobs_max_per_file - 1) / nprc_lcl + 1
1145  allocate( obrank_bufs(nobs_max_per_file_sub) )
1146  allocate( ri_bufs(nobs_max_per_file_sub) )
1147  allocate( rj_bufs(nobs_max_per_file_sub) )
1148 
1149  allocate( cntr(nprc_lcl) )
1150  allocate( dspr(nprc_lcl) )
1151 
1152  ! Use all processes to compute the basic obsevration information
1153  ! (locations in model grids and the subdomains they belong to)
1154  !-----------------------------------------------------------------------------
1155 
1156  do iof = 1, obs_in_num
1157  if( obs(iof)%nobs > 0 ) then ! Process basic obsevration information for all observations since this information is not saved in obsda files
1158  ! when using separate observation operators; ignore the 'OBSDA_RUN' setting for this section
1159  nsub = obs(iof)%nobs / nprc_lcl
1160  nmod = mod( obs(iof)%nobs, nprc_lcl )
1161  do ip = 1, nmod
1162  cntr(ip) = nsub + 1
1163  end do
1164  do ip = nmod+1, nprc_lcl
1165  cntr(ip) = nsub
1166  end do
1167  dspr(1) = 0
1168  do ip = 2, nprc_lcl
1169  dspr(ip) = dspr(ip-1) + cntr(ip-1)
1170  end do
1171 
1172  obrank_bufs(:) = -1
1173  do ibufs = 1, cntr(rank_lcl+1)
1174  n = dspr(rank_lcl+1) + ibufs
1175 
1176  call phys2ij(obs(iof)%lon(n), obs(iof)%lat(n), ri_bufs(ibufs), rj_bufs(ibufs))
1177  call rij_rank(ri_bufs(ibufs), rj_bufs(ibufs), obrank_bufs(ibufs))
1178  end do
1179 
1180  call mpi_allgatherv( obrank_bufs, cntr(rank_lcl+1), mpi_integer, obs(iof)%rank, cntr, dspr, mpi_integer, comm_lcl, ierr )
1181  call mpi_allgatherv( ri_bufs, cntr(rank_lcl+1), datatype, obs(iof)%ri, cntr, dspr, datatype, comm_lcl, ierr )
1182  call mpi_allgatherv( rj_bufs, cntr(rank_lcl+1), datatype, obs(iof)%rj, cntr, dspr, datatype, comm_lcl, ierr )
1183 
1184  end if
1185  end do
1186 
1187  deallocate(cntr, dspr)
1188  deallocate(obrank_bufs, ri_bufs, rj_bufs)
1189 
1190  ! Bucket sort of observation wrt. time slots and subdomains using the process rank 0
1191  !-----------------------------------------------------------------------------
1192 
1193  islot_time_out = slot_end + 1 ! slot = SLOT_END+1 for observation not in the assimilation time window
1194  islot_domain_out = slot_end + 2 ! slot = SLOT_END+2 for observation outside of the model domain
1195 
1196  allocate( bsn( slot_start :slot_end+2, 0:nprc_lcl-1 ) )
1197  allocate( bsna( slot_start-1:slot_end+2, 0:nprc_lcl-1 ) )
1198 
1199  if (rank_ens == 0) then
1200  allocate ( obset_bufs(nobs_all) )
1201  allocate ( obidx_bufs(nobs_all) )
1202  end if
1203 
1204  if (rank_ens == 0 .and. rank_lcl == 0) then
1205  allocate( bsnext( slot_start:slot_end+2, 0:nprc_lcl-1 ) )
1206  bsn(:,:) = 0
1207  bsna(:,:) = 0
1208  bsnext(:,:) = 0
1209 
1210  do iof = 1, obs_in_num
1211  if (obsda_run(iof) .and. obs(iof)%nobs > 0) then
1212  do n = 1, obs(iof)%nobs
1213  if (obs(iof)%rank(n) == -1) then
1214  ! process the observations outside of the model domain in process rank 0
1215  bsn(islot_domain_out, 0) = bsn(islot_domain_out, 0) + 1
1216  else
1217  islot = ceiling(obs(iof)%dif(n) / slot_tinterval - 0.5d0) + slot_base
1218  if (islot < slot_start .or. islot > slot_end) then
1219  islot = islot_time_out
1220  end if
1221  bsn(islot, obs(iof)%rank(n)) = bsn(islot, obs(iof)%rank(n)) + 1
1222  end if
1223  end do
1224  end if
1225  end do
1226 
1227  do ip = 0, nprc_lcl-1
1228  if (ip > 0) then
1229  bsna(slot_start-1, ip) = bsna(slot_end+2, ip-1)
1230  end if
1231  do islot = slot_start, slot_end+2
1232  bsna(islot, ip) = bsna(islot-1, ip) + bsn(islot, ip)
1233  end do
1234  bsnext(slot_start:slot_end+2, ip) = bsna(slot_start-1:slot_end+1, ip)
1235  end do
1236 
1237  do iof = 1, obs_in_num
1238  if (obsda_run(iof) .and. obs(iof)%nobs > 0) then
1239  do n = 1, obs(iof)%nobs
1240  if (obs(iof)%rank(n) == -1) then
1241  ! process the observations outside of the model domain in process rank 0
1242  bsnext(islot_domain_out, 0) = bsnext(islot_domain_out, 0) + 1
1243  obset_bufs(bsnext(islot_domain_out, 0)) = iof
1244  obidx_bufs(bsnext(islot_domain_out, 0)) = n
1245  else
1246  islot = ceiling(obs(iof)%dif(n) / slot_tinterval - 0.5d0) + slot_base
1247  if (islot < slot_start .or. islot > slot_end) then
1248  islot = islot_time_out
1249  end if
1250  bsnext(islot, obs(iof)%rank(n)) = bsnext(islot, obs(iof)%rank(n)) + 1
1251  obset_bufs(bsnext(islot, obs(iof)%rank(n))) = iof
1252  obidx_bufs(bsnext(islot, obs(iof)%rank(n))) = n
1253  end if
1254  end do
1255  end if
1256  end do
1257 
1258  deallocate( bsnext )
1259 
1260  end if
1261 
1262  ! Broadcast the bucket-sort observation numbers to all processes and print
1263  !-----------------------------------------------------------------------------
1264 
1265  if( rank_lcl == 0 ) then
1266  call mpi_bcast( bsn, (slot_end-slot_start+3)*nprc_lcl, mpi_integer, 0, comm_ens, ierr )
1267  call mpi_bcast( bsna, (slot_end-slot_start+4)*nprc_lcl, mpi_integer, 0, comm_ens, ierr )
1268  end if
1269  call mpi_bcast( bsn, (slot_end-slot_start+3)*nprc_lcl, mpi_integer, 0, comm_lcl, ierr )
1270  call mpi_bcast( bsna, (slot_end-slot_start+4)*nprc_lcl, mpi_integer, 0, comm_lcl, ierr )
1271 
1272  do islot = slot_start, slot_end
1273  slot_id(islot) = islot - slot_start + 1
1274  slot_lb(islot) = (real(islot - slot_base, rp) - 0.5d0) * slot_tinterval
1275  slot_ub(islot) = (real(islot - slot_base, rp) + 0.5d0) * slot_tinterval
1276  end do
1277 
1278  ! Scatter the basic obsevration information to processes group {myrank_e = 0},
1279  ! each of which only gets the data in its own subdomain
1280  !-----------------------------------------------------------------------------
1281 
1282  nobs = bsna( slot_end+2, rank_lcl ) - bsna( slot_start-1, rank_lcl )
1283 
1284  obsda_tmp%nobs = nobs
1285  call obs_da_value_allocate(obsda_tmp, 0)
1286 
1287  if (present(nobs_extern)) then
1288  obsda%nobs = nobs + nobs_extern
1289  else
1290  obsda%nobs = nobs
1291  end if
1292  call obs_da_value_allocate(obsda, nitmax)
1293 
1294  if (rank_ens == 0) then
1295  allocate (cnts(nprc_lcl))
1296  allocate (dsps(nprc_lcl))
1297  do ip = 0, nprc_lcl-1
1298  dsps(ip+1) = bsna(slot_start-1, ip)
1299  cnts(ip+1) = bsna(slot_end+2, ip) - dsps(ip+1)
1300  end do
1301 
1302  call mpi_scatterv(obset_bufs, cnts, dsps, mpi_integer, obsda_tmp%set, cnts(rank_lcl+1), mpi_integer, 0, comm_lcl, ierr)
1303  call mpi_scatterv(obidx_bufs, cnts, dsps, mpi_integer, obsda_tmp%idx, cnts(rank_lcl+1), mpi_integer, 0, comm_lcl, ierr)
1304 
1305  deallocate (cnts, dsps)
1306  deallocate (obset_bufs, obidx_bufs)
1307  end if
1308 
1309  ! Broadcast the basic obsevration information
1310  ! from processes group {myrank_e = 0} to all processes
1311  !-----------------------------------------------------------------------------
1312 
1313  call mpi_bcast(obsda_tmp%set, nobs, mpi_integer, 0, comm_ens, ierr)
1314  call mpi_bcast(obsda_tmp%idx, nobs, mpi_integer, 0, comm_ens, ierr)
1315 
1316  obsda%set(1:nobs) = obsda_tmp%set
1317  obsda%idx(1:nobs) = obsda_tmp%idx
1318 
1319  !-------------------------------------------------------------------------------
1320  ! Second scan of observation data in own subdomain: Compute H(x), QC, ... etc.
1321  !-------------------------------------------------------------------------------
1322 
1323  allocate ( v3dg(nlevh,nlonh,nlath,nv3dd) )
1324  allocate ( v2dg( nlonh,nlath,nv2dd) )
1325 
1326  do j = 1, nlath
1327  do i = 1, nlonh
1328  do k = 1, nlevh
1329  v3dg(k,i,j,iv3dd_u ) = u(k,i,j)
1330  v3dg(k,i,j,iv3dd_v ) = v(k,i,j)
1331  v3dg(k,i,j,iv3dd_w ) = w(k,i,j)
1332  v3dg(k,i,j,iv3dd_t ) = temp(k,i,j)
1333  v3dg(k,i,j,iv3dd_p ) = pres(k,i,j)
1334  v3dg(k,i,j,iv3dd_q ) = qv(k,i,j)
1335  v3dg(k,i,j,iv3dd_qc ) = qc(k,i,j)
1336  v3dg(k,i,j,iv3dd_qr ) = qr(k,i,j)
1337  v3dg(k,i,j,iv3dd_qi ) = qi(k,i,j)
1338  v3dg(k,i,j,iv3dd_qs ) = qs(k,i,j)
1339  v3dg(k,i,j,iv3dd_qg ) = qg(k,i,j)
1340  v3dg(k,i,j,iv3dd_rh ) = rh(k,i,j)
1341  v3dg(k,i,j,iv3dd_hgt) = hgt(k,i,j)
1342  end do
1343  end do
1344  end do
1345  do j = 1, nlath
1346  do i = 1, nlonh
1347  v2dg(i,j,iv2dd_topo) = topo(i,j)
1348  v2dg(i,j,iv2dd_ps ) = ps(i,j)
1349  v2dg(i,j,iv2dd_rain) = rain(i,j)
1350  v2dg(i,j,iv2dd_u10m) = u10m(i,j)
1351  v2dg(i,j,iv2dd_v10m) = v10m(i,j)
1352  v2dg(i,j,iv2dd_t2m ) = t2m(i,j)
1353  v2dg(i,j,iv2dd_q2m ) = q2m(i,j)
1354  end do
1355  end do
1356 
1357  do it = 1, nitmax
1358  if (nobs > 0) then
1359  obsda_tmp%qc(1:nobs) = iqc_undef
1360  end if
1361 
1362  ! Observations not in the assimilation time window
1363  !
1364  n1 = bsna(islot_time_out-1, rank_lcl) - bsna(slot_start-1, rank_lcl) + 1
1365  n2 = bsna(islot_time_out, rank_lcl) - bsna(slot_start-1, rank_lcl)
1366  if (n1 <= n2) then
1367  obsda_tmp%qc(n1:n2) = iqc_time
1368  end if
1369 
1370  ! Observations outside of the model domain
1371  !
1372  n1 = bsna(islot_domain_out-1, rank_lcl) - bsna(slot_start-1, rank_lcl) + 1
1373  n2 = bsna(islot_domain_out, rank_lcl) - bsna(slot_start-1, rank_lcl)
1374  if (n1 <= n2) then
1375  obsda_tmp%qc(n1:n2) = iqc_out_h
1376  end if
1377 
1378  ! Valid observations: loop over time slots
1379  !
1380  do islot = slot_start, slot_end
1381  n1 = bsna(islot-1, rank_lcl) - bsna(slot_start-1, rank_lcl) + 1
1382  n2 = bsna(islot, rank_lcl) - bsna(slot_start-1, rank_lcl)
1383  slot_nobsg = sum(bsn(islot, :))
1384 
1385  if (slot_nobsg <= 0) then
1386  cycle
1387  end if
1388 
1389  !call read_history( filename, it, islot, v3dg, v2dg )
1390 
1391  do nn = n1, n2
1392  iof = obsda_tmp%set(nn)
1393  n = obsda_tmp%idx(nn)
1394 
1395  call rij_g2l(rank_lcl, obs(iof)%ri(n), obs(iof)%rj(n), ril, rjl)
1396 
1397  if (.not. use_obs(obs(iof)%typ(n))) then
1398  obsda_tmp%qc(nn) = iqc_otype
1399  cycle
1400  end if
1401 
1402  select case (obs_in_format(iof))
1403  !=====================================================================
1404  case ( 'PREPBUFR' )
1405  !---------------------------------------------------------------------
1406  call phys2ijk(v3dg(:,:,:,iv3dd_p), obs(iof)%elm(n), ril, rjl, obs(iof)%lev(n), rk, obsda_tmp%qc(nn))
1407  if (obsda_tmp%qc(nn) == iqc_good) then
1408  call trans_xtoy(obs(iof)%elm(n), ril, rjl, rk, &
1409  obs(iof)%lon(n), obs(iof)%lat(n), v3dg, v2dg, obsda_tmp%val(nn), obsda_tmp%qc(nn))
1410  end if
1411  !=====================================================================
1412  case ( 'RADAR', 'PAWR_TOSHIBA', 'MP_PAWR_TOSHIBA', 'PAWR_JRC', 'HIMAWARI8' )
1413  !---------------------------------------------------------------------
1414  if ( obs(iof)%lev(n) > radar_zmax .or. obs(iof)%lev(n) < radar_zmin ) then
1415  obsda_tmp%qc(nn) = iqc_radar_vhi
1416  else
1417  call phys2ijkz(v3dg(:,:,:,iv3dd_hgt), ril, rjl, obs(iof)%lev(n), rkz, obsda_tmp%qc(nn))
1418  end if
1419  if (obsda_tmp%qc(nn) == iqc_good) then
1420  call trans_xtoy_radar(obs(iof)%elm(n), obs(iof)%meta(1), obs(iof)%meta(2), obs(iof)%meta(3), ril, rjl, rkz, &
1421  obs(iof)%lon(n), obs(iof)%lat(n), obs(iof)%lev(n), v3dg, v2dg, obsda_tmp%val(nn), obsda_tmp%qc(nn))
1422  !if (obsda_tmp%qc(nn) == iqc_ref_low) obsda_tmp%qc(nn) = iqc_good ! when process the observation operator, we don't care if reflectivity is too small
1423 
1424  call itpl_3d( v3dg(:,:,:,iv3dd_p), rkz, ril, rjl, obsda_tmp%pm(nn) )
1425  call itpl_3d( v3dg(:,:,:,iv3dd_t), rkz, ril, rjl, obsda_tmp%tm(nn) )
1426  call itpl_3d( v3dg(:,:,:,iv3dd_q), rkz, ril, rjl, obsda_tmp%qv(nn) )
1427 
1428  end if
1429  end select
1430  end do
1431  end do
1432 
1433  !!! (tentative) avoid a bug with GCC 5-8 compiler in L1421 (process to treat iqc_ref_low as iqc_good) !!!
1434  where( obsda_tmp%qc(:) == iqc_ref_low )
1435  obsda_tmp%qc(:) = iqc_good
1436  end where
1437 
1438  ! Prepare variables that will need to be communicated if obsda_return is given
1439  !
1440  call obs_da_value_partial_reduce_iter(obsda, it, 1, nobs, obsda_tmp%val, obsda_tmp%qc, &
1441  obsda_tmp%qv, obsda_tmp%tm, obsda_tmp%pm )
1442  end do
1443 
1444  deallocate ( v3dg, v2dg )
1445  deallocate ( bsn, bsna )
1446 
1447  call obs_da_value_deallocate( obsda_tmp )
1448 
1449  return

References scale_precision::rp.

Referenced by mod_da_driver::da_driver_update().

Here is the caller graph for this function:

◆ letkf_obs_initialize()

subroutine, public scale_letkf::letkf_obs_initialize ( integer, intent(in)  OBS_IN_NUM,
integer, intent(in), optional  nobs_extern 
)

Definition at line 1453 of file scale_letkf.F90.

1453  use scale_prc, only: &
1454  prc_abort
1455  use scale_const, only: &
1456  undef => const_undef, &
1457  tem00 => const_tem00
1458  IMPLICIT NONE
1459 
1460  integer, intent(in) :: OBS_IN_NUM
1461  integer, optional, intent(in) :: nobs_extern
1462 
1463  INTEGER :: n,i,j,ierr,im,iof,iidx
1464 
1465  integer :: n1, n2
1466 
1467  integer :: mem_ref
1468  real(RP) :: qvs, qdry
1469 
1470  integer :: it,ip
1471  integer :: ityp,ielm,ielm_u,ictype
1472  real(RP) :: target_grdspc
1473 
1474  integer :: myp_i,myp_j
1475  integer :: ip_i,ip_j
1476 
1477  integer :: nobs_sub(n_qc_steps),nobs_g(n_qc_steps)
1478 
1479  integer :: nobs_elms(nid_obs)
1480  integer :: nobs_elms_sum(nid_obs)
1481 
1482  integer :: nobs_intern
1483  integer :: nobs_extern_
1484 
1485 
1486  character(len=3) :: use_obs_print
1487  character(4) :: nstr
1488 
1489  integer :: cnts
1490  integer :: cntr(NPRC_LCL)
1491  integer :: dspr(NPRC_LCL)
1492  integer :: nensobs_div, nensobs_mod
1493  integer :: im_obs_1, im_obs_2, nensobs_part
1494 
1495  integer :: ns_ext, ne_ext, ns_bufr, ne_bufr
1496  integer :: ishift, jshift
1497 
1498  type(obs_da_value) :: obsbufs, obsbufr
1499  integer :: imin1,imax1,jmin1,jmax1,imin2,imax2,jmin2,jmax2
1500 
1501  real(RP),allocatable :: tmpelm(:)
1502  integer :: monit_nobs(nid_obs)
1503  real(RP) :: bias(nid_obs)
1504  real(RP) :: rmse(nid_obs)
1505 
1506  type(obs_da_value) :: obsda_ext
1507  logical :: ctype_use(nid_obs,nobtype)
1508  !-------------------------------------------------------------------------------
1509 
1510  log_progress(*) 'data-assimilation / LETKF / obs / initialize'
1511 
1512  if( present(nobs_extern) ) then
1513  nobs_extern_ = nobs_extern
1514  else
1515  nobs_extern_ = 0
1516  endif
1517 
1518  nobs_intern = obsda%nobs - nobs_extern_
1519  if( letkf_debug_log ) then
1520  log_info("LETKF_debug",'(1x,A,I8)') 'Internally processed observations: ', nobs_intern
1521  log_info("LETKF_debug",'(1x,A,I8)') 'Externally processed observations: ', nobs_extern_
1522  log_info("LETKF_debug",'(1x,A,I8)') 'Total observations: ', obsda%nobs
1523  endif
1524 
1525  !-------------------------------------------------------------------------------
1526  ! Read externally processed observations
1527  !-------------------------------------------------------------------------------
1528 
1529  !if( nobs_extern > 0 ) then
1530  ! n1 = nobs_intern + 1
1531  ! n2 = obsda%nobs
1532 
1533  ! obsda_ext%nobs = nobs_extern
1534  ! call obs_da_value_allocate(obsda_ext,0)
1535 
1536  ! do it = 1, nitmax
1537  ! im = myrank_to_mem(it)
1538  ! if ((im >= 1 .and. im <= MEMBER) .or. im == mmdetin) then
1539  ! if (im <= MEMBER) then
1540  ! obsdafile = OBSDA_IN_BASENAME
1541  ! call filename_replace_mem(obsdafile, im)
1542  ! else if (im == mmean) then
1543  ! obsdafile = OBSDA_MEAN_IN_BASENAME
1544  ! else if (im == mmdet) then
1545  ! obsdafile = OBSDA_MDET_IN_BASENAME
1546  ! end if
1547  ! call read_obs_da(trim(obsdafile)//obsda_suffix,obsda_ext,0)
1548 
1549  ! if (OBSDA_OUT) then
1550  ! if (im <= MEMBER) then
1551  ! obsdafile = OBSDA_OUT_BASENAME
1552  ! call filename_replace_mem(obsdafile, im)
1553  ! else if (im == mmean) then
1554  ! obsdafile = OBSDA_MEAN_OUT_BASENAME
1555  ! else if (im == mmdet) then
1556  ! obsdafile = OBSDA_MDET_OUT_BASENAME
1557  ! end if
1558  ! call write_obs_da(trim(obsdafile)//obsda_suffix,obsda_ext,0,append=.true.)
1559  ! end if
1560 
1561  ! ! variables without an ensemble dimension
1562  ! if (it == 1) then
1563  ! obsda%set(n1:n2) = obsda_ext%set
1564  ! obsda%idx(n1:n2) = obsda_ext%idx
1565  ! end if
1566 
1567  ! call obs_da_value_partial_reduce_iter(obsda, it, n1, n2, obsda_ext%val, obsda_ext%qc, &
1568  ! obsda_ext%qv, obsda_ext%tm, obsda_ext%pm )
1569 
1570  ! end if ! [ (im >= 1 .and. im <= MEMBER) .or. im == mmdetin ]
1571  ! end do ! [ it = 1, nitmax ]
1572 
1573  ! call obs_da_value_deallocate(obsda_ext)
1574 
1575  ! ! Broadcast the observation information shared by members (e.g., grid numbers)
1576  ! !---------------------------------------------------------------------------
1577 
1578  ! if (nprocs_e > MEMBER) then
1579  ! call MPI_BCAST(obsda%set(n1:n2), nobs_extern, MPI_INTEGER, 0, MPI_COMM_e, ierr)
1580  ! call MPI_BCAST(obsda%idx(n1:n2), nobs_extern, MPI_INTEGER, 0, MPI_COMM_e, ierr)
1581  ! end if
1582  !end if
1583 
1584  !-------------------------------------------------------------------------------
1585  ! Allreduce externally processed observations
1586  !---------------------------------------------------------------------------
1587 
1588  call obs_da_value_allreduce( obsda )
1589 
1590  !-------------------------------------------------------------------------------
1591  ! Process observations and quality control (QC)
1592  !-------------------------------------------------------------------------------
1593 
1594  ! Pre-process data
1595  !-----------------------------------------------------------------------------
1596 
1597  ctype_use(:,:) = .false.
1598  do iof = 1, obs_in_num
1599  do n = 1, obs(iof)%nobs
1600  select case (obs(iof)%elm(n))
1601  case (id_radar_ref_obs)
1602  if (obs(iof)%dat(n) >= 0.0d0 .and. obs(iof)%dat(n) < 1.0d10) then
1603  if (obs(iof)%dat(n) < min_radar_ref) then
1604  obs(iof)%elm(n) = id_radar_ref_zero_obs
1605  obs(iof)%dat(n) = min_radar_ref_dbz + low_ref_shift
1606  else
1607  obs(iof)%dat(n) = 10.0d0 * log10(obs(iof)%dat(n))
1608  end if
1609  else
1610  obs(iof)%dat(n) = undef
1611  end if
1612  if (use_obserr_radar_ref) then
1613  obs(iof)%err(n) = obserr_radar_ref
1614  end if
1615  case (id_radar_ref_zero_obs)
1616  obs(iof)%dat(n) = min_radar_ref_dbz + low_ref_shift
1617  if (use_obserr_radar_ref) then
1618  obs(iof)%err(n) = obserr_radar_ref
1619  end if
1620  case (id_radar_vr_obs)
1621  if (use_obserr_radar_vr) then
1622  obs(iof)%err(n) = obserr_radar_vr
1623  end if
1624  end select
1625 
1626  ! mark (elm, typ) combinations for which observations exist
1627  ctype_use(uid_obs(obs(iof)%elm(n)), obs(iof)%typ(n)) = .true.
1628  end do
1629  end do
1630 
1631  ! do this outside of the above obs loop, so these (ctype) arrays can be in ascending order
1632  nctype = count(ctype_use)
1633  if (allocated(elm_ctype )) deallocate(elm_ctype )
1634  if (allocated(elm_u_ctype )) deallocate(elm_u_ctype )
1635  if (allocated(typ_ctype )) deallocate(typ_ctype )
1636  if (allocated(hori_loc_ctype)) deallocate(hori_loc_ctype)
1637  if (allocated(vert_loc_ctype)) deallocate(vert_loc_ctype)
1638  allocate (elm_ctype(nctype))
1639  allocate (elm_u_ctype(nctype))
1640  allocate (typ_ctype(nctype))
1641  allocate (hori_loc_ctype(nctype))
1642  allocate (vert_loc_ctype(nctype))
1643  ictype = 0
1644  ctype_elmtyp(:,:) = 0
1645  do ityp = 1, nobtype
1646  do ielm_u = 1, nid_obs
1647  if (ctype_use(ielm_u, ityp)) then
1648  ictype = ictype + 1
1649  ctype_elmtyp(ielm_u, ityp) = ictype
1650 
1651  elm_ctype(ictype) = elem_uid(ielm_u)
1652  elm_u_ctype(ictype) = ielm_u
1653  typ_ctype(ictype) = ityp
1654 
1655  ! horizontal localization
1656  if (elm_ctype(ictype) == id_radar_ref_zero_obs) then
1657  hori_loc_ctype(ictype) = hori_local_radar_obsnoref
1658  else if (elm_ctype(ictype) == id_radar_vr_obs) then
1659  hori_loc_ctype(ictype) = hori_local_radar_vr
1660  else
1661  hori_loc_ctype(ictype) = hori_local(ityp)
1662  end if
1663  ! vertical localization
1664  if (elm_ctype(ictype) == id_radar_vr_obs) then
1665  vert_loc_ctype(ictype) = vert_local_radar_vr
1666  else
1667  vert_loc_ctype(ictype) = vert_local(ityp)
1668  end if
1669  end if ! [ ctype_use(ielm_u, ityp) ]
1670  end do ! [ ielm_u = 1, nid_obs ]
1671  end do ! [ ityp = 1, nobtype ]
1672 
1673  ! Compute perturbation and departure
1674  ! -- gross error check
1675  ! -- QC based on background (radar reflectivity)
1676  ! -- process Himawari-8 data
1677  !-----------------------------------------------------------------------------
1678 
1679  allocate(tmpelm(obsda%nobs))
1680 
1681  do n = 1, obsda%nobs
1682  IF(obsda%qc(n) > 0) cycle
1683 
1684  iof = obsda%set(n)
1685  iidx = obsda%idx(n)
1686 
1687  tmpelm(n) = obs(iof)%elm(iidx)
1688 
1689  !!! ###### RADAR assimilation ######
1690  if (obs(iof)%elm(iidx) == id_radar_ref_obs .or. obs(iof)%elm(iidx) == id_radar_ref_zero_obs) then
1691  if (.not. use_radar_ref) then
1692  obsda%qc(n) = iqc_otype
1693  cycle
1694  end if
1695 
1696  if (obs(iof)%dat(iidx) == undef) then
1697  obsda%qc(n) = iqc_obs_bad
1698  cycle
1699  end if
1700 
1701  !!! obsda%ensval: already converted to dBZ
1702  mem_ref = 0
1703  do i = 1, nmem
1704  if (obsda%ensval(i,n) > radar_ref_thres_dbz+1.0d-6 ) then
1705  mem_ref = mem_ref + 1
1706  end if
1707  end do
1708 
1709  ! Obs: Rain
1710  if (obs(iof)%dat(iidx) > radar_ref_thres_dbz+1.0d-6) then
1711  if (mem_ref < min_radar_ref_member_obsrain) then
1712  obsda%qc(n) = iqc_ref_mem
1713 
1714  if ( .not. radar_pqv ) cycle
1715  ! When RADAR_PQV=True, pseudo qv obs is assimilated even if mem_ref is
1716  ! too small
1717  end if
1718 
1719  else
1720  ! Obs: No rain
1721  if (mem_ref < min_radar_ref_member_obsnorain) then
1722  obsda%qc(n) = iqc_ref_mem
1723  cycle
1724  end if
1725  end if
1726 
1727  end if
1728 
1729  if (obs(iof)%elm(iidx) == id_radar_vr_obs) then
1730  if (.not. use_radar_vr) then
1731  obsda%qc(n) = iqc_otype
1732  cycle
1733  end if
1734  end if
1735  !!! ###### end RADAR assimilation ######
1736 
1737  obsda%val(n) = 0.0_rp
1738  do i = 1, nmem
1739  obsda%val(n) = obsda%val(n) + obsda%ensval(i,n)
1740  end do
1741  obsda%val(n) = obsda%val(n) / real(nmem,kind=rp)
1742 
1743  do i = 1, nmem
1744  obsda%ensval(i,n) = obsda%ensval(i,n) - obsda%val(n) ! Hdx
1745  end do
1746  obsda%val(n) = obs(iof)%dat(iidx) - obsda%val(n) ! y-Hx
1747  if( ens_with_mdet ) then
1748  obsda%ensval(mmdetobs,n) = obs(iof)%dat(iidx) - obsda%ensval(mmdetobs,n) ! y-Hx for deterministic run
1749  end if
1750 
1751  select case (obs(iof)%elm(iidx)) !gross error
1752  case (id_rain_obs)
1753  IF(abs(obsda%val(n)) > gross_error_rain * obs(iof)%err(iidx)) THEN
1754  obsda%qc(n) = iqc_gross_err
1755  END IF
1756  case (id_radar_ref_obs,id_radar_ref_zero_obs)
1757 
1758  if( radar_pqv .and. obsda%val(n) > radar_pqv_omb ) then
1759 
1760  ! pseudo qv
1761  obsda%val(n) = 0.0_rp
1762  do i = 1, nmem
1763  obsda%val(n) = obsda%val(n) + obsda%eqv(i,n)
1764  enddo
1765  obsda%val(n) = obsda%val(n) / real(nmem, kind=rp)
1766 
1767  do i = 1, nmem
1768  obsda%ensval(i,n) = obsda%eqv(i,n) - obsda%val(n) ! Hdx
1769  enddO
1770 
1771  ! Tetens equation es(Pa)
1772  qvs = 611.2d0*exp(17.67d0*(obsda%tm(n)-tem00)/(obsda%tm(n) - tem00 + 243.5d0))
1773 
1774  ! Saturtion mixing ratio
1775  qvs = 0.622d0*qvs / ( obsda%pm(n) - qvs )
1776 
1777  obsda%val(n) = qvs - obsda%val(n) ! y-Hx
1778 
1779  if (ens_with_mdet) then
1780  obsda%ensval(mmdetobs,n) = qvs - obsda%eqv(mmdetobs,n) ! y-Hx for deterministic run
1781  end if
1782 
1783  obsda%tm(n) = -1.0d0
1784 
1785  else
1786  IF(abs(obsda%val(n)) > gross_error_radar_ref * obs(iof)%err(iidx)) THEN
1787  obsda%qc(n) = iqc_gross_err
1788  END IF
1789  endif
1790  case (id_radar_vr_obs)
1791  IF(abs(obsda%val(n)) > gross_error_radar_vr * obs(iof)%err(iidx)) THEN
1792  obsda%qc(n) = iqc_gross_err
1793  END IF
1794  case (id_radar_prh_obs)
1795  IF(abs(obsda%val(n)) > gross_error_radar_prh * obs(iof)%err(iidx)) THEN
1796  obsda%qc(n) = iqc_gross_err
1797  END IF
1798  case default
1799  IF(abs(obsda%val(n)) > gross_error * obs(iof)%err(iidx)) THEN
1800  obsda%qc(n) = iqc_gross_err
1801  END IF
1802  end select
1803 
1804  END DO
1805 
1806  if( letkf_debug_log ) then
1807  log_info("LETKF_debug",'(1x,A,I6,A)') 'OBSERVATIONAL DEPARTURE STATISTICS (IN THIS SUBDOMAIN #', rank_lcl, '):'
1808 
1809  call monit_dep( obsda%nobs, tmpelm, obsda%val, obsda%qc, monit_nobs, bias, rmse )
1810  call monit_print( monit_nobs, bias, rmse )
1811  end if
1812 
1813  deallocate(tmpelm)
1814 
1815  !-------------------------------------------------------------------------------
1816  ! "Bucket sort" of observations of each combined type (with different sorting meshes)
1817  !-------------------------------------------------------------------------------
1818 
1819  if (allocated(obsgrd)) deallocate(obsgrd)
1820  allocate (obsgrd(nctype))
1821 
1822  ! Determine mesh size for bucket sort
1823  !-----------------------------------------------------------------------------
1824 
1825  do ictype = 1, nctype
1826  ityp = typ_ctype(ictype)
1827 
1828  if( obs_sort_grid_spacing(ityp) > 0.0_rp ) then
1829  target_grdspc = obs_sort_grid_spacing(ityp)
1830  else if( max_nobs_per_grid(ityp) > 0 ) then
1831  target_grdspc = 0.1_rp * sqrt( real(max_nobs_per_grid(ityp), kind=rp) ) * obs_min_spacing(ityp) ! need to be tuned
1832  else
1833  target_grdspc = hori_loc_ctype(ictype) * dist_zero_fac / 6.0_rp ! need to be tuned
1834  end if
1835  obsgrd(ictype)%ngrd_i = min( ceiling( dx * real( nlon, kind=rp ) / target_grdspc), nlon )
1836  obsgrd(ictype)%ngrd_j = min( ceiling( dy * real( nlat, kind=rp ) / target_grdspc), nlat )
1837  obsgrd(ictype)%grdspc_i = dx * real( nlon, kind=rp ) / real( obsgrd( ictype )%ngrd_i, kind=rp )
1838  obsgrd(ictype)%grdspc_j = dy * real( nlat, kind=rp ) / real( obsgrd( ictype )%ngrd_j, kind=rp )
1839  if( letkf_entire_grid_search_x ) then
1840  obsgrd(ictype)%ngrdsch_i = ceiling( dx * nlong / obsgrd( ictype )%grdspc_i )
1841  else
1842  obsgrd(ictype)%ngrdsch_i = ceiling( hori_loc_ctype( ictype ) * dist_zero_fac / obsgrd( ictype )%grdspc_i )
1843  endif
1844  if( letkf_entire_grid_search_y ) then
1845  obsgrd(ictype)%ngrdsch_j = ceiling( dy * nlatg / obsgrd( ictype )%grdspc_j )
1846  else
1847  obsgrd(ictype)%ngrdsch_j = ceiling( hori_loc_ctype( ictype ) * dist_zero_fac / obsgrd( ictype )%grdspc_j )
1848  endif
1849  obsgrd(ictype)%ngrdext_i = obsgrd( ictype )%ngrd_i + obsgrd( ictype )%ngrdsch_i * 2
1850  obsgrd(ictype)%ngrdext_j = obsgrd( ictype )%ngrd_j + obsgrd( ictype )%ngrdsch_j * 2
1851 
1852  allocate (obsgrd(ictype)%n ( obsgrd(ictype)%ngrd_i, obsgrd(ictype)%ngrd_j, 0:nprc_lcl-1))
1853  allocate (obsgrd(ictype)%ac(0:obsgrd(ictype)%ngrd_i, obsgrd(ictype)%ngrd_j, 0:nprc_lcl-1))
1854  allocate (obsgrd(ictype)%tot(0:nprc_lcl-1))
1855  allocate (obsgrd(ictype)%n_ext ( obsgrd(ictype)%ngrdext_i, obsgrd(ictype)%ngrdext_j))
1856  allocate (obsgrd(ictype)%ac_ext(0:obsgrd(ictype)%ngrdext_i, obsgrd(ictype)%ngrdext_j))
1857 
1858  obsgrd(ictype)%n (:,:,:) = 0
1859  obsgrd(ictype)%ac(:,:,:) = 0
1860  obsgrd(ictype)%tot(:) = 0
1861  obsgrd(ictype)%n_ext (:,:) = 0
1862  obsgrd(ictype)%ac_ext(:,:) = 0
1863  obsgrd(ictype)%tot_ext = 0
1864  obsgrd(ictype)%tot_sub(:) = 0
1865  obsgrd(ictype)%tot_g(:) = 0
1866 
1867  allocate (obsgrd(ictype)%next(obsgrd(ictype)%ngrd_i, obsgrd(ictype)%ngrd_j))
1868  end do
1869 
1870  ! First scan: count the observation numbers in each mesh (in each subdomian)
1871  !-----------------------------------------------------------------------------
1872 
1873  do n = 1, obsda%nobs
1874  iof = obsda%set(n)
1875  iidx = obsda%idx(n)
1876  ictype = ctype_elmtyp(uid_obs(obs(iof)%elm(iidx)), obs(iof)%typ(iidx))
1877 
1878  if (obsda%qc(n) == iqc_good) then
1879  call ij_obsgrd( ictype, obs(iof)%ri(iidx), obs(iof)%rj(iidx), i, j )
1880  if (i < 1) i = 1 ! Assume the process assignment was correct,
1881  if (i > obsgrd(ictype)%ngrd_i) i = obsgrd(ictype)%ngrd_i ! so this correction is only to remedy the round-off problem.
1882  if (j < 1) j = 1 !
1883  if (j > obsgrd(ictype)%ngrd_j) j = obsgrd(ictype)%ngrd_j !
1884 
1885  obsgrd(ictype)%n(i,j,rank_lcl) = obsgrd(ictype)%n(i,j,rank_lcl) + 1
1886  obsgrd(ictype)%tot_sub(i_after_qc) = obsgrd(ictype)%tot_sub(i_after_qc) + 1 ! only used for diagnostic print (obs number after qc)
1887  end if
1888 
1889  obsgrd(ictype)%tot_sub(i_before_qc) = obsgrd(ictype)%tot_sub(i_before_qc) + 1 ! only used for diagnostic print (obs number before qc)
1890  end do
1891 
1892  ! Compute the accumulated numbers in each mesh
1893  !-----------------------------------------------------------------------------
1894 
1895  do ictype = 1, nctype
1896  if (ictype > 1) then
1897  obsgrd(ictype)%ac(0,1,rank_lcl) = obsgrd(ictype-1)%ac(obsgrd(ictype-1)%ngrd_i,obsgrd(ictype-1)%ngrd_j,rank_lcl)
1898  end if
1899  do j = 1, obsgrd(ictype)%ngrd_j
1900  if (j > 1) then
1901  obsgrd(ictype)%ac(0,j,rank_lcl) = obsgrd(ictype)%ac(obsgrd(ictype)%ngrd_i,j-1,rank_lcl)
1902  end if
1903  do i = 1, obsgrd(ictype)%ngrd_i
1904  obsgrd(ictype)%ac(i,j,rank_lcl) = obsgrd(ictype)%ac(i-1,j,rank_lcl) + obsgrd(ictype)%n(i,j,rank_lcl)
1905  end do
1906  end do
1907  obsgrd(ictype)%next(1:obsgrd(ictype)%ngrd_i,:) = obsgrd(ictype)%ac(0:obsgrd(ictype)%ngrd_i-1,:,rank_lcl)
1908  end do
1909 
1910  ! Second scan: save the indices of bucket-sorted observations in obsda%keys(:)
1911  !-----------------------------------------------------------------------------
1912 
1913  do n = 1, obsda%nobs
1914  if (obsda%qc(n) == iqc_good) then
1915  iof = obsda%set(n)
1916  iidx = obsda%idx(n)
1917  ictype = ctype_elmtyp(uid_obs(obs(iof)%elm(iidx)), obs(iof)%typ(iidx))
1918 
1919  call ij_obsgrd( ictype, obs(iof)%ri(iidx), obs(iof)%rj(iidx), i, j )
1920  if (i < 1) i = 1 ! Assume the process assignment was correct,
1921  if (i > obsgrd(ictype)%ngrd_i) i = obsgrd(ictype)%ngrd_i ! so this correction is only to remedy the round-off problem.
1922  if (j < 1) j = 1 !
1923  if (j > obsgrd(ictype)%ngrd_j) j = obsgrd(ictype)%ngrd_j !
1924 
1925  obsgrd(ictype)%next(i,j) = obsgrd(ictype)%next(i,j) + 1
1926  obsda%key(obsgrd(ictype)%next(i,j)) = n
1927  end if
1928  end do
1929 
1930  ! ALLREDUCE observation number information from subdomains, and compute total numbers
1931  !-----------------------------------------------------------------------------
1932 
1933  nobs_sub(:) = 0
1934  nobs_g(:) = 0
1935  do ictype = 1, nctype
1936  if (nprc_lcl > 1) then
1937  call mpi_allreduce(mpi_in_place, obsgrd(ictype)%n, obsgrd(ictype)%ngrd_i*obsgrd(ictype)%ngrd_j*nprc_lcl, &
1938  mpi_integer, mpi_sum, comm_lcl, ierr)
1939  call mpi_allreduce(mpi_in_place, obsgrd(ictype)%ac(0:obsgrd(ictype)%ngrd_i,:,:), (obsgrd(ictype)%ngrd_i+1)*obsgrd(ictype)%ngrd_j*nprc_lcl, &
1940  mpi_integer, mpi_sum, comm_lcl, ierr)
1941  end if
1942  call mpi_allreduce(obsgrd(ictype)%tot_sub, obsgrd(ictype)%tot_g, n_qc_steps, mpi_integer, mpi_sum, comm_lcl, ierr)
1943 
1944  if (ictype == 1) then
1945  obsgrd(ictype)%tot(:) = obsgrd(ictype)%ac(obsgrd(ictype)%ngrd_i,obsgrd(ictype)%ngrd_j,:)
1946  else
1947  obsgrd(ictype)%tot(:) = obsgrd(ictype)%ac(obsgrd(ictype)%ngrd_i,obsgrd(ictype)%ngrd_j,:) &
1948  - obsgrd(ictype-1)%ac(obsgrd(ictype-1)%ngrd_i,obsgrd(ictype-1)%ngrd_j,:)
1949  end if
1950 
1951  nobs_sub(:) = nobs_sub(:) + obsgrd(ictype)%tot_sub(:)
1952  nobs_g(:) = nobs_g(:) + obsgrd(ictype)%tot_g(:)
1953 
1954  deallocate (obsgrd(ictype)%next)
1955  end do
1956 
1957  nobstotalg = nobs_g(i_after_qc) ! total obs number in the global domain (all types)
1958 
1959  ! Calculate observation numbers in the extended (localization) subdomain,
1960  ! in preparation for communicating obsetvations in the extended subdomain
1961  !-----------------------------------------------------------------------------
1962 
1963  call rank_1d_2d(rank_lcl, myp_i, myp_j)
1964 
1965  do ictype = 1, nctype
1966  imin1 = myp_i*obsgrd(ictype)%ngrd_i+1 - obsgrd(ictype)%ngrdsch_i
1967  imax1 = (myp_i+1)*obsgrd(ictype)%ngrd_i + obsgrd(ictype)%ngrdsch_i
1968  jmin1 = myp_j*obsgrd(ictype)%ngrd_j+1 - obsgrd(ictype)%ngrdsch_j
1969  jmax1 = (myp_j+1)*obsgrd(ictype)%ngrd_j + obsgrd(ictype)%ngrdsch_j
1970 
1971  do ip = 0, nprc_lcl-1
1972  call rank_1d_2d(ip, ip_i, ip_j)
1973  imin2 = max(1, imin1 - ip_i*obsgrd(ictype)%ngrd_i)
1974  imax2 = min(obsgrd(ictype)%ngrd_i, imax1 - ip_i*obsgrd(ictype)%ngrd_i)
1975  jmin2 = max(1, jmin1 - ip_j*obsgrd(ictype)%ngrd_j)
1976  jmax2 = min(obsgrd(ictype)%ngrd_j, jmax1 - ip_j*obsgrd(ictype)%ngrd_j)
1977  if (imin2 > imax2 .or. jmin2 > jmax2) cycle
1978 
1979  ishift = (ip_i - myp_i) * obsgrd(ictype)%ngrd_i + obsgrd(ictype)%ngrdsch_i
1980  jshift = (ip_j - myp_j) * obsgrd(ictype)%ngrd_j + obsgrd(ictype)%ngrdsch_j
1981  obsgrd(ictype)%n_ext(imin2+ishift:imax2+ishift, jmin2+jshift:jmax2+jshift) = obsgrd(ictype)%n(imin2:imax2, jmin2:jmax2, ip)
1982  end do
1983 
1984  if (ictype > 1) then
1985  obsgrd(ictype)%ac_ext(0,1) = obsgrd(ictype-1)%ac_ext(obsgrd(ictype-1)%ngrdext_i,obsgrd(ictype-1)%ngrdext_j)
1986  end if
1987  do j = 1, obsgrd(ictype)%ngrdext_j
1988  if (j > 1) then
1989  obsgrd(ictype)%ac_ext(0,j) = obsgrd(ictype)%ac_ext(obsgrd(ictype)%ngrdext_i,j-1)
1990  end if
1991  do i = 1, obsgrd(ictype)%ngrdext_i
1992  obsgrd(ictype)%ac_ext(i,j) = obsgrd(ictype)%ac_ext(i-1,j) + obsgrd(ictype)%n_ext(i,j)
1993  end do
1994  end do
1995 
1996  if (ictype == 1) then
1997  obsgrd(ictype)%tot_ext = obsgrd(ictype)%ac_ext(obsgrd(ictype)%ngrdext_i,obsgrd(ictype)%ngrdext_j)
1998  else
1999  obsgrd(ictype)%tot_ext = obsgrd(ictype)%ac_ext(obsgrd(ictype)%ngrdext_i,obsgrd(ictype)%ngrdext_j) &
2000  - obsgrd(ictype-1)%ac_ext(obsgrd(ictype-1)%ngrdext_i,obsgrd(ictype-1)%ngrdext_j)
2001  end if
2002  end do
2003 
2004  if (nctype > 0) then
2005  nobstotal = obsgrd(nctype)%ac_ext(obsgrd(nctype)%ngrdext_i,obsgrd(nctype)%ngrdext_j) ! total obs number in the extended subdomain (all types)
2006 
2007  maxnobs_per_ctype = obsgrd(1)%tot_ext
2008  do ictype = 2, nctype
2009  maxnobs_per_ctype = max(maxnobs_per_ctype, obsgrd(ictype)%tot_ext)
2010  end do
2011  else
2012  nobstotal = 0
2013  maxnobs_per_ctype = 0
2014  end if
2015 
2016  ! Construct sorted obsda_sort:
2017  !-----------------------------------------------------------------------------
2018 
2019  ! 1) Copy the observation data in own subdomains to send buffer (obsbufs) with
2020  ! sorted order
2021  !-----------------------------------------------------------------------------
2022  if (nctype > 0) then
2023  cntr(:) = obsgrd(nctype)%ac(obsgrd(nctype)%ngrd_i,obsgrd(nctype)%ngrd_j,:)
2024  cnts = cntr(rank_lcl+1)
2025  dspr(1) = 0
2026  do ip = 2, nprc_lcl
2027  dspr(ip) = dspr(ip-1) + cntr(ip-1)
2028  end do
2029 
2030  nensobs_mod = mod(nensobs, nprc_ens)
2031  nensobs_div = (nensobs - nensobs_mod) / nprc_ens
2032  if (rank_ens < nensobs_mod) then
2033  im_obs_1 = (nensobs_div+1) * rank_ens + 1
2034  im_obs_2 = (nensobs_div+1) * (rank_ens+1)
2035  nensobs_part = nensobs_div + 1
2036  else
2037  im_obs_1 = (nensobs_div+1) * nensobs_mod + nensobs_div * (rank_ens-nensobs_mod) + 1
2038  im_obs_2 = (nensobs_div+1) * nensobs_mod + nensobs_div * (rank_ens-nensobs_mod+1)
2039  nensobs_part = nensobs_div
2040  end if
2041 
2042  obsbufs%nobs = nobs_sub(i_after_qc)
2043  call obs_da_value_allocate(obsbufs, nensobs_part)
2044 
2045  do n = 1, nobs_sub(i_after_qc)
2046  obsbufs%set(n) = obsda%set(obsda%key(n))
2047  obsbufs%idx(n) = obsda%idx(obsda%key(n))
2048  obsbufs%val(n) = obsda%val(obsda%key(n))
2049  obsbufs%tm(n) = obsda%tm(obsda%key(n))
2050  if (nensobs_part > 0) then
2051  obsbufs%ensval(1:nensobs_part,n) = obsda%ensval(im_obs_1:im_obs_2,obsda%key(n))
2052  end if
2053  obsbufs%qc(n) = obsda%qc(obsda%key(n))
2054  end do
2055  end if
2056 
2057  call obs_da_value_deallocate(obsda)
2058 
2059  ! 2) Communicate to get global observations;
2060  ! for variables with an ensemble dimension (ensval),
2061  ! only obtain data from partial members (nensobs_part) to save memory usage
2062  !-----------------------------------------------------------------------------
2063  if (nctype > 0) then
2064  obsbufr%nobs = nobs_g(i_after_qc)
2065  call obs_da_value_allocate(obsbufr, nensobs_part)
2066 
2067  call mpi_allgatherv( obsbufs%set, cnts, mpi_integer, obsbufr%set, cntr, dspr, mpi_integer, comm_lcl, ierr )
2068  call mpi_allgatherv( obsbufs%idx, cnts, mpi_integer, obsbufr%idx, cntr, dspr, mpi_integer, comm_lcl, ierr )
2069  call mpi_allgatherv( obsbufs%val, cnts, datatype, obsbufr%val, cntr, dspr, datatype, comm_lcl, ierr )
2070  call mpi_allgatherv( obsbufs%tm, cnts, datatype, obsbufr%tm, cntr, dspr, datatype, comm_lcl, ierr )
2071  if (nensobs_part > 0) then
2072  call mpi_allgatherv(obsbufs%ensval, cnts*nensobs_part, datatype, obsbufr%ensval, cntr*nensobs_part, dspr*nensobs_part, datatype, comm_lcl, ierr)
2073  end if
2074  call mpi_allgatherv(obsbufs%qc, cnts, mpi_integer, obsbufr%qc, cntr, dspr, mpi_integer, comm_lcl, ierr)
2075 
2076  call obs_da_value_deallocate(obsbufs)
2077  end if
2078 
2079  ! 3) Copy observation data within the extended (localization) subdomains
2080  ! from receive buffer (obsbufr) to obsda_sort; rearrange with sorted order
2081  !-----------------------------------------------------------------------------
2082  obsda_sort%nobs = nobstotal
2083  call obs_da_value_allocate(obsda_sort, nensobs)
2084 
2085  do ip = 0, nprc_lcl-1
2086  call rank_1d_2d(ip, ip_i, ip_j)
2087 
2088  do ictype = 1, nctype
2089  imin1 = myp_i*obsgrd(ictype)%ngrd_i+1 - obsgrd(ictype)%ngrdsch_i
2090  imax1 = (myp_i+1)*obsgrd(ictype)%ngrd_i + obsgrd(ictype)%ngrdsch_i
2091  jmin1 = myp_j*obsgrd(ictype)%ngrd_j+1 - obsgrd(ictype)%ngrdsch_j
2092  jmax1 = (myp_j+1)*obsgrd(ictype)%ngrd_j + obsgrd(ictype)%ngrdsch_j
2093 
2094  imin2 = max(1, imin1 - ip_i*obsgrd(ictype)%ngrd_i)
2095  imax2 = min(obsgrd(ictype)%ngrd_i, imax1 - ip_i*obsgrd(ictype)%ngrd_i)
2096  jmin2 = max(1, jmin1 - ip_j*obsgrd(ictype)%ngrd_j)
2097  jmax2 = min(obsgrd(ictype)%ngrd_j, jmax1 - ip_j*obsgrd(ictype)%ngrd_j)
2098  if (imin2 > imax2 .or. jmin2 > jmax2) cycle
2099 
2100  ishift = (ip_i - myp_i) * obsgrd(ictype)%ngrd_i + obsgrd(ictype)%ngrdsch_i
2101  jshift = (ip_j - myp_j) * obsgrd(ictype)%ngrd_j + obsgrd(ictype)%ngrdsch_j
2102 
2103  do j = jmin2, jmax2
2104  ns_ext = obsgrd(ictype)%ac_ext(imin2+ishift-1,j+jshift) + 1
2105  ne_ext = obsgrd(ictype)%ac_ext(imax2+ishift ,j+jshift)
2106  if (ns_ext > ne_ext) cycle
2107 
2108  ns_bufr = dspr(ip+1) + obsgrd(ictype)%ac(imin2-1,j,ip) + 1
2109  ne_bufr = dspr(ip+1) + obsgrd(ictype)%ac(imax2 ,j,ip)
2110 
2111  obsda_sort%set(ns_ext:ne_ext) = obsbufr%set(ns_bufr:ne_bufr)
2112  obsda_sort%idx(ns_ext:ne_ext) = obsbufr%idx(ns_bufr:ne_bufr)
2113  obsda_sort%val(ns_ext:ne_ext) = obsbufr%val(ns_bufr:ne_bufr)
2114  obsda_sort%tm(ns_ext:ne_ext) = obsbufr%tm(ns_bufr:ne_bufr)
2115  if (nensobs_part > 0) then
2116  obsda_sort%ensval(im_obs_1:im_obs_2,ns_ext:ne_ext) = obsbufr%ensval(1:nensobs_part,ns_bufr:ne_bufr)
2117  end if
2118  obsda_sort%qc(ns_ext:ne_ext) = obsbufr%qc(ns_bufr:ne_bufr)
2119  end do
2120  end do
2121  end do
2122 
2123  ! Save the keys of observations within the subdomain (excluding the localization buffer area)
2124  obsda_sort%nobs_in_key = 0
2125  do ictype = 1, nctype
2126  imin1 = obsgrd(ictype)%ngrdsch_i + 1
2127  imax1 = obsgrd(ictype)%ngrdsch_i + obsgrd(ictype)%ngrd_i
2128  jmin1 = obsgrd(ictype)%ngrdsch_j + 1
2129  jmax1 = obsgrd(ictype)%ngrdsch_j + obsgrd(ictype)%ngrd_j
2130  call obs_choose_ext(ictype, imin1, imax1, jmin1, jmax1, obsda_sort%nobs_in_key, obsda_sort%key)
2131 
2132  deallocate (obsgrd(ictype)%n)
2133  deallocate (obsgrd(ictype)%ac)
2134  deallocate (obsgrd(ictype)%n_ext)
2135  end do
2136 
2137  if (nctype > 0) then
2138  call obs_da_value_deallocate(obsbufr)
2139  end if
2140 
2141  ! 4) For variables with an ensemble dimension (ensval),
2142  ! ALLREDUCE among the ensemble dimension to obtain data of all members
2143  !-----------------------------------------------------------------------------
2144  if (nprc_ens > 1) then
2145  call mpi_allreduce(mpi_in_place, obsda_sort%ensval, nensobs*nobstotal, datatype, mpi_sum, comm_ens, ierr)
2146  end if
2147 
2148  if( letkf_debug_log ) then
2149  log_info("LETKF_debug",'(1x,A,I6,A)') 'OBSERVATION COUNTS (GLOABL AND IN THIS SUBDOMAIN #', rank_lcl, '):'
2150  log_info("LETKF_debug",'(1x,A)') '====================================================================='
2151  log_info("LETKF_debug",'(1x,A)') 'TYPE VAR GLOBAL GLOBAL SUBDOMAIN SUBDOMAIN EXT_SUBDOMAIN'
2152  log_info("LETKF_debug",'(1x,A)') ' before QC after QC before QC after QC after QC'
2153  log_info("LETKF_debug",'(1x,A)') '---------------------------------------------------------------------'
2154  do ictype = 1, nctype
2155  ityp = typ_ctype(ictype)
2156  ielm_u = elm_u_ctype(ictype)
2157  log_info("LETKF_debug",'(1x,A6,1x,A3,1x,4I11,I14)') obtypelist(ityp), obelmlist(ielm_u), &
2158  obsgrd(ictype)%tot_g(i_before_qc), &
2159  obsgrd(ictype)%tot_g(i_after_qc), &
2160  obsgrd(ictype)%tot_sub(i_before_qc), &
2161  obsgrd(ictype)%tot_sub(i_after_qc), &
2162  obsgrd(ictype)%tot_ext
2163  end do
2164  log_info("LETKF_debug",'(1x,A)') '---------------------------------------------------------------------'
2165  log_info("LETKF_debug",'(1x,A,5x,4I11,I14)') 'TOTAL ', nobs_g(i_before_qc), nobs_g(i_after_qc), nobs_sub(i_before_qc), nobs_sub(i_after_qc)
2166  log_info("LETKF_debug",'(1x,A)') '====================================================================='
2167  end if
2168 
2169  return

References scale_const::const_tem00, scale_const::const_undef, i_after_qc, i_before_qc, n_qc_steps, nid_obs, nobtype, scale_prc::prc_abort(), and scale_precision::rp.

Referenced by mod_da_driver::da_driver_update().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ letkf_system()

subroutine, public scale_letkf::letkf_system ( integer, intent(in)  OBS_IN_NUM,
character(len=h_long), dimension(:), intent(in)  OBS_IN_FORMAT,
real(rp), dimension (nlev,nlon,nlat), intent(inout)  U,
real(rp), dimension (nlev,nlon,nlat), intent(inout)  V,
real(rp), dimension (nlev,nlon,nlat), intent(inout)  W,
real(rp), dimension(nlev,nlon,nlat), intent(inout)  TEMP,
real(rp), dimension(nlev,nlon,nlat), intent(inout)  PRES,
real(rp), dimension (nlev,nlon,nlat), intent(inout)  QV,
real(rp), dimension (nlev,nlon,nlat), intent(inout)  QC,
real(rp), dimension (nlev,nlon,nlat), intent(inout)  QR,
real(rp), dimension (nlev,nlon,nlat), intent(inout)  QI,
real(rp), dimension (nlev,nlon,nlat), intent(inout)  QS,
real(rp), dimension (nlev,nlon,nlat), intent(inout)  QG 
)

Definition at line 2189 of file scale_letkf.F90.

2189  use scale_prc, only: &
2190  prc_abort
2191  use scale_const, only: &
2192  undef => const_undef
2193  use scale_statistics, only: &
2194  average => statistics_average
2195  use scale_matrix, only: &
2196  matrix_solver_eigenvalue_decomposition
2197  use scale_random, only: &
2198  knuth_shuffle => random_knuth_shuffle
2199  implicit none
2200 
2201  integer, intent(in) :: OBS_IN_NUM
2202  character(len=H_LONG), intent(in) :: OBS_IN_FORMAT(:)
2203 
2204  real(RP), intent(inout) :: U (nlev,nlon,nlat)
2205  real(RP), intent(inout) :: V (nlev,nlon,nlat)
2206  real(RP), intent(inout) :: W (nlev,nlon,nlat)
2207  real(RP), intent(inout) :: TEMP(nlev,nlon,nlat)
2208  real(RP), intent(inout) :: PRES(nlev,nlon,nlat)
2209  real(RP), intent(inout) :: QV (nlev,nlon,nlat)
2210  real(RP), intent(inout) :: QC (nlev,nlon,nlat)
2211  real(RP), intent(inout) :: QR (nlev,nlon,nlat)
2212  real(RP), intent(inout) :: QI (nlev,nlon,nlat)
2213  real(RP), intent(inout) :: QS (nlev,nlon,nlat)
2214  real(RP), intent(inout) :: QG (nlev,nlon,nlat)
2215 
2216  real(RP) :: gues3d(nij1,nlev,nens+1,nv3d) ! background ensemble
2217  real(RP) :: gues2d(nij1, nens+1,nv2d) ! background ensemble
2218 
2219  real(RP) :: anal3d(nij1,nlev,nens+1,nv3d) ! analysis ensemble
2220  real(RP) :: anal2d(nij1, nens+1,nv2d) ! analysis ensemble
2221 
2222  real(RP) :: addi3d(nij1,nlev,nens+1,nv3d)
2223  real(RP) :: addi2d(nij1, nens+1,nv2d)
2224 
2225  real(RP) :: v3dg(nlev,nlon,nlat,nv3d)
2226  real(RP) :: v2dg( nlon,nlat,nv2d)
2227 
2228  real(RP) :: work3d(nij1,nlev,nv3d)
2229  real(RP) :: work2d(nij1, nv2d)
2230  real(RP), allocatable :: work3da(:,:,:) !GYL
2231  real(RP), allocatable :: work2da(:,:) !GYL
2232  real(RP), allocatable :: work3dn(:,:,:,:) !GYL
2233  real(RP), allocatable :: work2dn(:,:,:) !GYL
2234  real(RP), allocatable :: work3dg(:,:,:,:)
2235  real(RP), allocatable :: work2dg(:,:,:)
2236 
2237  real(RP), allocatable :: hdxf(:,:)
2238  real(RP), allocatable :: rdiag(:)
2239  real(RP), allocatable :: rloc(:)
2240  real(RP), allocatable :: dep(:)
2241  real(RP), allocatable :: depd(:) !GYL
2242 
2243  integer :: ctype_merge(nid_obs,nobtype)
2244 
2245  integer :: var_local_n2nc_max
2246  integer :: var_local_n2nc(nv3d+nv2d)
2247  integer :: var_local_n2n(nv3d+nv2d)
2248  logical :: var_local_n2n_found
2249  integer :: n2n, n2nc
2250 
2251  real(RP) :: parm
2252  real(RP), allocatable :: trans(:,:,:)
2253  real(RP), allocatable :: transm(:,:)
2254  real(RP), allocatable :: transmd(:,:)
2255  real(RP), allocatable :: pa(:,:,:)
2256  real(RP) :: transrlx(nmem,nmem)
2257  logical :: trans_done(nv3d+nv2d)
2258 
2259  integer :: ij,ilev,n,m,i,j,k,nobsl
2260  integer :: nobsl_t(nid_obs,nobtype) !GYL
2261  real(RP) :: cutd_t(nid_obs,nobtype) !GYL
2262  real(RP) :: beta !GYL
2263  real(RP) :: tmpinfl !GYL
2264  real(RP) :: q_mean,q_sprd !GYL
2265  real(RP) :: q_anal(nmem) !GYL
2266 
2267  integer :: mshuf,ierr !GYL
2268  integer :: ishuf(nmem) !GYL
2269  real(RP), allocatable :: addinfl_weight(:) !GYL
2270  real(RP) :: rdx,rdy,rdxy,ref_min_dist !GYL
2271  integer :: ic,ic2,iob !GYL
2272 
2273  integer, allocatable :: search_q0(:,:,:)
2274 
2275  log_progress(*) 'data-assimilation / LETKF / system'
2276 
2277  ! -- prepare the first-guess data
2278  do j = 1, nlat
2279  do i = 1, nlon
2280  do k = 1, nlev
2281  v3dg(k,i,j,iv3d_u ) = u(k,i,j)
2282  v3dg(k,i,j,iv3d_v ) = v(k,i,j)
2283  v3dg(k,i,j,iv3d_w ) = w(k,i,j)
2284  v3dg(k,i,j,iv3d_t ) = temp(k,i,j)
2285  v3dg(k,i,j,iv3d_p ) = pres(k,i,j)
2286  v3dg(k,i,j,iv3d_q ) = qv(k,i,j)
2287  v3dg(k,i,j,iv3d_qc) = qc(k,i,j)
2288  v3dg(k,i,j,iv3d_qr) = qr(k,i,j)
2289  v3dg(k,i,j,iv3d_qi) = qi(k,i,j)
2290  v3dg(k,i,j,iv3d_qs) = qs(k,i,j)
2291  v3dg(k,i,j,iv3d_qg) = qg(k,i,j)
2292  end do
2293  end do
2294  end do
2295  do j = 1, nlat
2296  do i = 1, nlon
2297  v2dg(i,j,:) = undef ! tentative
2298  end do
2299  end do
2300 
2301  if( letkf_debug_log ) then
2302  call monit_obs_mpi( obs_in_num, obs_in_format, v3dg, v2dg, monit_step=1 )
2303  endif
2304 
2305  call scatter_grd_mpi_all2all( 1, nens, v3dg, v2dg, gues3d(:,:,1:nens,:), gues2d(:,1:nens,:) )
2306 
2307  ! -- obtain the ensemble mean
2308  do n = 1, nv3d
2309  do k = 1, nlev
2310  do ij = 1, nij1
2311  gues3d(ij,k,mmean,n) = average( nmem, gues3d(ij,k,1:nmem,n), undef )
2312  end do
2313  end do
2314  end do
2315  do n = 1, nv2d
2316  do ij = 1, nij1
2317  gues2d(ij,mmean,n) = average( nmem, gues2d(ij,1:nmem,n), undef )
2318  end do
2319  end do
2320 
2321  !
2322  ! Variable localization
2323  !
2324  var_local(:,1) = var_local_uv(:)
2325  var_local(:,2) = var_local_t(:)
2326  var_local(:,3) = var_local_q(:)
2327  var_local(:,4) = var_local_ps(:)
2328  var_local(:,5) = var_local_rain(:)
2329  var_local(:,6) = var_local_tc(:)
2330  var_local(:,7) = var_local_radar_ref(:)
2331  var_local(:,8) = var_local_radar_vr(:)
2332 
2333  var_local_n2nc_max = 1
2334  var_local_n2nc(1) = 1
2335  var_local_n2n(1) = 1
2336 
2337  do n = 2, nv3d+nv2d
2338  var_local_n2n_found = .false.
2339  do i = 1, var_local_n2nc_max
2340  !if (maxval(abs(var_local(var_local_n2nc(i),:) - var_local(n,:))) < tiny(var_local(1,1))) then
2341  if (all(var_local(var_local_n2nc(i),:) == var_local(n,:))) then
2342  var_local_n2nc(n) = var_local_n2nc(i)
2343  var_local_n2n(n) = var_local_n2n(var_local_n2nc(n))
2344  var_local_n2n_found = .true.
2345  exit
2346  end if
2347  end do
2348  if( .NOT. var_local_n2n_found ) then
2349  var_local_n2nc_max = var_local_n2nc_max + 1
2350  var_local_n2nc(n) = var_local_n2nc_max
2351  var_local_n2n(n) = n
2352  end if
2353  end do
2354 
2355  !
2356  ! Observation number limit (*to be moved to namelist*)
2357  !
2358  ctype_merge(:,:) = 0
2359  ctype_merge(uid_obs(id_radar_ref_obs),22) = 1
2360  ctype_merge(uid_obs(id_radar_ref_zero_obs),22) = 1
2361 
2362  allocate (n_merge(nctype))
2363  allocate (ic_merge(nid_obs*nobtype,nctype))
2364  n_merge(:) = 1
2365  do ic = 1, nctype
2366  if (n_merge(ic) > 0) then
2367  ic_merge(1,ic) = ic
2368  if (ctype_merge(elm_u_ctype(ic),typ_ctype(ic)) > 0) then
2369  do ic2 = ic+1, nctype
2370  if (ctype_merge(elm_u_ctype(ic2),typ_ctype(ic2)) == ctype_merge(elm_u_ctype(ic),typ_ctype(ic))) then
2371  n_merge(ic) = n_merge(ic) + 1
2372  ic_merge(n_merge(ic),ic) = ic2
2373  n_merge(ic2) = 0
2374  end if
2375  end do
2376  end if
2377  end if
2378  end do
2379  n_merge_max = maxval(n_merge)
2380 
2381  allocate (search_q0(nctype,nv3d+1,nij1))
2382  search_q0(:,:,:) = 1
2383 
2384  radar_only = .true.
2385  do ic = 1, nctype
2386  if (obtypelist(typ_ctype(ic)) /= 'PHARAD') then
2387  radar_only = .false.
2388  exit
2389  end if
2390  end do
2391  !
2392  ! FCST PERTURBATIONS
2393  !
2394  ! .... this has been done by write_ensmean in letkf.f90
2395  ! CALL ensmean_grd(nens,nij1,gues3d,gues2d,mean3d,mean2d)
2396  DO n=1,nv3d
2397  DO m=1,nmem
2398  DO k=1,nlev
2399  DO i=1,nij1
2400  gues3d(i,k,m,n) = gues3d(i,k,m,n) - gues3d(i,k,mmean,n)
2401  END DO
2402  END DO
2403  END DO
2404  END DO
2405  DO n=1,nv2d
2406  DO m=1,nmem
2407  DO i=1,nij1
2408  gues2d(i,m,n) = gues2d(i,m,n) - gues2d(i,mmean,n)
2409  END DO
2410  END DO
2411  END DO
2412 
2413  !
2414  ! multiplicative inflation
2415  !
2416  IF(infl_mul > 0.0d0) THEN ! fixed multiplicative inflation parameter
2417  work3d = infl_mul
2418  work2d = infl_mul
2419  ELSE ! 3D parameter values are read-in
2420  log_error("LETKF_system",*) 'This system has not been implemented yet. INFL_MUL must be greather than 0.0.'
2421  call prc_abort
2422  ! not used now (TODO)
2423  ! allocate (work3dg(nlon,nlat,nlev,nv3d))
2424  ! allocate (work2dg(nlon,nlat,nv2d))
2425  ! IF(myrank_e == mmean_rank_e) THEN
2426  ! call read_restart(INFL_MUL_IN_BASENAME,work3dg,work2dg)
2427  ! END IF
2428  !
2429  ! CALL scatter_grd_mpi(mmean_rank_e,work3dg,work2dg,work3d,work2d)
2430  !
2431  END IF
2432  IF(infl_mul_min > 0.0d0) THEN
2433  work3d = max(work3d, infl_mul_min)
2434  work2d = max(work2d, infl_mul_min)
2435  END IF
2436 
2437  ! This loop cannot use OpenMP on FUGAKU (T. Honda, as of 10/16/2020)
2438  allocate( hdxf( nobstotal, nmem ) )
2439  allocate (rdiag(nobstotal))
2440  allocate (rloc(nobstotal))
2441  allocate (dep(nobstotal))
2442  if( ens_with_mdet ) then
2443  allocate (depd(nobstotal))
2444  end if
2445  allocate (trans(nmem,nmem,var_local_n2nc_max))
2446  allocate (transm(nmem, var_local_n2nc_max))
2447  allocate (transmd(nmem, var_local_n2nc_max))
2448  allocate (pa(nmem,nmem,var_local_n2nc_max))
2449 
2450  !
2451  ! MAIN ASSIMILATION LOOP
2452  !
2453  DO ilev=1,nlev
2454  DO ij=1,nij1
2455 
2456  trans_done(:) = .false. !GYL
2457 
2458  ! weight parameter based on grid locations (not for covariance inflation purpose)
2459  ! if the weight is zero, no analysis update is needed
2460  call relax_beta(rig1(ij),rjg1(ij),hgt1(ij,ilev),beta)
2461 
2462  if (beta == 0.0d0) then
2463  do n = 1, nv3d
2464  do m = 1, nmem
2465  anal3d(ij,ilev,m,n) = gues3d(ij,ilev,mmean,n) + gues3d(ij,ilev,m,n)
2466  end do
2467  if (ens_with_mdet) then
2468  anal3d(ij,ilev,mmdet,n) = gues3d(ij,ilev,mmdet,n)
2469  end if
2470  end do
2471  if (ilev == 1) then
2472  do n = 1, nv2d
2473  do m = 1, nmem
2474  anal2d(ij,m,n) = gues2d(ij,mmean,n) + gues2d(ij,m,n)
2475  end do
2476  if (ens_with_mdet) then
2477  anal2d(ij,mmdet,n) = gues2d(ij,mmdet,n)
2478  end if
2479  end do
2480  end if
2481 
2482  cycle
2483  end if
2484 
2485  ! update 3D variables
2486  DO n=1,nv3d
2487 
2488  n2nc = var_local_n2nc(n)
2489  n2n = var_local_n2n(n)
2490 
2491  if (gues3d(ij,ilev,mmean,iv3d_p) < q_update_top .and. n >= iv3d_q .and. n <= iv3d_qg) then !GYL - Upper bound of Q update levels
2492  do m = 1, nmem !GYL
2493  anal3d(ij,ilev,m,n) = gues3d(ij,ilev,mmean,n) + gues3d(ij,ilev,m,n) !GYL
2494  end do !GYL
2495  if (ens_with_mdet) then !GYL
2496  anal3d(ij,ilev,mmdet,n) = gues3d(ij,ilev,mmdet,n) !GYL
2497  end if !GYL
2498 
2499  cycle !GYL
2500  end if !GYL
2501 
2502  if (relax_to_inflated_prior) then
2503  parm = work3d(ij,ilev,n)
2504  else
2505  parm = 1.0d0
2506  end if
2507 
2508  ! calculate mean and perturbation weights
2509  if (trans_done(n2nc)) then
2510  ! if weights already computed for other variables can be re-used(no variable localization), do not need to compute again
2511  if (infl_mul_adaptive) then
2512  work3d(ij,ilev,n) = work3d(ij,ilev,n2n)
2513  end if
2514  if (nobs_out) then
2515  work3dn(:,ij,ilev,n) = work3dn(:,ij,ilev,n2n)
2516  end if
2517 
2518  ELSE
2519  ! compute weights with localized observations
2520  if (ens_with_mdet) then !GYL
2521  CALL obs_local(obs(:),rig1(ij),rjg1(ij),gues3d(ij,ilev,mmean,iv3d_p),hgt1(ij,ilev),n,& !GYL
2522  hdxf,rdiag,rloc,dep,nobsl,depd=depd,nobsl_t=nobsl_t,cutd_t=cutd_t,srch_q0=search_q0(:,n,ij)) !GYL
2523  else !GYL
2524  CALL obs_local(obs(:),rig1(ij),rjg1(ij),gues3d(ij,ilev,mmean,iv3d_p),hgt1(ij,ilev),n,& !GYL
2525  hdxf,rdiag,rloc,dep,nobsl,nobsl_t=nobsl_t,cutd_t=cutd_t,srch_q0=search_q0(:,n,ij)) !GYL
2526  end if !GYL
2527  IF(relax_alpha_spread /= 0.0d0) THEN !GYL
2528  if (ens_with_mdet) then !GYL
2529  CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work3d(ij,ilev,n), & !GYL
2530  trans(:,:,n2nc),transm=transm(:,n2nc),pao=pa(:,:,n2nc), & !GYL
2531  rdiag_wloc=.true.,infl_update=infl_mul_adaptive, & !GYL
2532  depd=depd,transmd=transmd(:,n2nc)) !GYL
2533  else !GYL
2534  CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work3d(ij,ilev,n), & !GYL
2535  trans(:,:,n2nc),transm=transm(:,n2nc),pao=pa(:,:,n2nc), & !GYL
2536  rdiag_wloc=.true.,infl_update=infl_mul_adaptive) !GYL
2537  end if !GYL
2538  ELSE !GYL
2539  if (ens_with_mdet) then !GYL
2540  CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work3d(ij,ilev,n), & !GYL
2541  trans(:,:,n2nc),transm=transm(:,n2nc), & !GYL
2542  rdiag_wloc=.true.,infl_update=infl_mul_adaptive, & !GYL
2543  depd=depd,transmd=transmd(:,n2nc)) !GYL
2544  else !GYL
2545  CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work3d(ij,ilev,n), & !GYL
2546  trans(:,:,n2nc),transm=transm(:,n2nc), & !GYL
2547  rdiag_wloc=.true.,infl_update=infl_mul_adaptive) !GYL
2548  end if !GYL
2549  END IF !GYL
2550  trans_done(n2nc) = .true. !GYL
2551  IF(nobs_out) THEN !GYL
2552  work3dn(:,ij,ilev,n) = real(sum(nobsl_t, dim=1),rp) !GYL !!! NOBS: sum over all variables for each report type
2553  work3dn(nobtype+1,ij,ilev,n) = real(nobsl_t(9,22),rp) !GYL !!! NOBS: ref
2554  work3dn(nobtype+2,ij,ilev,n) = real(nobsl_t(10,22),rp) !GYL !!! NOBS: re0
2555  work3dn(nobtype+3,ij,ilev,n) = real(nobsl_t(11,22),rp) !GYL !!! NOBS: vr
2556  work3dn(nobtype+4,ij,ilev,n) = real(cutd_t(9,22),rp) !GYL !!! CUTOFF_DIST: ref
2557  work3dn(nobtype+5,ij,ilev,n) = real(cutd_t(10,22),rp) !GYL !!! CUTOFF_DIST: re0
2558  work3dn(nobtype+6,ij,ilev,n) = real(cutd_t(11,22),rp) !GYL !!! CUTOFF_DIST: vr
2559  END IF !GYL
2560 
2561  END IF
2562 
2563  ! relaxation via LETKF weight
2564  IF(relax_alpha /= 0.0d0) THEN !GYL - RTPP method (Zhang et al. 2004)
2565  CALL weight_rtpp(trans(:,:,n2nc),parm,transrlx) !GYL
2566  ELSE IF(relax_alpha_spread /= 0.0d0) THEN !GYL - RTPS method (Whitaker and Hamill 2012)
2567  IF(relax_spread_out) THEN !GYL
2568  CALL weight_rtps(trans(:,:,n2nc),pa(:,:,n2nc),gues3d(ij,ilev,:,n), & !GYL
2569  parm,transrlx,work3da(ij,ilev,n)) !GYL
2570  ELSE !GYL
2571  CALL weight_rtps(trans(:,:,n2nc),pa(:,:,n2nc),gues3d(ij,ilev,:,n), & !GYL
2572  parm,transrlx,tmpinfl) !GYL
2573  END IF !GYL
2574  ELSE !GYL
2575  transrlx = trans(:,:,n2nc) !GYL - No relaxation
2576  END IF !GYL
2577 
2578  ! total weight matrix
2579  do m = 1, nmem
2580  do k = 1, nmem
2581  transrlx(k,m) = (transrlx(k,m) + transm(k,n2nc)) * beta
2582  end do
2583  transrlx(m,m) = transrlx(m,m) + (1.0_rp-beta)
2584  end do
2585  ! analysis update of members
2586  do m = 1, nmem
2587  anal3d(ij,ilev,m,n) = gues3d(ij,ilev,mmean,n)
2588  do k = 1, nmem
2589  anal3d(ij,ilev,m,n) = anal3d(ij,ilev,m,n) + gues3d(ij,ilev,k,n) * transrlx(k,m)
2590  end do
2591  end do
2592 
2593  if( ens_with_mdet ) then
2594  ! analysis update of deterministic run
2595  anal3d(ij,ilev,mmdet,n) = 0.0_rp
2596  do k = 1, nmem
2597  anal3d(ij,ilev,mmdet,n) = anal3d(ij,ilev,mmdet,n) + gues3d(ij,ilev,k,n) * transmd(k,n2nc)
2598  end do
2599  anal3d(ij,ilev,mmdet,n) = gues3d(ij,ilev,mmdet,n) + anal3d(ij,ilev,mmdet,n) * beta
2600  end if
2601 
2602  ! limit q spread
2603  IF(q_sprd_max > 0.0d0 .and. n == iv3d_q) THEN !GYL
2604  q_mean = sum(anal3d(ij,ilev,1:nmem,n)) / real(nmem,rp) !GYL
2605  q_sprd = 0.0d0 !GYL
2606  DO m=1,nmem !GYL
2607  q_anal(m) = anal3d(ij,ilev,m,n) - q_mean !GYL
2608  q_sprd = q_sprd + q_anal(m)**2 !GYL
2609  END DO !GYL
2610 
2611  if ( q_mean > 0.0_rp ) then
2612  q_sprd = sqrt(q_sprd / real(nmem-1,rp)) / q_mean !GYL
2613  IF(q_sprd > q_sprd_max) THEN !GYL
2614  DO m=1,nmem !GYL
2615  anal3d(ij,ilev,m,n) = q_mean + q_anal(m) * q_sprd_max / q_sprd !GYL
2616  END DO !GYL
2617  END IF !GYL
2618  endif
2619  END IF !GYL
2620 
2621  END DO ! [ n=1,nv3d ]
2622 
2623  ! update 2D variables at ilev = 1
2624  IF(ilev == 1) THEN
2625 
2626  DO n=1,nv2d
2627 
2628  n2nc = var_local_n2nc(nv3d+n)
2629  n2n = var_local_n2n(nv3d+n)
2630 
2631  if (relax_to_inflated_prior) then
2632  parm = work2d(ij,n)
2633  else
2634  parm = 1.0d0
2635  end if
2636 
2637  ! calculate mean and perturbation weights
2638  if (trans_done(n2nc)) then
2639  ! if weights already computed for other variables can be re-used(no variable localization), do not need to compute again
2640  IF(n2n <= nv3d) then
2641  if (infl_mul_adaptive) then
2642  work2d(ij,n) = work3d(ij,ilev,n2n)
2643  end if
2644  if (nobs_out) then
2645  work2dn(:,ij,n) = work3dn(:,ij,ilev,n2n)
2646  end if
2647  else
2648  if (infl_mul_adaptive) then
2649  work2d(ij,n) = work2d(ij,n2n-nv3d)
2650  end if
2651  if (nobs_out) then
2652  work2dn(:,ij,n) = work2dn(:,ij,n2n-nv3d)
2653  end if
2654  end if
2655 
2656  ELSE
2657  ! compute weights with localized observations
2658  if (ens_with_mdet) then !GYL
2659  CALL obs_local(obs(:),rig1(ij),rjg1(ij),gues3d(ij,ilev,mmean,iv3d_p),hgt1(ij,ilev),nv3d+n,hdxf,rdiag,rloc,dep,nobsl,depd=depd,nobsl_t=nobsl_t,cutd_t=cutd_t,srch_q0=search_q0(:,nv3d+1,ij))
2660  else !GYL
2661  CALL obs_local(obs(:),rig1(ij),rjg1(ij),gues3d(ij,ilev,mmean,iv3d_p),hgt1(ij,ilev),nv3d+n,hdxf,rdiag,rloc,dep,nobsl,nobsl_t=nobsl_t,cutd_t=cutd_t,srch_q0=search_q0(:,nv3d+1,ij))
2662  end if !GYL
2663  IF(relax_alpha_spread /= 0.0d0) THEN !GYL
2664  if (ens_with_mdet) then !GYL
2665  CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work2d(ij,n), & !GYL
2666  trans(:,:,n2nc),transm=transm(:,n2nc),pao=pa(:,:,n2nc), & !GYL
2667  rdiag_wloc=.true.,infl_update=infl_mul_adaptive, & !GYL
2668  depd=depd,transmd=transmd(:,n2nc)) !GYL
2669  else !GYL
2670  CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work2d(ij,n), & !GYL
2671  trans(:,:,n2nc),transm=transm(:,n2nc),pao=pa(:,:,n2nc), & !GYL
2672  rdiag_wloc=.true.,infl_update=infl_mul_adaptive) !GYL
2673  end if !GYL
2674  ELSE !GYL
2675  if (ens_with_mdet) then !GYL
2676  CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work2d(ij,n), & !GYL
2677  trans(:,:,n2nc),transm=transm(:,n2nc), & !GYL
2678  rdiag_wloc=.true.,infl_update=infl_mul_adaptive, & !GYL
2679  depd=depd,transmd=transmd(:,n2nc)) !GYL
2680  else !GYL
2681  CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work2d(ij,n), & !GYL
2682  trans(:,:,n2nc),transm=transm(:,n2nc), & !GYL
2683  rdiag_wloc=.true.,infl_update=infl_mul_adaptive) !GYL
2684  end if !GYL
2685  END IF !GYL
2686  trans_done(n2nc) = .true. !GYL
2687  IF(nobs_out) THEN !GYL
2688  work2dn(:,ij,n) = real(sum(nobsl_t,dim=1),rp) !GYL !!! NOBS: sum over all variables for each report type
2689  END IF !GYL
2690 
2691  END IF
2692 
2693  ! relaxation via LETKF weight
2694  IF(relax_alpha /= 0.0d0) THEN !GYL - RTPP method (Zhang et al. 2004)
2695  CALL weight_rtpp(trans(:,:,n2nc),parm,transrlx) !GYL
2696  ELSE IF(relax_alpha_spread /= 0.0d0) THEN !GYL - RTPS method (Whitaker and Hamill 2012)
2697  IF(relax_spread_out) THEN !GYL
2698  CALL weight_rtps(trans(:,:,n2nc),pa(:,:,n2nc),gues2d(ij,:,n), & !GYL
2699  parm,transrlx,work2da(ij,n)) !GYL
2700  ELSE !GYL
2701  CALL weight_rtps(trans(:,:,n2nc),pa(:,:,n2nc),gues2d(ij,:,n), & !GYL
2702  parm,transrlx,tmpinfl) !GYL
2703  END IF !GYL
2704  ELSE !GYL
2705  transrlx = trans(:,:,n2nc) !GYL - No relaxation
2706  END IF !GYL
2707 
2708  ! total weight matrix
2709  do m = 1, nmem
2710  do k = 1, nmem
2711  transrlx(k,m) = (transrlx(k,m) + transm(k,n2nc)) * beta
2712  end do
2713  transrlx(m,m) = transrlx(m,m) + (1.0_rp-beta)
2714  end do
2715 
2716  ! analysis update of members
2717  do m = 1, nmem
2718  anal2d(ij,m,n) = gues2d(ij,mmean,n)
2719  do k = 1, nmem
2720  anal2d(ij,m,n) = anal2d(ij,m,n) + gues2d(ij,k,n) * transrlx(k,m)
2721  end do
2722  end do
2723 
2724  ! analysis update of deterministic run
2725  if (ens_with_mdet) then
2726  anal2d(ij,mmdet,n) = 0.0_rp
2727  do k = 1, nmem
2728  anal2d(ij,mmdet,n) = anal2d(ij,mmdet,n) + gues2d(ij,k,n) * transmd(k,n2nc)
2729  end do
2730  anal2d(ij,mmdet,n) = gues2d(ij,mmdet,n) + anal2d(ij,mmdet,n) * beta
2731  end if
2732 
2733  END DO
2734  END IF
2735  END DO
2736  END DO ! [ ilev=1,nlev ]
2737 
2738  deallocate (hdxf,rdiag,rloc,dep)
2739  if (ens_with_mdet) then
2740  deallocate (depd)
2741  end if
2742  deallocate (trans,transm,transmd,pa)
2743 
2744  deallocate (n_merge,ic_merge)
2745  deallocate (search_q0)
2746 
2747  IF (allocated(work3dg)) deallocate (work3dg)
2748  IF (allocated(work2dg)) deallocate (work2dg)
2749 
2750  !
2751  ! Additive inflation
2752  !
2753  IF(infl_add > 0.0d0) THEN
2754 
2755  if (infl_add_q_ratio) then
2756  work3d(:,:,:) = gues3d(:,:,mmean,:)
2757  else
2758  work3d(:,:,:) = 1.0d0
2759  end if
2760 
2761  allocate (addinfl_weight(nij1))
2762  if (infl_add_ref_only) then
2763  addinfl_weight(:) = 0.0d0
2764  ic = ctype_elmtyp(uid_obs(id_radar_ref_obs), 22)
2765  if (ic > 0) then
2766  do ij = 1, nij1
2767  ref_min_dist = 1.0d33
2768  !!!!!! save this (ref_min_dist) information when doing DA
2769  do iob = obsgrd(ic)%ac_ext(0, 1) + 1, obsgrd(ic)%ac_ext(obsgrd(ic)%ngrdext_i, obsgrd(ic)%ngrdext_j)
2770  rdx = (rig1(ij) - obs(obsda_sort%set(iob))%ri(obsda_sort%idx(iob))) * dx
2771  rdy = (rjg1(ij) - obs(obsda_sort%set(iob))%rj(obsda_sort%idx(iob))) * dy
2772  rdxy = rdx*rdx + rdy*rdy
2773  if (rdxy < ref_min_dist) then
2774  ref_min_dist = rdxy
2775  end if
2776  end do
2777 
2778  ref_min_dist = ref_min_dist / (hori_loc_ctype(ic) * hori_loc_ctype(ic))
2779  if (ref_min_dist <= dist_zero_fac_square) then
2780  addinfl_weight(ij) = exp(-0.5d0 * ref_min_dist)
2781  end if
2782  end do
2783  end if
2784  else
2785  addinfl_weight(:) = 1.0d0
2786  end if
2787 
2788  log_info("LETKF_system",*) 'Additive covariance inflation, parameter:', infl_add
2789  if (infl_add_shuffle) then
2790  if (rank_ens== 0) then
2791  call knuth_shuffle(nens, ishuf)
2792  end if
2793  call mpi_bcast(ishuf, nens, mpi_integer, 0, comm_ens, ierr)
2794  log_info("LETKF_system",*) 'shuffle members: on'
2795  log_info("LETKF_system",*) 'shuffle sequence:', ishuf
2796  end if
2797 
2798  DO n=1,nv3d
2799  DO m=1,nens
2800  if (infl_add_shuffle) then
2801  mshuf = ishuf(m)
2802  else
2803  mshuf = m
2804  end if
2805  if (n == iv3d_q .or. n == iv3d_qc .or. n == iv3d_qr .or. n == iv3d_qi .or. n == iv3d_qs .or. n == iv3d_qg) then
2806  DO k=1,nlev
2807  DO i=1,nij1
2808  anal3d(i,k,m,n) = anal3d(i,k,m,n) &
2809  & + addi3d(i,k,mshuf,n) * infl_add * addinfl_weight(i) * work3d(i,k,n)
2810  END DO
2811  END DO
2812  else
2813  DO k=1,nlev
2814  DO i=1,nij1
2815  anal3d(i,k,m,n) = anal3d(i,k,m,n) &
2816  & + addi3d(i,k,mshuf,n) * infl_add * addinfl_weight(i)
2817  END DO
2818  END DO
2819  end if
2820  END DO
2821  END DO
2822  DO n=1,nv2d
2823  DO m=1,nens
2824  if (infl_add_shuffle) then
2825  mshuf = ishuf(m)
2826  else
2827  mshuf = m
2828  end if
2829  DO i=1,nij1
2830  anal2d(i,m,n) = anal2d(i,m,n) + addi2d(i,mshuf,n) * infl_add * addinfl_weight(i)
2831  END DO
2832  END DO
2833  END DO
2834 
2835  deallocate (addinfl_weight)
2836 
2837  END IF
2838 
2839  ! -- prepare the output values
2840  call gather_grd_mpi_all2all( 1, nens, anal3d(:,:,1:nens,:), anal2d(:,1:nens,:), v3dg, v2dg )
2841 
2842  do j = 1, nlat
2843  do i = 1, nlon
2844  do k = 1, nlev
2845  u(k,i,j) = v3dg(k,i,j,iv3d_u )
2846  v(k,i,j) = v3dg(k,i,j,iv3d_v )
2847  w(k,i,j) = v3dg(k,i,j,iv3d_w )
2848  temp(k,i,j) = v3dg(k,i,j,iv3d_t )
2849  pres(k,i,j) = v3dg(k,i,j,iv3d_p )
2850  qv(k,i,j) = v3dg(k,i,j,iv3d_q )
2851  qc(k,i,j) = v3dg(k,i,j,iv3d_qc)
2852  qr(k,i,j) = v3dg(k,i,j,iv3d_qr)
2853  qi(k,i,j) = v3dg(k,i,j,iv3d_qi)
2854  qs(k,i,j) = v3dg(k,i,j,iv3d_qs)
2855  qg(k,i,j) = v3dg(k,i,j,iv3d_qg)
2856  end do
2857  end do
2858  end do
2859 
2860  if( letkf_debug_log ) then
2861  call monit_obs_mpi( obs_in_num, obs_in_format, v3dg, v2dg, monit_step=2 )
2862  endif
2863 
2864  return

References scale_const::const_undef, scale_matrix::matrix_solver_eigenvalue_decomposition(), nid_obs, nobtype, nv2d, nv3d, scale_prc::prc_abort(), scale_random::random_knuth_shuffle(), and scale_precision::rp.

Referenced by mod_da_driver::da_driver_update().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ letkf_param_estimation_system()

subroutine, public scale_letkf::letkf_param_estimation_system ( integer, intent(in)  PEST_PMAX,
real(rp), dimension(nens,pest_pmax), intent(inout)  PEST_VAR0 
)

Definition at line 2873 of file scale_letkf.F90.

2873  use scale_const, only: &
2874  undef => const_undef
2875  use scale_statistics, only: &
2876  average => statistics_average
2877  implicit none
2878 
2879  integer, intent(in) :: PEST_PMAX
2880 
2881  real(RP), intent(inout) :: PEST_VAR0(nens,PEST_PMAX)
2882 
2883  real(RP) :: gues0d(nens+1,PEST_PMAX)
2884  real(RP) :: anal0d(nens+1,PEST_PMAX)
2885 
2886  real(RP) :: work0d(PEST_PMAX)
2887 
2888  real(RP),allocatable :: hdxf(:,:)
2889  real(RP),allocatable :: rdiag(:)
2890  real(RP),allocatable :: rloc(:)
2891  real(RP),allocatable :: dep(:)
2892  real(RP),allocatable :: depd(:) !gyl
2893 
2894  real(RP) :: parm
2895  real(RP),allocatable :: trans(:,:)
2896  real(RP),allocatable :: transm(:)
2897  real(RP),allocatable :: transmd(:)
2898  real(RP),allocatable :: pa(:,:)
2899  real(RP) :: transrlx(nmem,nmem)
2900 
2901  integer :: n,m,k,nobsl
2902  real(RP) :: beta !gyl
2903  real(RP) :: tmpinfl !gyl
2904  integer :: ierr
2905 
2906  log_progress(*) 'data-assimilation / LETKF / parameter estimation'
2907 
2908  ! -- prepare the first-guess data
2909  do n = 1, pest_pmax
2910  do m = 1, nens
2911  gues0d(m,n) = pest_var0(m,n)
2912  end do
2913  end do
2914 
2915  ! -- obtain the ensemble mean
2916  do n = 1, pest_pmax
2917  gues0d(mmean,n) = average( nmem, gues0d(1:nmem,n), undef )
2918  end do
2919 
2920  !
2921  ! No variable localization
2922  !
2923 
2924  !
2925  ! FCST PERTURBATIONS
2926  !
2927 
2928  do n = 1, pest_pmax
2929  do m = 1, nmem
2930  gues0d(m,n) = gues0d(m,n) - gues0d(mmean,n)
2931  end do
2932  end do
2933 
2934  allocate (hdxf(nobstotal,nmem))
2935  allocate (rdiag(nobstotal))
2936  allocate (rloc(nobstotal))
2937  allocate (dep(nobstotal))
2938  if (ens_with_mdet) then
2939  allocate (depd(nobstotal))
2940  allocate (transmd(nmem))
2941  end if
2942  allocate (trans(nmem,nmem))
2943  allocate (transm(nmem))
2944  allocate (pa(nmem,nmem))
2945 
2946  !
2947  ! MAIN ASSIMILATION LOOP
2948  !
2949  parm = 1.0_rp
2950  beta = 1.0_rp
2951 
2952  work0d = 1.0_rp
2953 
2954  DO n = 1, pest_pmax
2955  ! calculate mean and perturbation weights
2956  ! compute weights with localized observations
2957  if (ens_with_mdet) then
2958  CALL obs_pest_etkf(hdxf, rdiag, rloc, dep, nobsl, depd=depd)
2959  else
2960  CALL obs_pest_etkf(hdxf, rdiag, rloc, dep, nobsl)
2961  end if
2962 
2963  if (ens_with_mdet) then
2964  CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work0d(n), &
2965  trans(:,:),transm=transm(:),pao=pa(:,:), &
2966  rdiag_wloc=.true.,&
2967  depd=depd,transmd=transmd(:))
2968  else
2969  CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work0d(n), &
2970  trans(:,:),transm=transm(:),pao=pa(:,:), &
2971  rdiag_wloc=.true.)
2972  end if
2973 
2974  ! relaxation via LETKF weight
2975  CALL weight_rtps_const(trans(:,:),pa(:,:),gues0d(:,n),transrlx,tmpinfl)
2976 
2977  ! total weight matrix
2978  do m = 1, nmem
2979  do k = 1, nmem
2980  transrlx(k,m) = (transrlx(k,m) + transm(k)) * beta
2981  end do
2982  transrlx(m,m) = transrlx(m,m) + (1.0_rp-beta)
2983  end do
2984 
2985  ! analysis update of members
2986  do m = 1, nmem
2987  anal0d(m,n) = gues0d(mmean,n)
2988  do k = 1, nmem
2989  anal0d(m,n) = anal0d(m,n) + gues0d(k,n) * transrlx(k,m)
2990  end do
2991  end do
2992 
2993  ! analysis update of deterministic run
2994  if (ens_with_mdet) then
2995  anal0d(mmdet,n) = 0.0_rp
2996  do k = 1, nmem
2997  anal0d(mmdet,n) = anal0d(mmdet,n) + gues0d(k,n) * transmd(k)
2998  end do
2999  anal0d(mmdet,n) = gues0d(mmdet,n) + anal0d(mmdet,n) * beta
3000  end if
3001 
3002  enddo
3003 
3004  do n = 1, pest_pmax
3005  do m = 1, nens
3006  pest_var0(m,n) = anal0d(m,n)
3007  end do
3008  end do
3009 
3010  deallocate (hdxf,rdiag,rloc,dep)
3011  if (ens_with_mdet) then
3012  deallocate (depd,transmd)
3013  end if
3014  deallocate (trans,transm,pa)
3015 
3016  return

References scale_const::const_undef.

Referenced by mod_da_param_estimation::da_param_estimation_update().

Here is the caller graph for this function:

◆ letkf_add_inflation_setup()

subroutine, public scale_letkf::letkf_add_inflation_setup ( real(rp), dimension(:,:,:,:), intent(out)  addi3d,
real(rp), dimension(:,:,:), intent(out)  addi2d 
)

Definition at line 3025 of file scale_letkf.F90.

3025  implicit none
3026 
3027  real(RP), intent(out) :: addi3d(:,:,:,:)
3028  real(RP), intent(out) :: addi2d(:,:,:)
3029  integer :: i, k, m, n
3030 
3031  call read_ens_mpi_addiinfl(addi3d, addi2d)
3032 
3033  call ensmean_grd(nens, nens, nij1, addi3d, addi2d)
3034 
3035  do n = 1, nv3d
3036  do m = 1, nens
3037  do k = 1, nlev
3038  do i = 1, nij1
3039  addi3d(i,k,m,n) = addi3d(i,k,m,n) - addi3d(i,k,mmean,n)
3040  end do
3041  end do
3042  end do
3043  end do
3044 
3045  do n = 1, nv2d
3046  do m = 1, nens
3047  do i = 1, nij1
3048  addi2d(i,m,n) = addi2d(i,m,n) - addi2d(i,mmean,n)
3049  end do
3050  end do
3051  end do
3052 
3053  return

References scale_da_read_pawr_toshiba::azdim, scale_da_read_mp_pawr_toshiba::azdim, scale_calendar::calendar_adjust_daysec(), scale_calendar::calendar_date2daysec(), scale_calendar::calendar_daysec2date(), scale_const::const_d2r, scale_const::const_pi, scale_const::const_r2d, scale_const::const_radius, scale_da_read_mp_pawr_toshiba::da_read_mp_pawr_toshiba(), scale_da_read_pawr_toshiba::da_read_pawr_toshiba(), scale_precision::dp, scale_da_read_mp_pawr_toshiba::eldim, scale_da_read_pawr_toshiba::eldim, scale_io::io_get_available_fid(), scale_matrix::matrix_solver_eigenvalue_decomposition(), nv2d, nv3d, scale_prc::prc_abort(), radar_superobing(), scale_da_read_mp_pawr_toshiba::rdim, scale_da_read_pawr_toshiba::rdim, read_toshiba(), read_toshiba_mpr(), scale_precision::rp, scale_sort::sort_quickselect_arg(), scale_sort::sort_quickselect_desc_arg(), and scale_precision::sp.

Here is the call graph for this function:

◆ radar_superobing()

subroutine scale_letkf::radar_superobing ( integer, intent(in)  na,
integer, intent(in)  nr,
integer, intent(in)  ne,
real(rp), dimension(na, nr, ne), intent(in)  radlon,
real(rp), dimension(na, nr, ne), intent(in)  radlat,
real(rp), dimension(na, nr, ne), intent(in)  radz,
real(rp), dimension(na, nr, ne), intent(in)  ze,
real(rp), dimension(na, nr, ne), intent(in)  vr,
real(rp), dimension(na, nr, ne), intent(in)  qcflag,
real(rp), dimension(na, nr, ne), intent(in)  attenuation,
integer, intent(in)  nlon,
integer, intent(in)  nlat,
integer, intent(in)  nlev,
real(rp), dimension(nlon), intent(in)  lon,
real(rp), dimension(nlat), intent(in)  lat,
real(rp), dimension(nlev), intent(in)  z,
real(rp), intent(in)  dlon,
real(rp), intent(in)  dlat,
real(rp), intent(in)  dz,
real(rp), intent(in)  missing,
logical, intent(in)  input_is_dbz,
real(rp), intent(in)  lon0,
real(rp), intent(in)  lat0,
integer(8), intent(out)  nobs_sp,
integer(8), dimension(:), intent(out), allocatable  grid_index,
real(rp), dimension(:), intent(out), allocatable  grid_ref,
real(rp), dimension(:), intent(out), allocatable  grid_lon_ref,
real(rp), dimension(:), intent(out), allocatable  grid_lat_ref,
real(rp), dimension(:), intent(out), allocatable  grid_z_ref,
integer(8), dimension(:), intent(out), allocatable  grid_count_ref,
real(rp), dimension(:), intent(out), allocatable  grid_vr,
real(rp), dimension(:), intent(out), allocatable  grid_lon_vr,
real(rp), dimension(:), intent(out), allocatable  grid_lat_vr,
real(rp), dimension(:), intent(out), allocatable  grid_z_vr,
integer(8), dimension(:), intent(out), allocatable  grid_count_vr 
)

Definition at line 4404 of file scale_letkf.F90.

4404  use scale_sort, only: &
4405  sort_uniq_int_sorted
4406  implicit none
4407 
4408  integer, intent(in) :: na, nr, ne ! array size of the spherical grid
4409  real(RP), intent(in), dimension(na, nr, ne) :: radlon, radlat, radz ! georeference
4410  real(RP), intent(in) :: ze(na, nr, ne), vr(na, nr, ne) ! main data
4411  real(RP), intent(in) :: qcflag(na, nr, ne), attenuation(na, nr, ne) ! additional info
4412  integer, intent(in) :: nlon, nlat, nlev ! array size of the cartesian grid
4413  real(RP), intent(in) :: lon(nlon), lat(nlat), z(nlev)
4414  real(RP), intent(in) :: dlon, dlat, dz, missing
4415  logical, intent(in) :: input_is_dbz
4416  real(RP), intent(in) :: lon0, lat0
4417 
4418  integer(8), intent(out) :: nobs_sp
4419  integer(8), allocatable, intent(out) :: grid_index(:), grid_count_ref(:), grid_count_vr(:)
4420  REAL(RP), allocatable, intent(out) :: grid_ref(:), grid_vr(:)
4421  REAL(RP), allocatable, intent(out) :: grid_lon_ref(:), grid_lat_ref(:), grid_z_ref(:) !Lat, Lon and Z weighted by the observations.
4422  REAL(RP), allocatable, intent(out) :: grid_lon_vr(:), grid_lat_vr(:), grid_z_vr(:) !Lat, Lon and Z weighted by the observations.
4423 
4424  REAL(RP), ALLOCATABLE :: grid_w_vr(:)
4425  REAL(RP), ALLOCATABLE :: grid_meanvr(:), grid_stdvr(:) !Non weighted radial velocity and its dispersion within each box.
4426 
4427  integer(8), allocatable :: packed_grid(:), pack2uniq(:), nobs_each_elev(:)
4428  real(RP), allocatable :: packed_data(:, :)
4429  logical, allocatable :: packed_attn(:)
4430 
4431  REAL(RP) :: count_inv, tmpstdvr, tmpmeanvr
4432  integer(8) :: idx, jdx, nobs, sidx(ne), eidx(ne)
4433 
4434  integer :: i
4435  !integer i, e0, e1, ne_mpi
4436  !integer, allocatable :: j_mpi(:, :), sendcount(:), recvcount(:), recvoffset(:)
4437  !integer(8), allocatable :: sendbuf(:), nobs_each_mpi(:)
4438  !integer(8), allocatable :: packed_grid_mpi(:)
4439  !real(RP), allocatable :: packed_data_mpi(:, :)
4440  !logical, allocatable :: packed_attn_mpi(:)
4441 
4442  integer :: ierr
4443 
4444  !allocate(sendcount(0:(mpiprocs - 1)), recvcount(0:(mpiprocs - 1)), recvoffset(0:(mpiprocs - 1)))
4445  !allocate(j_mpi(2, 0:(mpiprocs - 1)))
4446  !call set_mpi_div(j_mpi, mpiprocs, int(ne, 8))
4447  !e0 = j_mpi(1, myrank)
4448  !e1 = j_mpi(2, myrank)
4449  !ne_mpi = e1 - e0 + 1
4450 
4451  ! AVERAGE DATA AND INCLUDE OBSERVATIONA ERROR.
4452  ! We will compute the i,j,k for each radar grid point and box average the
4453  ! data.
4454 
4455  ! QC, Indexing, and packing simultaneously
4456  allocate(nobs_each_elev(ne))
4457  !if(mpiprocs > 1) then
4458  ! allocate(packed_grid_mpi(na * nr * ne_mpi), &
4459  ! & packed_data_mpi(5, na * nr * ne_mpi), &
4460  ! & packed_attn_mpi(na * nr * ne_mpi), &
4461  ! & nobs_each_mpi(e0:e1))
4462  ! call qc_indexing_and_packing( &
4463  ! & na, nr, ne_mpi, ze(:, :, e0:e1), vr(:, :, e0:e1), & ! input spherical
4464  ! & radlon(:, :, e0:e1), radlat(:, :, e0:e1), radz(:, :, e0:e1), & ! input spherical
4465  ! & qcflag(:, :, e0:e1), input_is_dbz, attenuation(:, :, e0:e1), & ! input spherical
4466  ! & nlon, nlat, nlev, lon, lat, z, dlon, dlat, dz, & ! input cartesian
4467  ! & missing, & ! input param
4468  ! & lon0, lat0, &
4469  ! & nobs_each_mpi, packed_grid_mpi, packed_data_mpi, packed_attn_mpi) ! output
4470 
4471  ! sendcount(:) = ne_mpi
4472  ! recvcount(0:(mpiprocs - 1)) = j_mpi(2, 0:(mpiprocs - 1)) - j_mpi(1, 0:(mpiprocs - 1)) + 1
4473  ! recvoffset = j_mpi(1, :) - 1 !START FROM 0
4474  ! call MPI_Allgatherv(nobs_each_mpi, sendcount(0), MPI_INTEGER8, &
4475  ! & nobs_each_elev, recvcount, recvoffset, MPI_INTEGER8, &
4476  ! & comm, ierr)
4477  ! deallocate(nobs_each_mpi)
4478  !else
4479  allocate(packed_grid(na * nr * ne), &
4480  & packed_data(5, na * nr * ne), &
4481  & packed_attn(na * nr * ne))
4482  call qc_indexing_and_packing( &
4483  & na, nr, ne, ze, vr, radlon, radlat, radz, & ! input spherical
4484  & qcflag, input_is_dbz, attenuation, & ! input spherical
4485  & nlon, nlat, nlev, lon, lat, z, dlon, dlat, dz, & ! input cartesian
4486  & missing, & ! input param
4487  & lon0, lat0, &
4488  & nobs_each_elev, packed_grid, packed_data, packed_attn) ! output
4489  !end if
4490  nobs = sum(nobs_each_elev)
4491  sidx(1) = 1
4492  do i = 1, ne - 1
4493  sidx(i + 1) = sidx(i) + nobs_each_elev(i)
4494  end do !i
4495  eidx = sidx + nobs_each_elev - 1
4496  deallocate(nobs_each_elev)
4497  !! MPI packed data
4498  !if(mpiprocs > 1) then
4499  ! sendcount = eidx(e1) - sidx(e0) + 1
4500  ! do i = 0, mpiprocs - 1
4501  ! recvoffset(i) = sidx(j_mpi(1, i)) - 1 !START FROM 0
4502  ! recvcount(i) = eidx(j_mpi(2, i)) - sidx(j_mpi(1, i)) + 1
4503  ! end do
4504  ! allocate(packed_grid(nobs), packed_data(5, nobs), packed_attn(nobs))
4505  ! ! NEEDED BY ALL
4506  ! call MPI_Allgatherv(packed_grid_mpi, sendcount(0), MPI_INTEGER8, &
4507  ! & packed_grid, recvcount, recvoffset, MPI_INTEGER8, &
4508  ! & comm, ierr)
4509  ! ! ONLY NEEDED BY ROOT
4510  ! call MPI_Gatherv(packed_attn_mpi, sendcount(0), MPI_LOGICAL, &
4511  ! & packed_attn, recvcount, recvoffset, MPI_LOGICAL, &
4512  ! & 0, comm, ierr)
4513  ! call MPI_Gatherv(packed_data_mpi, sendcount(0) * 5, datatype, &
4514  ! & packed_data, recvcount * 5, recvoffset * 5, datatype, &
4515  ! & 0, comm, ierr)
4516  ! deallocate(packed_grid_mpi, packed_data_mpi, packed_attn_mpi)
4517 
4518  !end if
4519 
4520  ! Sort index array
4521  allocate(grid_index(nobs))
4522  do idx = 1, nobs
4523  grid_index(idx) = packed_grid(idx)
4524  end do
4525  !if(mpiprocs > 1) then
4526  ! call merge_sort_mpi(nobs, grid_index, comm) !ONLY RANK 0 RETURNS DATA
4527  ! call MPI_Bcast(grid_index, int(nobs), MPI_INTEGER8, 0, comm, ierr)
4528  !else
4529  call merge_sort_parallel(nobs, grid_index)
4530  !end if
4531 
4532  ! Unique superobs (nobs_sp)
4533  call sort_uniq_int_sorted(nobs, grid_index, nobs_sp)
4534 
4535  ! Inverted indexing
4536  allocate(pack2uniq(nobs))
4537  !call set_mpi_div(j_mpi, mpiprocs, nobs)
4538  !allocate(sendbuf(j_mpi(2, myrank) - j_mpi(1, myrank) + 1))
4539  !do idx = j_mpi(1, myrank), j_mpi(2, myrank)
4540  ! sendbuf(idx - j_mpi(1, myrank) + 1) = binary_search_i8(nobs_sp, grid_index, packed_grid(idx))
4541  do idx = 1, nobs
4542  pack2uniq(idx) = binary_search_i8(nobs_sp, grid_index, packed_grid(idx))
4543  end do
4544  !sendcount(:) = j_mpi(2, myrank) - j_mpi(1, myrank) + 1
4545  !recvcount(0:(mpiprocs - 1)) = j_mpi(2, 0:(mpiprocs - 1)) - j_mpi(1, 0:(mpiprocs - 1)) + 1
4546  !recvoffset = j_mpi(1, :) - 1
4547  !! ONLY NEEDED BY ROOT
4548  !call MPI_Gatherv(sendbuf, sendcount(0), MPI_INTEGER8, &
4549  ! & pack2uniq, recvcount, recvoffset, MPI_INTEGER8, &
4550  ! & 0, comm, ierr)
4551  !deallocate(j_mpi, sendbuf, sendcount, recvcount, recvoffset) !END OF MPI
4552 
4553  ! Allocate output arrays
4554  allocate(grid_ref(nobs_sp), grid_vr(nobs_sp))
4555  allocate(grid_count_ref(nobs_sp), grid_count_vr(nobs_sp))
4556  allocate(grid_lon_ref(nobs_sp), grid_lat_ref(nobs_sp), grid_z_ref(nobs_sp))
4557  allocate(grid_lon_vr(nobs_sp) , grid_lat_vr(nobs_sp) , grid_z_vr(nobs_sp))
4558  allocate(grid_w_vr(nobs_sp))
4559  allocate(grid_meanvr(nobs_sp), grid_stdvr(nobs_sp))
4560 
4561  !if(myrank > 0) return !ONLY RANK 0 DOES THE REMAINING WORK
4562 
4563  !Initialize arrays
4564  do idx = 1, nobs_sp
4565  grid_count_ref(idx) = 0
4566  grid_count_vr(idx) = 0
4567  grid_ref(idx) = 0.0d0
4568  grid_vr(idx) = 0.0d0
4569  grid_w_vr(idx) = 0.0d0
4570  grid_lon_ref(idx) = 0.0d0
4571  grid_lat_ref(idx) = 0.0d0
4572  grid_z_ref(idx) = 0.0d0
4573  grid_lon_vr(idx) = 0.0d0
4574  grid_lat_vr(idx) = 0.0d0
4575  grid_z_vr(idx) = 0.0d0
4576  grid_meanvr(idx) = 0.0d0
4577  grid_stdvr(idx) = 0.0d0
4578  end do
4579 
4580  do jdx = 1, nobs
4581  idx = pack2uniq(jdx)
4582  ! PROCESS REFLECITIVITY
4583  ! use attenuation estimates / ignore estimates
4584  if(packed_attn(jdx)) then
4585  grid_ref(idx) = grid_ref(idx) + packed_data(1, jdx)
4586  grid_count_ref(idx) = grid_count_ref(idx) + 1
4587  grid_lon_ref(idx) = grid_lon_ref(idx) + packed_data(3, jdx)
4588  grid_lat_ref(idx) = grid_lat_ref(idx) + packed_data(4, jdx)
4589  grid_z_ref(idx) = grid_z_ref(idx) + packed_data(5, jdx)
4590  ENDIF
4591 
4592  ! CONSIDER THE RADIAL VELOCITY
4593  ! Wind will be averaged using an average weighted by returned power.
4594  ! (this should reduce noise).
4595  IF(packed_data(2, jdx) .GT. missing) THEN !PROCESS WIND
4596  grid_w_vr(idx) = grid_w_vr(idx) + packed_data(1, jdx)
4597  grid_count_vr(idx) = grid_count_vr(idx) + 1
4598  grid_meanvr(idx) = grid_meanvr(idx) + packed_data(2, jdx)
4599  grid_stdvr(idx) = grid_stdvr(idx) + packed_data(2, jdx) ** 2
4600  grid_vr(idx) = grid_vr(idx) + packed_data(2, jdx) * packed_data(1, jdx)
4601  grid_lon_vr(idx) = grid_lon_vr(idx) + packed_data(3, jdx) * packed_data(1, jdx)
4602  grid_lat_vr(idx) = grid_lat_vr(idx) + packed_data(4, jdx) * packed_data(1, jdx)
4603  grid_z_vr(idx) = grid_z_vr(idx) + packed_data(5, jdx) * packed_data(1, jdx)
4604  ENDIF
4605  ENDDO !jdx
4606 
4607  ! Average data and write observation file (FOR LETKF)
4608  DO idx = 1, nobs_sp
4609  IF(grid_count_ref(idx) > 0) THEN !Process reflectivity
4610  count_inv = 1.0d0 / dble(grid_count_ref(idx))
4611  grid_ref(idx) = grid_ref(idx) * count_inv
4612  grid_lon_ref(idx) = grid_lon_ref(idx) * count_inv
4613  grid_lat_ref(idx) = grid_lat_ref(idx) * count_inv
4614  grid_z_ref(idx) = grid_z_ref(idx) * count_inv
4615  ENDIF
4616 
4617  IF(grid_count_vr(idx) > 0) THEN
4618  count_inv = 1.0d0 / grid_w_vr(idx)
4619  grid_vr(idx) = grid_vr(idx) * count_inv
4620  grid_lon_vr(idx) = grid_lon_vr(idx) * count_inv
4621  grid_lat_vr(idx) = grid_lat_vr(idx) * count_inv
4622  grid_z_vr(idx) = grid_z_vr(idx) * count_inv
4623 
4624  ! If variability within a box is big then we may have:
4625  ! -small scale strong phenomena (tornado!)
4626  ! -noise in the data.
4627  ! In any case averaging the data is not a god idea so this data
4628  ! can be rejected a priori.
4629  IF( radar_use_vr_std ) THEN
4630  count_inv = 1.0d0 / dble(grid_count_vr(idx))
4631  tmpmeanvr = grid_meanvr(idx) * count_inv
4632  tmpstdvr = grid_stdvr(idx) * count_inv
4633  tmpstdvr = sqrt(tmpstdvr - (tmpmeanvr ** 2))
4634  IF(tmpstdvr > vr_std_threshold) grid_count_vr(idx) = 0 !Reject the observation.
4635  ENDIF
4636  end IF
4637  end do
4638 
4639  return

References com_gamma(), scale_const::const_d2r, scale_const::const_grav, scale_const::const_pi, scale_const::const_r2d, scale_const::const_radius, scale_const::const_rdry, and scale_sort::sort_uniq_int_sorted().

Referenced by letkf_add_inflation_setup().

Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ nv3d

integer, parameter, public scale_letkf::nv3d = 11

Definition at line 47 of file scale_letkf.F90.

47  integer, public, parameter :: nv3d = 11 ! number of 3D prognostic variables

Referenced by letkf_add_inflation_setup(), letkf_setup(), and letkf_system().

◆ nv2d

integer, parameter, public scale_letkf::nv2d = 0

Definition at line 48 of file scale_letkf.F90.

48  integer, public, parameter :: nv2d = 0 ! number of 2D prognostic variables

Referenced by letkf_add_inflation_setup(), letkf_setup(), and letkf_system().

◆ nid_obs

integer, parameter, public scale_letkf::nid_obs = 16

Definition at line 49 of file scale_letkf.F90.

49  integer, public, parameter :: nid_obs = 16 ! number of variable types

Referenced by letkf_obs_initialize(), and letkf_system().

◆ nobtype

integer, parameter, public scale_letkf::nobtype = 24

Definition at line 50 of file scale_letkf.F90.

50  integer, public, parameter :: nobtype = 24 ! number of observation report types

Referenced by letkf_obs_initialize(), and letkf_system().

◆ max_obs_info_meta

integer, parameter, public scale_letkf::max_obs_info_meta = 3

Definition at line 52 of file scale_letkf.F90.

52  integer, public, parameter :: max_obs_info_meta = 3 ! maximum array size for type(obs_info)%meta

Referenced by letkf_obs_readfile().

◆ n_qc_steps

integer, parameter, public scale_letkf::n_qc_steps = 2

Definition at line 53 of file scale_letkf.F90.

53  integer, public, parameter :: n_qc_steps = 2

Referenced by letkf_obs_initialize().

◆ i_before_qc

integer, parameter, public scale_letkf::i_before_qc = 1

Definition at line 54 of file scale_letkf.F90.

54  integer, public, parameter :: i_before_qc = 1

Referenced by letkf_obs_initialize().

◆ i_after_qc

integer, parameter, public scale_letkf::i_after_qc = 2

Definition at line 55 of file scale_letkf.F90.

55  integer, public, parameter :: i_after_qc = 2

Referenced by letkf_obs_initialize().

◆ n_search_incr

integer, parameter scale_letkf::n_search_incr = 8

Definition at line 183 of file scale_letkf.F90.

183  integer, parameter :: n_search_incr = 8

Referenced by com_gamma().

scale_statistics
module Statistics
Definition: scale_statistics.F90:11
scale_precision::sp
integer, parameter, public sp
Definition: scale_precision.F90:31
scale_atmos_grid_cartesc_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:52
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:350
scale_atmos_grid_cartesc_index::ihalo
integer, public ihalo
Definition: scale_atmos_grid_cartesC_index.F90:44
scale_atmos_grid_cartesc_index::ka
integer, public ka
Definition: scale_atmos_grid_cartesC_index.F90:47
scale_sort
module SORT
Definition: scale_sort.F90:11
scale_random::random_knuth_shuffle
subroutine, public random_knuth_shuffle(num, a)
Definition: scale_random.F90:240
scale_random
module RANDOM
Definition: scale_random.F90:11
scale_atmos_grid_cartesc_index::khalo
integer, parameter, public khalo
Definition: scale_atmos_grid_cartesC_index.F90:43
scale_atmos_grid_cartesc_index::imax
integer, public imax
Definition: scale_atmos_grid_cartesC_index.F90:37
scale_matrix::matrix_solver_eigenvalue_decomposition
subroutine, public matrix_solver_eigenvalue_decomposition(n, a, eival, eivec, simdlen)
Definition: scale_matrix.F90:879
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_atmos_grid_cartesc_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:54
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_grid_cartesc_index::ia
integer, public ia
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_prc_cartesc::prc_2drank
integer, dimension(:,:), allocatable, public prc_2drank
node index in 2D topology
Definition: scale_prc_cartesC.F90:45
scale_atmos_grid_cartesc_index::kmax
integer, public kmax
Definition: scale_atmos_grid_cartesC_index.F90:36
scale_prc_cartesc
module process / cartesC
Definition: scale_prc_cartesC.F90:11
scale_atmos_grid_cartesc_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:53
scale_atmos_grid_cartesc_index::jhalo
integer, public jhalo
Definition: scale_atmos_grid_cartesC_index.F90:45
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
scale_atmos_grid_cartesc_index::ja
integer, public ja
Definition: scale_atmos_grid_cartesC_index.F90:49
scale_time
module TIME
Definition: scale_time.F90:11
scale_atmos_grid_cartesc_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:51
scale_const::const_tem00
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
Definition: scale_const.F90:99
scale_time::time_gettimelabel
subroutine, public time_gettimelabel(timelabel)
generate time label
Definition: scale_time.F90:93
scale_atmos_grid_cartesc_index::jmax
integer, public jmax
Definition: scale_atmos_grid_cartesC_index.F90:38
scale_atmos_grid_cartesc_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:55
scale_matrix
module MATRIX
Definition: scale_matrix.F90:17
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
scale_atmos_grid_cartesc_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:56
scale_sort::sort_uniq_int_sorted
subroutine, public sort_uniq_int_sorted(n, ary, c)
Definition: scale_sort.F90:200
Generated by   doxygen 1.8.17