text:
! This Fortran95 Code/Module performs Machine Learning Force prediction based on XYZ structural Files;
! Sorting/Selecting Algorithms are used to construct sub training database during heavy machine learning task;
! Database can be updated according to an error_threshold, e.g. pred_err;
! Internal Vector representation is used for describing the atomic structures;
! Any issues or bugs regarding this code please contact Zhenwei Li (zhenwei.li@unibas.ch)
! Reference article : Molecular Dynamics with On-the-Fly Machine Learning of Quantum-Mechanical Forces, Zhenwei Li, James Kermode, Alessandro De Vita, PRL, 114, 096405 (2015)
program force_gaussian_prediction

   use potential_module
   use libatoms_module
   use force_machine_learning_module

   implicit none

   type(Atoms)         :: at_in
   type(Potential)     :: pot
   type(CInOutput)     :: in
   type(Dictionary)    :: params
   real(dp)            :: r_cut, feature_len, thresh, sigma_error, error_ls, error_gp, error_gp_sq, time, error_ls_max
   real(dp)            ::  sigma_covariance, sigma_cov, dist_shift_factor, dim_tol, error_frame, force_magntd_pred, dist
   real(dp), dimension(:), allocatable           :: r_grid, m_grid, sigma, covariance_pred,force_magntd_data,force_magntd_data_select
   real(dp), dimension(:), allocatable           :: force_proj_test, force_proj_ivs_pred, distance_confs_unsorted
   real(dp), parameter                           :: TOL_REAL=1e-7_dp, SCALE_IVS=100.0_dp
   real(dp)                                      :: force(3), variance_xyz(3), feature_inner_matrix(3,3), feature_inv(3,3), kappa, tmp
   real(dp), dimension(:,:,:), allocatable       :: ivs_set_norm, ivs_set, ivs_set_pred, nearest_data
   real(dp), dimension(:,:), allocatable         :: force_proj_ivs, force_proj_ivs_select, inv_covariance, out_u, out_vt, geometry_matrix
   real(dp), dimension(:,:), allocatable, target :: covariance, distance_matrix, distance_matrix_dij, force_covariance_matrix, covariance_tiny
   real(dp), dimension(:,:), allocatable         :: ivs_set_norm_pred, ivs_set_norm_pred_t
   real(dp), dimension(:,:), pointer             :: force_ptr, force_ptr_mm, force_ptr_tff, pred_err
   integer                                       :: i,j, k, n, n_0, t,  n_center_atom, n_loop, error
   integer                                       :: add_vector, n_teach, func_type, weight_index, temp_integer, ii, do_atz_index
   integer	                 		 :: local_ml_optim_size, dimension_mark(3)
   integer, dimension(:), allocatable            :: distance_index, mark_zero_ivs, species
   logical                         ::  top_teach, top_pred, do_gp, do_sigma, least_sq, print_at, do_svd, print_data,do_teach, do_read, print_verbosity, update_data, do_insert, do_dynamic_x, write_data_dij, read_data_dij
   character(STRING_LENGTH)                      :: teach_file, test_file, out_file, data_dir

   call system_initialise(enable_timing=.true.)

   call initialise(params)
   call param_register(params, 'r_cut',  '8.0', r_cut, "the cutoff radius for the spherical atomic environment")
   call param_register(params, 'thresh', '10.0', thresh, "the threshold for doing the Singular Value Decompostion of the Covariance Matrix")
   call param_register(params, 'weight_index', '0', weight_index, "Type of weight function in deriving IVs, zero stands for the form of $exp(r/rcut)^m$")
   call param_register(params, 'sigma_error', '0.05', sigma_error, "the noise assumed on the teaching data")
   call param_register(params, 'do_atz_index', '0', do_atz_index, "0: charge-weight disabled; 1: using at%Z, 2: using ionic charge instead")
   call param_register(params, 'func_type', '0', func_type, "which kernel function is used to build the covariance matrix")
   call param_register(params, 'do_gp',  'T', do_gp, "true for doing a gaussian processes")
   call param_register(params, 'n_teach', '1000', n_teach, "the number of relevant confs you would like to do machine learning with")
   call param_register(params, 'do_sigma', 'F', do_sigma, "set true to print out the distance on every single dimension of the IVs space")
   call param_register(params, 'do_dynamic_x', 'F', do_dynamic_x, "true for updating sigma.dat every test configuration based on nearest data points")
   call param_register(params, 'top_teach', 'T', top_teach, "only the first configuration is considered when doing teaching")
   call param_register(params, 'top_pred', 'T', top_pred, "only the first configuration is considered when doing predicting")
   call param_register(params, 'least_sq',    'T', least_sq, "if true, the internal force components will be tranformed to real force using least squares")
   call param_register(params, 'do_svd', 'F', do_svd, "if true, doing inverting the feature_inner_matrix by SVD")
   call param_register(params, 'sigma_cov', '1.0', sigma_cov, "correlation length for pairs of data points")
   call param_register(params, 'dist_shift_factor', '0.2', dist_shift_factor, "distance shifted with a given factor")
   call param_register(params, 'dim_tol', '0.005', dim_tol, "threshold for determing the dimensionality of the group of IVs")
   call param_register(params, 'teach_file', 'data.xyz', teach_file, "file to read teaching configurations from")
   call param_register(params, 'test_file', 'test.xyz', test_file, "file to read the testing configurations from")
   call param_register(params, 'print_data', 'F', print_data, "if true, print out the information for the database" )
   call param_register(params, 'data_dir', './info/', data_dir, "the directory where data files locate")
   call param_register(params, 'do_teach', 'T', do_teach, "if false, the teaching part will be skipped")
   call param_register(params, 'do_read', 'F', do_read, "if true, reading data file from the database")
   call param_register(params, 'do_insert', 'F', do_insert, "INsert any addiitional vector as force_plus")
   call param_register(params, 'print_verbosity', 'F', print_verbosity, "if true, print out verbosity stuff")
   call param_register(params, 'out_file', 'out.xyz', out_file, "output configuration file containing the predicted forces")
   call param_register(params, 'print_at', 'F', print_at, "true for print out the testing configuration with predicted forces")
   call param_register(params, 'update_data', 'F', update_data, "whether or not updating the database")
   call param_register(params, 'read_data_dij', 'F', read_data_dij, 'whether or not writing the distance file')
   call param_register(params, 'write_data_dij', 'F', write_data_dij, 'whether or not reading the pair distnace matrix')
   call param_register(params, 'local_ml_optim_size', '0', local_ml_optim_size, "Optimise sigma_error and sigma_cov using local_ml_optim_size LS confs. If 0, no optimisation is performed")

   if (.not. param_read_args(params, task="fgp command line params")) then
      call print("Usage: fgp [options]")
      call system_abort("Confused by command line arguments")
   end if

   call print_title('params')
   call param_print(params)
   call finalise(params)

if (do_teach) then
  !
  call load_iv_params(data_dir, r_grid, m_grid, species, k, add_vector, n_data=n_0)
  ! the current k is number of defined internal vectors
  if (add_vector > 0)   k=k+add_vector
  call print('The number of valid Feature Vectors (inc. f_mm): '//k)
  allocate(sigma(k)) ! knowing k, sigma canbe assigned

  ! to know the total number of teaching confs: n
  call initialise(in, teach_file, INPUT)
  n = 0
  do i=1,  in%n_frame
        call read(in, at_in)
        if (top_teach) then
            n = n + 1
        else
            n = n +  at_in%N
        endif
  enddo
  call finalise(in)
  call print("The total number of teaching configurations:    "//n)

  call print_title('starting the teaching process')
  call system_timer('Preparing Information from the DATABASE')

  ! the main work starts
  call initialise(in, teach_file, INPUT)
  call print('Got '//in%n_frame//' input frames.')

  allocate(ivs_set(k,3,n))
  allocate(ivs_set_norm(k,3,n))
  allocate(force_proj_ivs(k, n))
  allocate(force_magntd_data(n))

 t=0 ! counter for loop the atomic configurations
 do i=1, in%n_frame
        call read(in, at_in)
        call assign_property_pointer(at_in, 'force', ptr=force_ptr)
        if (add_vector > 0) then
           if (do_insert) then
              call assign_property_pointer(at_in, 'force_plus', ptr=force_ptr_mm) !treat the plus Force as mm
           else
             call Potential_Filename_Initialise(pot, args_str='IP SW', param_filename='SW.xml')
             call calc(pot, at_in, args_str='force=force_mm',  error=error)
             call assign_property_pointer(at_in, 'force_mm', ptr=force_ptr_mm)
             call finalise(pot)
           endif
        endif
        ! the first vector is labelled force_mm in the input file, and is usually the SW force
        if (add_vector > 1) call assign_property_pointer(at_in, 'force_i', ptr=force_ptr_tff)
        ! the second vector is labelled force_i
        ! call print('the frame: '//i)

        ! Including the symmetry images
        call set_cutoff(at_in, r_cut)
        call calc_connect(at_in)

     if (top_teach) then
        n_loop = 1
     else
        n_loop = at_in%N
     endif

    do n_center_atom =1,  n_loop
          t = t + 1
          force = force_ptr(:,n_center_atom)
          !call print("atomic force in real space : "//norm(force))
          !call print("The atom:"//n_center_atom)

        do j=1, k-add_vector
           ivs_set(j,:,t) = internal_vector(at_in, r_grid(j), m_grid(j), species(j), n_center_atom, weight_index, do_atz_index) *SCALE_IVS
           feature_len = norm(ivs_set(j,:, t))
           if (feature_len < TOL_REAL)  then
                   feature_len=1.0_dp
                   call print("WARNING: TEACHING, encountered the numerical limit in getting the unit direction of IVs")
           endif

           ivs_set_norm(j,:,t) = ivs_set(j,:, t)/feature_len
           if (print_verbosity) call print("internal vectors ( "//r_grid(j)//" "//m_grid(j)//" "//species(j)//" )"//ivs_set(j,1,t)//" "//ivs_set(j,2,t)//" "//ivs_set(j,3,t))
        enddo

        if (add_vector > 0 ) then
           do j=k-add_vector+1, k
              temp_integer=k-j
              select case (temp_integer)
              case (0)
                 ivs_set(j,:,t) = force_ptr_mm(:,n_center_atom)
              case (1)
                 ivs_set(j,:,t) = force_ptr_tff(:,n_center_atom)
              end select

              feature_len = norm(ivs_set(j,:,t))
              if (feature_len < TOL_REAL)  then
                 feature_len=1.0_dp
                 call print("WARNING: TEACHING, encountered the numerical limit in getting the unit direction of IVs")
              endif

              ivs_set_norm(j,:,t) = ivs_set(j,:,t)/feature_len
              !call print("added internal vectors :"//ivs_set(j,1,t)//"  "//ivs_set(j,2,t)//" "//ivs_set(j,3,t))
           enddo
        endif

        ! store the force info for every configuration
        force_magntd_data(t) = norm(force)
        do j=1, k
             force_proj_ivs(j,t) = dot_product(ivs_set_norm(j,:,t), force)
        enddo

     enddo   ! loop over the atoms
enddo       ! loop over the frames

call finalise(in)
call system_timer('Preparing Information from the DATABASE')

if (do_sigma) then
   call sigma_factor_from_single_dimension(ivs_set(:,:,:), sigma, print_verbosity, TOL_REAL, data_dir)
endif

!---------------------------------------------------------------------------------------------------------------------------------
!  	WRITTING OUT THE DATAFILES
!---------------------------------------------------------------------------------------------------------------------------------

if (print_data) then
  call system_timer('Writing Data File: Grid, Force and IV')
  ! for updating database, initial value taken into account
  call  write_iv_params(data_dir, r_grid, m_grid, species, k-add_vector, add_vector, n)
  if (.not. update_data ) call write_distance(0.0_dp, data_dir)
  !
  OPEN(2,status='replace',file=trim(data_dir)//'Force.dat', form='UNFORMATTED')
  OPEN(3,status='replace',file=trim(data_dir)//'IV.dat', form='UNFORMATTED')

	do t=1, n
	    write(2)  force_proj_ivs(:, t), force_magntd_data(t)
	    do i=1,3
	      write(3)  ivs_set(:,i,t)
	    enddo
	enddo ! for n configuration to augument the existing DATABASE
        close(2)
        close(3)
        call system_timer('Writing Data File: Grid, Force and IV')

endif  ! print_data or not
end if ! do_teach, otherwise skipped

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! if do_read, teaching information to be collected from data files
if (do_read) then

    call print_title('READING TEACHING INFORMATION')
    call system_timer('Collecting Information From Files')
    call load_iv_params(data_dir, r_grid, m_grid, species, k, add_vector, n_data=n)
    k = k + add_vector  ! adding the Classical Force Vectors

    allocate(force_proj_ivs(k,n))
    allocate(ivs_set(k,3,n))
    allocate(ivs_set_norm(k,3,n))
    allocate(force_magntd_data(n))
    allocate(sigma(k))

    OPEN(2,status='old',file=trim(data_dir)//'Force.dat',form='UNFORMATTED')
    OPEN(3,status='old',file=trim(data_dir)//'IV.dat',form='UNFORMATTED')

    do t=1, n
       read(2) force_proj_ivs(:, t), force_magntd_data(t)
       !write(*,*) "%", force_proj_ivs(:, t), force_magntd_data(t)
       do i=1, 3
          read(3) ivs_set(:,i,t)
       enddo
       !write(*,*) "$", ivs_set(:,:,t)
    enddo
    close(2)
    close(3)

  ! calculating the normalised vectors
  do t=1, n
    do j=1, k
     feature_len = norm(ivs_set(j,:, t))
     if (feature_len < TOL_REAL)  then
             feature_len=1.0_dp
             call print("WARNING: TEACHING, encountered the numerical limit in getting the unit direction of IVs")
     endif
     ivs_set_norm(j,:,t)=ivs_set(j,:, t)/feature_len
    enddo
 enddo

    ! reading sigma file
    open(unit=1, status='old', file=trim(data_dir)//'sigma.dat', form='formatted')
    read(1, *) sigma
    close(unit=1)
    call system_timer('Collecting Information From Files')
    call print('sigma is:    '//sigma)
endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

if ((do_teach) .and. (.not. do_read)) then
    open(unit=1, status='old', file=trim(data_dir)//'sigma.dat', form='formatted')
    read(1, *) sigma
    close(unit=1)
    call system_timer('Collecting Information From Files')
    call print('sigma is:    '//sigma)
endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! the following part can be changed to a subroutine
call print_title('starting the predicting process')
!
if (n_teach > n) then
  n_teach = n
  call print("the Actual Number of Configs for Prediction: "//n_teach)
  call print("WARNNING: Trying to select more configurations than the database has")
endif
!
allocate(covariance_pred(n_teach))
allocate(force_proj_ivs_pred(k))
allocate(mark_zero_ivs(k))
allocate(geometry_matrix(3,k))
allocate(covariance(n_teach,n_teach))
allocate(distance_confs_unsorted(n), distance_index(n))
if (local_ml_optim_size > 0) then
allocate(distance_matrix(local_ml_optim_size,local_ml_optim_size), force_covariance_matrix(local_ml_optim_size,local_ml_optim_size), covariance_tiny(local_ml_optim_size,local_ml_optim_size))
endif
allocate(inv_covariance(n_teach,n_teach))
allocate(force_proj_ivs_select(k, n_teach))
allocate(force_magntd_data_select(n_teach))
allocate(ivs_set_norm_pred(k,3))
if (read_data_dij) allocate(distance_matrix_dij(n, n))

if (read_data_dij) then
          call system_timer('Read Pair Distance and Bulid Covariance')
          call read_distance_file(distance_matrix_dij(:,:), data_dir)
          !write(*,*)  'Distance Read', distance_matrix_dij
          call system_timer('Read Pair Distance and Bulid Covariance')
endif

call initialise(in, test_file, INPUT)
do i=1, in%n_frame
   call read(in, at_in)
   if (.not. get_value(at_in%params, 'time', time))  then
     call print("TIME of FRAME : UNKOWN")
   else
     call print("TIME of FRAME :"//time)
   endif
   call add_property(at_in, 'force', 0.0_dp, n_cols=3, overwrite=.false.) ! how do you know dont want overwrite
   call assign_property_pointer(at_in, 'force', ptr=force_ptr)

   ! for writing out the predicted error
   call add_property(at_in, 'pred_err', 0.0_dp, n_cols=3, overwrite=.true.)
   call assign_property_pointer(at_in, 'pred_err', ptr=pred_err)

   if (add_vector > 0) then
           if (do_insert) then
              call assign_property_pointer(at_in, 'force_plus', ptr=force_ptr_mm)
           else
             call Potential_Filename_Initialise(pot, args_str='IP SW', param_filename='SW.xml')
             call calc(pot, at_in, args_str='force=force_mm',  error=error)
             call assign_property_pointer(at_in, 'force_mm', ptr=force_ptr_mm)
             call finalise(pot)
           endif
   endif
   if (add_vector>1) call assign_property_pointer(at_in, 'force_i', ptr=force_ptr_tff)

   call set_cutoff(at_in, r_cut)
   call calc_connect(at_in)
   if (top_pred) then
      n_loop = 1
   else
      n_loop = at_in%N
   endif

   allocate(ivs_set_pred(k,3,n_loop))

   error_frame=0.0_dp           ! for calculating the average error of the system.

  do n_center_atom=1, n_loop    ! loop over atoms in the Frame, make prediction (incl. error)

     do j=1, k                  ! reset the mark for zero internal vectors
         mark_zero_ivs(j)=1
     enddo

     do j= 1, k-add_vector
        ivs_set_pred(j,:,n_center_atom) = internal_vector(at_in, r_grid(j), m_grid(j), species(j), n_center_atom, weight_index, do_atz_index) *SCALE_IVS
        if (print_verbosity) call print("internal vectors ( "//r_grid(j)//" "//m_grid(j)//" "//species(j)//" ):  "//ivs_set_pred(j,1,n_center_atom)//"  "//ivs_set_pred(j,2,n_center_atom)//"  "//ivs_set_pred(j,3,n_center_atom))
       feature_len = norm(ivs_set_pred(j,:,n_center_atom))

        if (feature_len < TOL_REAL)  then
            feature_len=1.0_dp
            mark_zero_ivs(j)=0
            call print("WARNING: PREDICTION, encountered the numerical limit in getting the unit direction of IVs")
        endif

        ivs_set_norm_pred(j,:) = ivs_set_pred(j,:,n_center_atom)/feature_len
     enddo

     if (add_vector>0) then
         do j = k-add_vector+1, k
            temp_integer=k-j
            select case (temp_integer)
            case (0)
               ivs_set_pred(j,:,n_center_atom)=force_ptr_mm(:, n_center_atom)
            case (1)
               ivs_set_pred(j,:,n_center_atom)=force_ptr_tff(:, n_center_atom)
            end select
            feature_len = norm(ivs_set_pred(j,:,n_center_atom))

            if (feature_len < TOL_REAL) then
               feature_len=1.0_dp
               mark_zero_ivs(j)=0
               call print("WARNING: PREDICTION, encountered the numerical limit in getting the unit direction of IVs")
            endif

            ivs_set_norm_pred(j,:) = ivs_set_pred(j,:,n_center_atom)/feature_len
            if (print_verbosity) call print("added internal vectors : "//ivs_set_pred(j,1,n_center_atom)//"  "//ivs_set_pred(j,2,n_center_atom)//"  "//ivs_set_pred(j,3,n_center_atom))
         enddo
      endif

     if (print_verbosity) then
       do j=1, k
          do ii=1, k
           if(print_verbosity) then
              call print("ATOM : "//n_center_atom//" ivs_set: ("//j//" "//ii//") "//dot_product(ivs_set_pred(j,:,n_center_atom), ivs_set_norm_pred(ii,:)) )
           endif
          enddo
       enddo
     endif

      ! initialise distance_index(:)
      do t=1, n
         distance_index(t) =t
      enddo

        ! do the sorting and selection
      call sorting_configuration(ivs_set_pred(:,:,n_center_atom), ivs_set, ivs_set_norm_pred, ivs_set_norm, sigma, distance_confs_unsorted, distance_index)
      call print("The First 20 Closest data Configurations are :"//distance_index(:20))
      call print("Min and Max DISTANCE with Index after Sorting: "//distance_confs_unsorted(distance_index(1))//" and "// &
	                                    distance_confs_unsorted(distance_index(n_teach))//"  the INDEX: "//distance_index(1)//" and "//distance_index(n_teach))
       !!!!!!!!!!!!!
    if (do_dynamic_x) then
        allocate(nearest_data(k,3,n_teach))
        do ii=1, n_teach
           nearest_data(:,:,ii)=ivs_set(:,:,distance_index(ii))
        enddo
        call sigma_factor_from_single_dimension(nearest_data, sigma, print_verbosity, TOL_REAL, data_dir)
        ! update 'sigma.dat'
        deallocate(nearest_data)
     endif

     if ( distance_confs_unsorted( distance_index(1) ) > dist_shift_factor) then
                ! distance_confs_unsorted: the pair distance betwen test and database.
                ! sigma_cov: the original input for covariance length
                sigma_covariance = sigma_cov*distance_confs_unsorted(distance_index(1))/dist_shift_factor
                ! this intrinsically changes the value of SIGMA_COVARIANCE and enters next loop
                call print("sigma_covariance are shifted by a factor of :"//distance_confs_unsorted(distance_index(1))/dist_shift_factor)
                call print("Modified sigma_covariance :"//sigma_covariance)
     else
               sigma_covariance = sigma_cov
               ! the sigma_covariance used in GP
     endif

     if (.false.) then
         ! massive output, only for testing use
         call print("Detail on the teaching information")
         do ii = 1, n_teach
            call print("Index: "//distance_index(ii)//" Distance: "//distance_confs_unsorted(ii))
            call print("Force in IVs space: "//force_proj_ivs(:,distance_index(ii)))
         enddo
     endif

     ! max marginal likelihood
      if (local_ml_optim_size > 0) call print("Optimise: "//local_ml_optim_size)
      !!! ML optimisation here: first calculates the necessary input matrices, then optimises
      if (local_ml_optim_size > 0) then
         call print("Local Maximum Likelihood selected")
         do t = 1, local_ml_optim_size  ! loop to generate the force_cov_mat and the distance_matrix of the restricted subset of the LS
            do j = t, local_ml_optim_size
               force_covariance_matrix(t,j) = dot_product(force_proj_ivs(:,distance_index(t)), force_proj_ivs(:,distance_index(j)))
               ! check if possible, may be a wreck because IVs make absolute representation of the force impossible
               force_covariance_matrix(j,t) = force_covariance_matrix(t,j)
               covariance_tiny(t,j) = cov(ivs_set(:,:,distance_index(t)), ivs_set(:,:,distance_index(j)), &
                    ivs_set_norm(:,:,distance_index(t)), ivs_set_norm(:,:,distance_index(j)), sigma, sigma_covariance, func_type=func_type, &
                    distance=distance_matrix(t,j))
               covariance_tiny(j,t) = covariance_tiny(t,j)
               distance_matrix(j,t) = distance_matrix(t,j)
            enddo
         enddo
         call do_optimise_likelihood(sigma_error, sigma_covariance)
      endif

      if (read_data_dij) then
          do t = 1, n_teach  ! loop to generate the covariance matrix of the learning set
                 dist=distance_matrix_dij(distance_index(t), distance_index(t))
                 covariance(t,t) = cov_dij(sigma_covariance, dist, func_type)
           !      write(*,*) 'diagoanl value', covariance(t,t)
                 if (do_gp) covariance(t,t) = covariance(t,t) + sigma_error**2
                 do j = t+1, n_teach ! above diagonal
                    dist=distance_matrix_dij(distance_index(t),distance_index(j))
                    covariance(t,j) = cov_dij(sigma_covariance, dist, func_type)
                    covariance(j,t) = covariance(t,j)
                 enddo
          enddo
      else
	      call system_timer('Evaluating the Cov')
	      !Generate covariance matrix
	      !$omp parallel default(none) shared(covariance, n_teach, ivs_set, ivs_set_norm, distance_index, sigma, sigma_error, sigma_covariance, do_gp, func_type), private(j)
	      !$omp do
	      do t = 1, n_teach  ! loop to generate the covariance matrix of the learning set
		 covariance(t,t) = cov(ivs_set(:,:,distance_index(t)), ivs_set(:,:,distance_index(t)), &
			ivs_set_norm(:,:,distance_index(t)), ivs_set_norm(:,:,distance_index(t)), sigma, sigma_covariance, func_type=func_type)
		 if (do_gp) covariance(t,t) = covariance(t,t) + sigma_error**2
		 do j = t+1, n_teach ! above diagonal
		    covariance(t,j) = cov(ivs_set(:,:,distance_index(t)), ivs_set(:,:,distance_index(j)), &
			 ivs_set_norm(:,:,distance_index(t)), ivs_set_norm(:,:,distance_index(j)), sigma, sigma_covariance, func_type=func_type)
		    covariance(j,t) = covariance(t,j)
		enddo
	       enddo
	      !$omp end parallel
	       call system_timer('Evaluating the Cov')
      endif ! not read_data_dij

     call system_timer('Inverting the Covariance Matrix')
      if (do_gp) then
         call inverse(covariance, inv_covariance)
      else
         ! To Do Sigular Value Decomposition (SVD): A = U*SIGMA*VT
         call inverse_svd_threshold(covariance, thresh, result_inv=inv_covariance)
      endif
      call system_timer('Inverting the Covariance Matrix')

      do t=1, n_teach
         force_proj_ivs_select(:,t)=force_proj_ivs(:,distance_index(t))
         force_magntd_data_select(t) = force_magntd_data(distance_index(t))
      enddo

      ! the test configuration covariance vector
      do t= 1, n_teach
        covariance_pred(t) = cov(ivs_set_pred(:,:,n_center_atom), ivs_set(:,:,distance_index(t)), ivs_set_norm_pred, &
                                         ivs_set_norm(:,:,distance_index(t)), sigma, sigma_covariance, func_type=func_type)
      enddo

      ! the predicted force components on each of the internal directions
      force_proj_ivs_pred(:) = matmul(covariance_pred, matmul(inv_covariance, transpose(force_proj_ivs_select(:,:)) ))

      ! Predict the Force Magnitude
      force_magntd_pred =  dot_product(covariance_pred, matmul(inv_covariance, force_magntd_data_select))

      ! we need to mannually set the force components to be zero, projected on the zero internal-vector directions
      do j=1, k
         force_proj_ivs_pred(j)=force_proj_ivs_pred(j)*real(mark_zero_ivs(j), dp)
      enddo

     if (print_verbosity) then
         do j=1, k
               call print("ivs_set_norm_pred"//ivs_set_norm_pred(j,:))
         enddo
     endif

     allocate(out_u(3,3), out_vt(3,3))      ! to obtain the transformation matrix with the principal axis
     call inverse_svd_threshold(transpose(ivs_set_norm_pred) .mult. ivs_set_norm_pred, thresh, u_out=out_u, vt_out=out_vt)

     allocate(ivs_set_norm_pred_t(k,3))

     call internal_dimension_mark(ivs_set_pred(:,:,n_center_atom) .mult. out_u, dim_tol, dimension_mark)

     ivs_set_norm_pred_t = ivs_set_norm_pred .mult. out_u

  ! ivs_set_pred_t contains the Internal-directions after PCA transformation.
  if (print_verbosity) then
       do j=1, k
           write(*,*) "ivs_set_norm_pred_t : ", ivs_set_norm_pred_t(j, :)
       enddo
  endif

  select case (sum(dimension_mark(:)))
    case(0)
         force= 0.0_dp
    case(1)
         force=(transpose(ivs_set_norm_pred) .mult. force_proj_ivs_pred )/real(sum(mark_zero_ivs),dp)
    case(2)
         feature_inner_matrix=transpose(ivs_set_norm_pred_t) .mult. ivs_set_norm_pred_t
         call rank_lowered_inverse(feature_inner_matrix, feature_inv)
         force = out_u .mult. feature_inv .mult. transpose(ivs_set_norm_pred_t) .mult. force_proj_ivs_pred
    case(3)

         feature_inner_matrix=transpose(ivs_set_norm_pred) .mult. ivs_set_norm_pred
         if (do_svd) then
             if (print_verbosity) write(*,*)  "feature inner matrix :", feature_inner_matrix
             call inverse_svd_threshold(feature_inner_matrix, thresh, result_inv=feature_inv)
             if (print_verbosity) write(*,*)  "inverse feature inner matrix :", feature_inv
         else
             if (print_verbosity) write(*,*)  "feature_inner_matrix :", feature_inner_matrix
             call inverse(feature_inner_matrix, feature_inv)
             if (print_verbosity) write(*,*)  "inverse feature inner matrix :", feature_inv
         endif

         if (least_sq) then
            force = feature_inv .mult. transpose(ivs_set_norm_pred) .mult. force_proj_ivs_pred
         else
            call real_expection_sampling_components(ivs_set_norm_pred, force_proj_ivs_pred, force)
         endif
   end select

  deallocate(out_u)
  deallocate(out_vt)
  deallocate(ivs_set_norm_pred_t)

   call print("Predicted Force Magnitude: "//force_magntd_pred)
   call print("force in external space : "//force//" Magnitude: "//norm(force))
   call print("the original force:"//force_ptr(:, n_center_atom)//" Magnitude: "//norm(force_ptr(:, n_center_atom)))
   call print("max error :    "//maxval(abs(force_ptr(:,n_center_atom)-force))//" norm  error :  "//norm(force_ptr(:,n_center_atom)-force))
   if (norm(force_ptr(:,n_center_atom)) > TOL_REAL) then
           call print("Relative Error : "//norm(force_ptr(:,n_center_atom)-force)/norm(force_ptr(:,n_center_atom)) )
   else
           call print("Relative Error : "//" ZERO Magnititude ")
   endif
   ! measure the consistence between the least-square inference and the components obtained from GP processes
   if (least_sq) then
       error_ls=0.0_dp
       do j=1, k
            error_ls = error_ls + (dot_product(ivs_set_norm_pred(j,:), force(:)) - force_proj_ivs_pred(j))**2
       enddo
       error_ls = sqrt(error_ls / real(k,dp))
       call print("deviation of least-square inference : "//error_ls)
   error_ls_max=maxval(abs( matmul(ivs_set_norm_pred, force) - force_proj_ivs_pred(:) ))
   call print("Max deviation Component of LS : "//error_ls_max)
   endif

   kappa = cov(ivs_set_pred(:,:,n_center_atom), ivs_set_pred(:,:,n_center_atom), ivs_set_norm_pred, ivs_set_norm_pred, sigma, sigma_covariance, func_type=func_type) + sigma_error**2
   ! the passing of  uncertainty from  Gaussian Processes to the Force Vector
   error_gp_sq = kappa - covariance_pred .dot. matmul(inv_covariance, covariance_pred)
   error_gp = sqrt(abs(error_gp_sq))
   geometry_matrix = feature_inv .mult. transpose(ivs_set_norm_pred)

   do j=1, 3  ! Deviation in External Cartesian Coordinate
      variance_xyz(j)= sqrt(dot_product(geometry_matrix(j,:), geometry_matrix(j,:))) * (error_gp)
      !call print("Predicted Error Factor :"//sum(abs(geometry_matrix(j,:)))//" Dimension :"//j)
   enddo

   call print("GP Uncertainty in Cartesian Coordinate: "//variance_xyz//" GP variance: "//error_gp//" variance Squre :"//error_gp_sq)
   pred_err(1, n_center_atom)=norm(variance_xyz)
   pred_err(2, n_center_atom)=abs(force_magntd_pred - norm(force))
   pred_err(3, n_center_atom)=error_ls
  if (update_data) then
          call system_timer('Updating the DATABSE : dist.dat force.dat ivs.dat')
          allocate(force_proj_test(k)) ! the projection of real force on test conf
          do j = 1, k
               force_proj_test(j) = dot_product(force_ptr(:,n_center_atom), ivs_set_norm_pred(j,:) )
          enddo
          call  write_iv_params(data_dir, r_grid, m_grid, species, k-add_vector, add_vector, n+1)
          if (write_data_dij) call write_distance(distance_confs_unsorted, data_dir)
          ! append distance file
          if (print_verbosity) then
                  write(*,*) 'Distance', distance_confs_unsorted(:)
          endif
          OPEN(2,status='old',file=trim(data_dir)//'Force.dat', form='UNFORMATTED', position='APPEND')
          OPEN(3,status='old',file=trim(data_dir)//'IV.dat', form='UNFORMATTED',  position='APPEND')
          write(2)  force_proj_test(:), norm(force_ptr(:,n_center_atom))
          do j=1, 3
                write(3)  ivs_set_pred(:,j,n_center_atom)
          enddo
	  close(2)
	  close(3)
          deallocate(force_proj_test)
          call system_timer('Updating the DATABSE : dist.dat force.dat ivs.dat')
   endif
   if (print_at) force_ptr(:,n_center_atom)=force
 enddo  ! loop over positions

  deallocate(ivs_set_pred)

   error_frame = maxval( pred_err)
   call print("Error Bar of the Frame :"//error_frame)

   if (print_at) then
       ! write the output configurations with predicted force, we dont want appending practically
       call write(at_in, out_file, append=.False.)
   endif

enddo                         ! loop over frames

call finalise(in)

deallocate(ivs_set_norm_pred)   ! because it is defined within the loop

deallocate(geometry_matrix)
deallocate(r_grid, m_grid, species, sigma)
deallocate(force_proj_ivs_pred)
deallocate(ivs_set_norm)
deallocate(force_magntd_data, force_magntd_data_select)
deallocate(ivs_set)
deallocate(mark_zero_ivs)
deallocate(force_proj_ivs)
deallocate(force_proj_ivs_select)
deallocate(covariance, inv_covariance)
deallocate(covariance_pred)
deallocate(distance_index, distance_confs_unsorted)
if(allocated(force_covariance_matrix)) deallocate(force_covariance_matrix)
if(allocated(covariance_tiny)) deallocate(covariance_tiny)
if(allocated(distance_matrix)) deallocate(distance_matrix)
if (allocated(distance_matrix_dij)) deallocate(distance_matrix_dij)
call system_finalise()
call print('-------------END OF FGP----------------')

end program force_gaussian_prediction

max_stars_repo_name: albapa/QUIP

max_issues_repo_name: albapa/QUIP

max_forks_repo_name: albapa/QUIP


text:
module vegetables_input_test_case_m
    use iso_varying_string, only: varying_string, operator(//), put_line, var_str
    use strff, only: operator(.includes.), to_string
    use vegetables_command_line_m, only: DEBUG
    use vegetables_input_m, only: input_t
    use vegetables_test_m, only: &
            filter_result_t, test_t, filter_failed, filter_matched
    use vegetables_test_case_result_m, only: test_case_result_t
    use vegetables_test_interfaces_m, only: computation_i, input_test_i
    use vegetables_result_m, only: fail
    use vegetables_test_result_item_m, only: test_result_item_t

    implicit none
    private
    public :: input_test_case_t

    type, extends(test_t) :: input_test_case_t
        private
        type(varying_string) :: description_
        procedure(input_test_i), nopass, pointer :: test
        logical :: has_setup_and_teardown
        procedure(computation_i), nopass, pointer :: setup
        procedure(computation_i), nopass, pointer :: teardown
    contains
        private
        procedure, public :: description
        procedure, public :: filter
        procedure, public :: num_cases
        procedure, public :: run_with_input
        procedure, public :: run_without_input
    end type

    interface input_test_case_t
        module procedure constructor_basic
        module procedure constructor_bracketed
    end interface
contains
    function constructor_basic(description, test) result(input_test_case)
        type(varying_string), intent(in) :: description
        procedure(input_test_i) :: test
        type(input_test_case_t) :: input_test_case

        input_test_case%description_ = description
        input_test_case%test => test
        input_test_case%has_setup_and_teardown = .false.
    end function

    function constructor_bracketed( &
            description, test, setup, teardown) result(input_test_case)
        type(varying_string), intent(in) :: description
        procedure(input_test_i) :: test
        procedure(computation_i) :: setup
        procedure(computation_i) :: teardown
        type(input_test_case_t) :: input_test_case

        input_test_case%description_ = description
        input_test_case%test => test
        input_test_case%has_setup_and_teardown = .true.
        input_test_case%setup => setup
        input_test_case%teardown => teardown
    end function

    pure function description(self)
        class(input_test_case_t), intent(in) :: self
        type(varying_string) :: description

        description = self%description_
    end function

    function filter(self, filter_string) result(filter_result)
        class(input_test_case_t), intent(in) :: self
        type(varying_string), intent(in) :: filter_string
        type(filter_result_t) :: filter_result

        if (self%description_.includes.filter_string) then
            filter_result = filter_matched(self)
        else
            filter_result = filter_failed()
        end if
    end function

    pure function num_cases(self)
        class(input_test_case_t), intent(in) :: self
        integer :: num_cases

        associate(unused => self)
        end associate

        num_cases = 1
    end function

    recursive function run_with_input(self, input) result(result_)
        class(input_test_case_t), intent(in) :: self
        class(input_t), intent(in) :: input
        type(test_result_item_t) :: result_

        if (DEBUG) call put_line( &
                "Beginning execution of: " // self%description_ &
                // merge(" on image " // to_string(this_image()), var_str(""), num_images() > 1))
        if (self%has_setup_and_teardown) call self%setup
        result_ = test_result_item_t(test_case_result_t( &
                self%description_, self%test(input)))
        if (self%has_setup_and_teardown) call self%teardown
        if (DEBUG) call put_line( &
                "Completed execution of: " // self%description_&
                // merge(" on image " // to_string(this_image()), var_str(""), num_images() > 1))
    end function

    function run_without_input(self) result(result_)
        class(input_test_case_t), intent(in) :: self
        type(test_result_item_t) :: result_

        if (DEBUG) call put_line( &
                "Improper use of: " // self%description_&
                // merge(" on image " // to_string(this_image()), var_str(""), num_images() > 1))
        result_ = test_result_item_t(test_case_result_t( &
                self%description_, fail("No input provided")))
    end function
end module

max_stars_repo_name: everythingfunctional/vegetables

max_issues_repo_name: everythingfunctional/vegetables

max_forks_repo_name: everythingfunctional/vegetables


subroutine read_goesglm(nread,ndata,nodata,infile,obstype,lunout,twindin,sis)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:  read_goesglm                reads lightning obs from a BUFR file
!   prgmmr:      k apodaca <karina.apodaca@colostate.edu>
!      org:      CSU/CIRA, Data Assimilation Group
!     date:      2015-03-12
!
! abstract:  This routine reads lightning data (Earth-relative location, frequency) from
!            file and prepares it for assimiliation in the form of lightning flash rate
!            (#hits km-2 hr-1).
!
!            Note: when running GSI in regional mode, the code only
!            retains those observations that fall within the regional
!            domain
!
! program history log:
!   2015-06-12  zupanski  - include the convert_to_flash_rate subroutine to convert
!                           lightning strike observations into lightning flash rate
!                           (#hits km-2 hr-1).
!   2015-11-18  apodaca   - include the convert_time subroutine to deal with computer
!                           precision dependencies.
!   2018-02-07  apodaca   - add further documentation
!
!   input argument list:
!     infile   - unit from which to read BUFR data
!     obstype  - observation type to process
!     lunout   - unit to which to write data for further processing
!
!   output argument list:
!     nread    - number of type "obstype" observations read
!     nodata   - number of individual "obstype" observations read
!     ndata    - number of type "obstype" observations retained for further processing
!     twindin  - input group time window (hours)
!     sis      - satellite/instrument/sensor indicator
!
! attributes:
!   language: Fortran 90 and/or above
!   machine:
!
!$$$
  use kinds, only: r_single,r_kind,r_double,i_kind
  use constants, only: zero,one_tenth,one,deg2rad,&
      three,rad2deg,&
      r60inv,ten
  use gridmod, only: diagnostic_reg,wrf_mass_regional,regional,nlon,nlat,&
      tll2xy,txy2ll,&
      rlats,rlons
  use lightinfo, only: iuse_light,nlighttype
  use obsmod, only: iadate
  use obsmod, only: offtime_data
  use gsi_4dvar, only: l4dvar,l4densvar,time_4dvar,winlen

  implicit none

! Declare passed variables
  character(len=*)                      ,intent(in   ) :: infile,obstype
  character(len=*)                      ,intent(in   ) :: sis
  integer(i_kind)                       ,intent(in   ) :: lunout
  integer(i_kind)                       ,intent(inout) :: nread,ndata
  real(r_kind)                          ,intent(in   ) :: twindin
  integer(i_kind)                                      :: nodata

! Declare local parameters
  real(r_kind),parameter:: r6   = 6.0_r_kind
  real(r_kind),parameter:: r90  = 90.0_r_kind
  real(r_kind),parameter:: r180 = 180.0_r_kind
  real(r_kind),parameter:: r360 = 360.0_r_kind

!--- Declare local variables
  logical lob
  logical outside

  character(40) hdstr,oestr,qcstr
  character(10) date
  character(8)  subset
  character(1)  sidchr(8)

  integer(i_kind) ireadmg,ireadsb,icntpnt,icount
  integer(i_kind) lunin,i
  integer(i_kind) itx
  integer(i_kind) ihh,idd,idate,iret,im,iy,k
  integer(i_kind) nchanl,nreal,ilat,ilon
  integer(i_kind) lqm
  integer(i_kind) iout
  integer(i_kind) ntest,nvtest
  integer(i_kind) minobs,minan
  integer(i_kind) ntb
  integer(i_kind) nmsg                ! message index
  integer(i_kind),parameter :: maxobs=2000000
  integer(i_kind),dimension(5):: idate5
  integer(i_kind),allocatable,dimension(:):: isort,iloc

  real(r_kind) time
  real(r_kind) usage
  real(r_kind) loe,lmerr
  real(r_kind) time_correction
  real(r_kind) dlat,dlon,dlat_earth,dlon_earth
  real(r_kind) cdist,disterr,disterrmax,rlon00,rlat00
  real(r_kind) vdisterrmax
  real(r_kind) timex,timeobs,toff,t4dv,zeps
  real(r_kind),allocatable,dimension(:,:):: cdata_all,cdata_out
!--- flash rate
  real(r_kind),allocatable,dimension(:,:):: cdata_flash,cdata_flash_h
  integer(i_kind)                        :: ndata_flash,ndata_flash_h

  real(r_double) rstation_id

  real(r_double),dimension(3):: hdr
  real(r_double),dimension(1,1):: qcmark,obserr

!  equivalence to handle character names
  equivalence(rstation_id,sidchr)

!--- data statements
  data hdstr  /'XOB YOB DHR'/
  data oestr  /'LOE'/
  data qcstr  /'LQM'/

  data lunin / 13 /

  nreal=13
  lob = obstype == 'goes_glm'

!                .      .    .                                       .

! Open, then read date from BUFR file

  call closbf(lunin)
  open(lunin,file=infile,form='unformatted')
  call openbf(lunin,'IN',lunin)
  call datelen(10)

! Initialization

  ntb = 0
  nmsg = 0
  disterrmax=-9999.0_r_kind

  allocate(cdata_all(nreal,maxobs),isort(maxobs))
  isort = 0
  cdata_all=zero
  nread=0
  ntest=0
  nvtest=0
  nchanl=0
  ilon=2
  ilat=3
  icntpnt=0

!                .      .    .                                       .
! Big loop over glmbufr file : READING THE BUFR FILE

  loop_msg: do while (ireadmg(lunin,subset,idate)== 0)
     nmsg = nmsg+1
     write(*,'(3a,i10)') 'subset=',subset,' cycle time =',idate

     loop_readsb: do while(ireadsb(lunin) == 0)
!       use msg lookup table to decide which messages to skip
!       use report id lookup table to only process matching reports
        ntb = ntb+1

!       Extract location and date information
        call ufbint(lunin,hdr,3,1,iret,hdstr)

        if(abs(hdr(2))>r90 .or. abs(hdr(1))>r360) cycle loop_readsb
        if(hdr(1) > r180)hdr(1)=hdr(1)-r360
        if(hdr(1) < zero)hdr(1)=hdr(1)+r360
        dlon_earth=hdr(1)*deg2rad
        dlat_earth=hdr(2)*deg2rad

        if (regional) then

!-- WRF-ARW

           if (wrf_mass_regional) then

              call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)    ! convert to rotated coordinates
              if(diagnostic_reg) then
                 call txy2ll(dlon,dlat,rlon00,rlat00)
                 ntest=ntest+1
                 cdist=sin(dlat_earth)*sin(rlat00)+cos(dlat_earth)*cos(rlat00)* &
                       (sin(dlon_earth)*sin(rlon00)+cos(dlon_earth)*cos(rlon00))
                 cdist=max(-one,min(cdist,one))
                 disterr=acos(cdist)*rad2deg
                 disterrmax=max(disterrmax,disterr)
              end if
              if(outside) cycle loop_readsb   ! check to see if outside regional domain

           endif ! wrf_mass_regional

        endif !if (regional) then

! Global

        if (.not. regional) then
            dlat = dlat_earth
            dlon = dlon_earth
            call grdcrd1(dlat,rlats,nlat,1)
            call grdcrd1(dlon,rlons,nlon,1)
        endif !  end global block

        if (offtime_data) then

!           in time correction for observations to account for analysis
!           time being different from obs file time.
            write(date,'( i10)') idate
            read (date,'(i4,3i2)') iy,im,idd,ihh
            idate5(1)=iy
            idate5(2)=im
            idate5(3)=idd
            idate5(4)=ihh
            idate5(5)=0
            call w3fs21(idate5,minobs)    !  obs ref time in seconds relative to historic date
            idate5(1)=iadate(1)
            idate5(2)=iadate(2)
            idate5(3)=iadate(3)
            idate5(4)=iadate(4)
            idate5(5)=0
            call w3fs21(idate5,minan)    !  analysis ref time in seconds relative to historic date

!           Add obs reference time, then subtract analysis time to get obs time relative to analysis
            time_correction=float(minobs-minan)*r60inv
        else
            time_correction=zero
        end if

        timeobs=real(real(hdr(3),r_single),r_double)
        t4dv=timeobs + toff
        zeps=1.0e-8_r_kind
        if (t4dv<zero  .and.t4dv>      -zeps) t4dv=zero
        if (t4dv>winlen.and.t4dv<winlen+zeps) t4dv=winlen
        t4dv=t4dv + time_correction
        time=timeobs + time_correction

        if (l4dvar.or.l4densvar) then
           if (t4dv<zero.or.t4dv>winlen) cycle loop_readsb ! outside time window
        else
           if((real(abs(time)) > real(twindin)))cycle loop_readsb ! outside time window
        endif

        timex=time

!       Note: An assessment of the GLM detection error is still a work in
!       progress and at present, there are no clear metrics to assign a
!       "quality mark" for an effective quality control procedure for
!       lightning observations from this sensor. Nonetheless, the infrastructure
!       for passing observations errors and quality control information to other
!       routines within the GSI source code has been left in place. In addition, a
!       temporarty undefined value (-.9999) has been assigned to the mnemonics
!       corresponding to these variables in the lightning BUFR file. Developments
!       for the calculation of observation errors and quality control (sanity checks) are
!       expected in future upgrades to the "GOES/GLM variational lightning
!       assimilation package."

!       Extract observation error informatiom
        call ufbint(lunin,obserr,1,1,iret,oestr)

        loe=obserr(1,1)
        lmerr=loe

!       Extract quality control information
        call ufbint(lunin,qcmark,1,1,iret,qcstr)

        lqm=qcmark(1,1)

!       Data counter

        nread=nread+1
        icntpnt=icntpnt+1

        ndata=ndata+1
        nodata=nodata+1
        iout=ndata
        isort(icntpnt)=iout

        if (ndata > maxobs) then
           write(6,*)'READ_GOESGLM:  ***WARNING*** ndata > maxobs for ',obstype
           ndata = maxobs
        end if

!       Set usage variable
        usage = zero

        if (iuse_light(nlighttype) <= 0)usage=100._r_kind
        if (lob) then
           cdata_all(1,iout) =loe                  ! lightning observation error
           cdata_all(2,iout) =dlon                 ! grid relative longitude
           cdata_all(3,iout) =dlat                 ! grid relative latitude
           cdata_all(4,iout) =iout                 ! lightning obs
           cdata_all(5,iout) =rstation_id          ! station id
           cdata_all(6,iout) =t4dv                 ! analysis time
           cdata_all(7,iout) =nlighttype           ! type
           cdata_all(8,iout) =lmerr                ! lightning max error
           cdata_all(9,iout) =lqm                  ! quality mark
           cdata_all(10,iout)=loe                  ! original lightning obs error loe
           cdata_all(11,iout)=usage                ! usage parameter
           cdata_all(12,iout)=dlon_earth*rad2deg   ! earth relative lon (degrees)
           cdata_all(13,iout)=dlat_earth*rad2deg   ! earth relative lat (degrees)
        end if

! end loop on read line BUFR

     end do loop_readsb

! end of BUFR read loop

  enddo loop_msg   !Uncomment if reading messeges in a loop

!                .      .    .                                       .

! Close unit to bufr file

  call closbf(lunin)

! Write header record and data to output file for further processing
  allocate(iloc(ndata))
  icount=0.
  do i=1,maxobs
     if(isort(i) > 0)then
       icount=icount+1
       iloc(icount)=isort(i)
     end if
  end do
  if(ndata /= icount)then
     write(6,*) ' READ_GOESGLM: mix up in read_goesglm ,ndata,icount ',ndata,icount
     call stop2(50)
  end if

  allocate(cdata_out(nreal,ndata))
  do i=1,ndata
     itx=iloc(i)
     do k=1,nreal
        cdata_out(k,i)=cdata_all(k,itx)
     end do
  end do

  deallocate(iloc,isort,cdata_all)

!                .      .    .                                       .

! Call to the subroutine that transforms lightning strikes into lightning flash rate

  if(ndata /= 0) then

      !! count flash rate data and alocate temporary domain
      !! begin with the current number of strikes as the theoretical upper limit

     ndata_flash_h=ndata

     allocate(cdata_flash_h(nreal,ndata_flash_h))

     call convert_to_flash_rate   &
              (nreal,ndata,cdata_out,ndata_flash_h,cdata_flash_h,ndata_flash)

     deallocate(cdata_out)
     ndata=ndata_flash
     allocate(cdata_flash(nreal,ndata))

     do i=1,ndata
        do k=1,nreal
           cdata_flash(k,i)=cdata_flash_h(k,i)
        end do
     end do

     deallocate(cdata_flash_h)

     write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
     write(lunout) cdata_flash

!!!!!!!!!
!!! Write lat, lon, time, lightning flash rate "superobs" into a file:
!!! cdata_flash_h(2,iout),cdata_flash_h(3,iout),cdata_flash_h(4,iout), &
!!! cdata_flash_h(6,iout)cdata_flash(4,:)

     deallocate(cdata_flash)

  else  ! ndata=0

     write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
     write(lunout) cdata_out
     deallocate(cdata_out)

  end if  !!  if(ndata =/ 0) then

900 continue
  if(diagnostic_reg .and. ntest>0) write(6,*)'READ_GOESGLM:  ',&
     'ntest,disterrmax=',ntest,disterrmax
  if(diagnostic_reg .and. nvtest>0) write(6,*)'READ_GOESGLM:  ',&
     'nvtest,vdisterrmax=',ntest,vdisterrmax

  if (ndata == 0) then
     call closbf(lunin)
     write(6,*)'READ_GOESGLM:  closbf(',lunin,')'
  endif

  close(lunin)

  close(55)

! End of routine
  return

end subroutine read_goesglm

!!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

subroutine convert_to_flash_rate   &
   (nreal,ndata_strike,cdata_strike,ndata_flash_h,cdata_flash_h,ndata_flash)

!$$$  documentation block
!                .      .    .                                       .
! subroutine:  convert_to_flash_rate      converts geo-located lightning strikes into
!                                         lightning flash rate (#hits km-2 hr-1)

!   prgmmr: zupanski <milija.zupanski@colostate.edu>
!      org: CSU/CIRA, Data Assimilation group
!     date: 2015-06-12
!
! abstract:  This subroutine does the following:
!----
!- 1- counts the number of hits surrounding a GSI analysis grid point
!- 2- calculates the flash rate averaged over time and area (hr*km**2)
!- 3- finds the center of mass in terms of lon,lat and glon,glat
!- 4- assigns a center of mass to be the flash rate observation point
!----
!
! program history log:
!   2015-07-01  apodaca  - several updates in the calculation of lightning flashrate

  use kinds, only: r_single,r_kind,r_double,i_kind
  use constants, only: zero,one_tenth,one,deg2rad,&
      three,half,zero,&
      r10,r100,ten, r1000, rearth
  use gridmod, only: wrf_mass_regional,regional,nlon,nlat,&
      tll2xy,txy2ll
  use gridmod, only: lat2, lon2
  use wrf_mass_guess_mod, only: ges_xlon, ges_xlat
  use gsi_4dvar, only: nhr_assimilation

  implicit none

  integer(i_kind),intent(in   )               :: nreal,ndata_strike,ndata_flash_h
  integer(i_kind),intent(inout)               :: ndata_flash
  real(r_kind),intent(inout),dimension(nreal,ndata_strike)  :: cdata_strike
  real(r_kind),intent(inout),dimension(nreal,ndata_flash_h) :: cdata_flash_h

  real(r_kind),allocatable,dimension(:) :: gtim_central
  real(r_kind),allocatable,dimension(:) :: glon_central
  real(r_kind),allocatable,dimension(:) :: glat_central
  real(r_kind),allocatable,dimension(:) ::  lon_central
  real(r_kind),allocatable,dimension(:) ::  lat_central
  integer(i_kind),allocatable,dimension(:) ::  ind_min
  integer(i_kind),allocatable,dimension(:) ::  lcount

  real(r_kind)      :: rearth2
  real(r_kind)      :: dtime,darea,cosine
  real(r_kind)      :: darea_sum
  real(r_kind)      :: delta_lon,delta_lat
  real(r_kind)      :: lat_ref
  real(r_kind)      :: xx,yy
  real(r_kind)      :: dist2,dist_min
  integer(i_kind)   :: ii0,jj0
  integer(i_kind)   :: ngridh
  integer(i_kind)   :: index
  integer(i_kind)   :: iobs,usage
  logical           :: xflag,yflag

  real(r_kind)      :: xbound,ybound
  integer(i_kind)   :: nxdim,nydim
  integer(i_kind)   :: icount

  !  Output files
!----

!- unit grid lat-lon distances (Earth difference divided by grid difference)
!- calculated from all averaged grid areas
!!- note: this can be relaxed if the unit grid area (km) is known (i.e. darea)
!! darea = (r*cos(lat)*dlon)*(dlat)

  rearth2=(rearth/r1000)**2  !! squared earth radius in km (need for darea calculation)

  if (ndata_strike>0) then

     darea_sum=0._r_kind
     do iobs=1,ndata_strike

        ii0=INT(cdata_strike(2,iobs))
        jj0=INT(cdata_strike(3,iobs))

        delta_lon=ges_xlon(jj0,ii0+1,1)-ges_xlon(jj0,ii0,1)
        delta_lat=ges_xlat(jj0+1,ii0,1)-ges_xlat(jj0,ii0,1)

        lat_ref  =half*(cdata_strike(13,iobs)+cdata_strike(13,iobs-1))

        cosine=cos(lat_ref*deg2rad)

        darea_sum=darea_sum+rearth2*cos(lat_ref*deg2rad)*(abs(delta_lon)*deg2rad)*&
                  (abs(delta_lat)*deg2rad)

     end do  !! do iobs=2,ndata_strike

     darea=darea_sum/float(ndata_strike)
  else   !! ndata_strike=0
     darea=zero

  end if  !! if(ndata_strike>0) then

  dtime=float(nhr_assimilation)

  ! Regional

  if (regional) then

!-- WRF-ARW

     if (wrf_mass_regional) then

        nxdim=lon2
        nydim=lat2

     endif ! wrf_mass_regional

  endif !if (regional) then

! Global

  if (.not. regional) then

     nxdim=nlon
     nydim=nlat

  endif !  end global block

!!! Allocate new var for flash rate here (nxdim,nydim)
! asign zero to all points
! update the relevant points in the loop

  ngridh=nxdim*nydim

  allocate(gtim_central(1:ngridh))
  allocate(glon_central(1:ngridh))
  allocate(glat_central(1:ngridh))
  allocate( lon_central(1:ngridh))
  allocate( lat_central(1:ngridh))
  allocate(      lcount(1:ngridh))

  lcount(:)=0
  glon_central(:)=zero
  glat_central(:)=zero
  lon_central(:) =zero
  lat_central(:) =zero
  gtim_central(:)=zero

  do iobs=1,ndata_strike

     xx=cdata_strike(2,iobs)   !! glon
     yy=cdata_strike(3,iobs)   !! glat
     ii0=INT(cdata_strike(2,iobs))
     jj0=INT(cdata_strike(3,iobs))

!!   find lightning strikes near the (ii0,jj0) point

     xbound=float(ii0)
     ybound=float(jj0)

     xflag=(xx>xbound) .AND. (xx<xbound+1.)
     yflag=(yy>ybound) .AND. (yy<ybound+1.)

     if (xflag .AND. yflag ) then

        index=(jj0-1)*nxdim+ii0
        lcount(index)=lcount(index)+1

        glon_central(index)=glon_central(index)+cdata_strike(2,iobs)
        glat_central(index)=glat_central(index)+cdata_strike(3,iobs)
        lon_central(index)= lon_central(index)+cdata_strike(12,iobs)
        lat_central(index)= lat_central(index)+cdata_strike(13,iobs)

        if (lcount(index)<2) then
            gtim_central(index)=cdata_strike(6,iobs)
        else
            call convert_time (gtim_central(index),cdata_strike(6,iobs),lcount(index))
        end if

     end if  !! if(xflag .AND. yflag )

  enddo !! do iobs=1,ndata_strike

!-- find the center of mass

  do index=1,ngridh
     if (lcount(index)>0) then
        glon_central(index)=glon_central(index)/float(lcount(index))
        glat_central(index)=glat_central(index)/float(lcount(index))
        lon_central(index)= lon_central(index)/float(lcount(index))
        lat_central(index)= lat_central(index)/float(lcount(index))
     endif  !! if(lcount(index)>0) then
  enddo  !! do index=1,ngridh

!-- find the original index of the nearest strike (need for transfer of input obs)

  allocate(ind_min(1:ngridh))
  ind_min(:)=-99

  dist_min=1.e10_r_kind
  do iobs=1,ndata_strike

     xx=cdata_strike(2,iobs)   !! glon
     yy=cdata_strike(3,iobs)   !! glat
     ii0=INT(cdata_strike(2,iobs))
     jj0=INT(cdata_strike(3,iobs))
     index=(jj0-1)*nxdim+ii0

     if (lcount(index)>0) then
        dist2=(xx-glon_central(index))**2+(yy-glat_central(index))**2
        if (dist2<dist_min) then
           dist_min=dist2
           ind_min(index)=iobs
        end if
     end if  !! if(lcount(index)>0) then

  enddo  !! do iobs=1,ndata_strike

!----
!---- Output
!----

!-- count the non-zero flash rates and assign a temporary domain cdata_flash_h
!-- Note: it is assumed that only non-zero flash rates are true observations

  icount=0
  do index=1,ngridh
     if (lcount(index)>0) then
        icount=icount+1
        cdata_flash_h( 1,icount)=cdata_strike( 1,ind_min(index))
        cdata_flash_h( 2,icount)=glon_central(index)
        cdata_flash_h( 3,icount)=glat_central(index)

        if (darea>0._r_kind) then
           cdata_flash_h( 4,icount)=float(lcount(index))/(darea*dtime)
        else
           cdata_flash_h( 4,icount)=0.
        end if

        cdata_flash_h( 5,icount)=cdata_strike( 5,ind_min(index))
        cdata_flash_h( 6,icount)=gtim_central(index)
        cdata_flash_h( 7,icount)=cdata_strike( 7,ind_min(index))
        cdata_flash_h( 8,icount)=cdata_strike( 8,ind_min(index))
        cdata_flash_h( 9,icount)=cdata_strike( 9,ind_min(index))
        cdata_flash_h(10,icount)=cdata_strike(10,ind_min(index))
        cdata_flash_h(11,icount)=usage
        cdata_flash_h(12,icount)=lon_central(index)
        cdata_flash_h(13,icount)=lat_central(index)

     endif  !! if(lcount(index)>0) then

  enddo  !! do index=1,ngridh

  ndata_flash=icount

  deallocate(ind_min)
  deallocate(lcount)
  deallocate(gtim_central)
  deallocate(glon_central)
  deallocate(glat_central)
  deallocate(lon_central)
  deallocate(lat_central)

!-----
! End of routine
  return

end subroutine convert_to_flash_rate

!-----
subroutine convert_time (date_old,date_new,nmax)

!$$$  documentation block
!                .      .    .                                       .
! subroutine:  convert_time
! prgmmr: k apodaca               date: 2015-11-18
!
! abstract:  This subroutine performs a date conversion to deal with
!            computer precission dependencies associated with a 10-digit
!            float for the analysis date/time.
!--
  use kinds, only: r_kind,i_kind

  implicit none

  integer(i_kind), intent(in) :: nmax
  real(r_kind), intent(inout) :: date_old
  real(r_kind), intent(in)  :: date_new
  integer(i_kind) :: i,sumidd
  integer(i_kind) :: idd,jdd,kdd
  real(r_kind), allocatable :: xdate(:)
  real(r_kind) :: dd,hh,ysumidd,xsumidd
  real(r_kind) :: xdd,xhh,ydate
  real(r_kind) :: xccyy

  allocate(xdate(1:nmax))

  xdate(1:nmax-1) = date_old
  xdate(nmax) = date_new

  sumidd=0._r_kind
  do i=1,nmax
     xccyy = INT(1.0e-8_r_kind*xdate(i))*1.0e8_r_kind
     xdate(i) = INT(xdate(i))-xccyy

     jdd=INT(0.0001_r_kind*xdate(i))
     idd=INT(xdate(i))-jdd*10000

     ysumidd=float(idd)
     dd=float(INT(0.01_r_kind*ysumidd))
     hh=ysumidd-dd*100._r_kind

     sumidd=sumidd+dd*24._r_kind+hh

  enddo  !! do i=1,nmax

  xsumidd=float(sumidd)/nmax
  ysumidd=float(INT(xsumidd))

  kdd=INT(xsumidd/24._r_kind)
  xdd=float(kdd)
  xhh=ysumidd-float(kdd)*24._r_kind

  ydate=float(jdd)*10000._r_kind+xdd*100._r_kind+xhh+xccyy

  date_old=ydate

  deallocate(xdate)

end subroutine convert_time


text:
      SUBROUTINE T2DED(SCR,LSCR)

!                             PROCESS T2D  (TIME TO DEPTH CONVERSION)
!                             ------- ---
!
!  DOCUMENT DATE: 19 November 1990
!
!     T2D TRANSFORMS TIME DOMAIN TRACES TO DEPTH DOMAIN TRACES USING USER
!  DEFINED VELOCITY FUNCTIONS.  THE VELOCITIES MAY BE VARIED TEMPORALLY AND
!  SPATIALLY.  THE VELOCITIES MAY BE AVERAGE VELOCITIES OR INTERVAL VELOCITIES.
!  INTERVAL VELOCITIES ARE CONVERTED TO AVERAGE VELOCITIES.  DEPTH
!  CONVERSION UTILIZES THE FORMULA     D=T/2*V
!  INTERVAL VELOCITY DEPTH CONVERSION UTILIZES THE FORMULA
!  D(N+1)=(T(N+1)-T(N))/2*V(N+1).
!     THE VELOCITY FUNCTION MAY BE VARIED SPATIALLY BETWEEN VELOCITY CONTROL
!  POINTS.  VELOCITY CONTROL POINTS ARE DEFINED AS EITHER SHOT POINT NUMBERS
!  OR RP NUMBERS. EACH PARAMETER LIST HAS A START AND END CONTROL POINT WITH THE
!  PARAMETERS BEING CONSTANT FOR ALL POINTS BETWEEN THE FIRST AND LAST CONTROL
!  POINT OF THE LIST.  E.G.  IF THE FIRST CONTROL POINT=100 AND THE LAST CONTROL
!  POINT=110, THE POINTS 101 TO 109 WILL ALSO HAVE THE SAME PARAMETERS APPLIED.
!                 E.G. FNO 1 VTP 2100 1.1 3100 2.1 END
!                      FNO 3 VTP 2300 1.3 3300 2.3 END
!                      FNO 4 VTP 2200 1.0 2500 1.5 3000 2.0B END
!      RESULTS IN FNO 2 VTP 2200 1.2 3200 2.2
!      THE AVERAGE VELOCITY FUNCTION IS EVALUATED FOR EVERY TIME SAMPLE FOR THE
!  CONTROL POINTS.  TRACES BETWEEN VELOCITY CONTROL POINTS RECEIVE A VELOCITY
!  FUNCTION THAT IS THE LINEAR INTERPOLATION OF THE VELOCITIES OF THE
!  CORRESPONDING TIMES ON THE ADJACENT CONTROL POINTS.
!     THE OUTPUT SAMPLE INTERVAL IS CHANGED BY PROCESS T2D TO 1 MIL (.001 SECOND).
!  THUS, IN ORDER TO PLOT DEPTH DATA, USE PLOT PARAMETERS FOR 1 MIL TIME DATA.
!     SPATIAL VARIATION IS BY SHOT IF THE DATA IS SORTED BY SHOT AND IS
!  VARIED BY RP IF THE DATA IS SORTED BY RP.
!     The output units of T2D are meters.  Whenever the SEGY header
!  uses milliseconds, T2D converts it to meters.  Whenever the SEGY header
!  uses seconds, T2D converts it to kilometers.
!     EACH PARAMETER LIST MUST BE TERMINATED WITH THE WORD END.  THE ENTIRE SET
!  OF T2D PARAMETERS MUST BE TERMINATED BY THE WORD END.
!
!    EXAMPLE:
!       T2D
!            SDEPTH 0 EDEPTH 20000 OSI 10 VTYPE AVE
!           FNO 1 VTP 1500 2.2 2200 4.5 3500 6.7 4500 7.0 END
!           FNO 223 VTP 1500 2.9 2100 4.2 3500 6.6 4500 7.0 END
!           FNO 253 VTP 1500 3.0 2100 4.9 3500 6.5 4500 7.0 END
!       END
!
!  THE PARAMETER DICTIONARY
!  --- --------- ----------
!  FNO    - THE FIRST SHOT (OR RP) TO APPLY THE VELOCITIES TO.  SHOT (RP) NUMBERS
!           MUST INCREASE MONOTONICALLY.
!           PRESET=1
!  LNO    - THE LAST SHOT (RP) NUMBER TO APPLY THE VELOCITIES TO.  LNO MUST BE
!           LARGER THAN FNO IN EACH LIST AND MUST INCREASE LIST TO LIST.
!           DEFAULT=FNO
!  VTYPE  - THE TYPE OF VELOCITY FUNCTION GIVEN IN THE VTP PARAMETER.
!         = AVE,  AVERAGE VELOCITIES.
!         = INT,  INTERVAL VELOCITIES.
!           PRESET = INT.  E.G. VTYPE AVE
!  VTP    - VELOCITY-TIME-PAIRS.  A LIST OF VELOCITY AND TWO-WAY TRAVE TIME
!           (IN SECONDS) PAIRS.  VTP MUST BE GIVEN IN EACH T2D PARAMETER
!           LIST.  A MAXIMUM OF 25 PAIRS MAY BE GIVEN.  DATA TIMES BEFORE
!           THE FIRST GIVEN TIME IN VTP WILL RECEIVE THE FIRST GIVEN VELOCITY.
!           LIKEWISE, DATA TIMES EXCEEDING THE LAST GIVEN TIME IN VTP WILL
!           RECEIVE THE LAST VELOCITY GIVEN IN VTP.
!           DEFAULT=NONE. E.G. VTP 1490 1.0 2000 2.0
!  SDEPTH - THE START DEPTH OF THE OUTPUT DEPTH SECTION.  When this is,
!           the SEGY header is changed so that the deep water delay = SDEPTH
!           PRESET = 0.  E.G. SDEPTH 1000.
!  EDEPTH - THE END DEPTH OF THE OUTPUT DEPTH SECTION.  ALL OUTPUT TRACES WILL
!           HAVE THE SAME END DEPTH.
!           PRESET = 10000.  E.G. EDEPTH 8000
!  OSI    - THE OUTPUT SAMPLE INTERVAL EXPRESSED AS THE DISTANCE BETWEEN SAMPLES.
!           The output sample interval should be similar to time sample
!           intervals.  e.g.  osi 4 gets written in the SEGY as 4 meters
!           whereas in time it would be 4 mils.
!           PRESET = 5. (METERS).  E.G. OSI 10.
!  ADDWB  - WHEN GIVEN A VALUE OF YES, THE WATER BOTTOM TIME IS ADDED TO THE
!           VELOCITY FUNCTION AFTER SPATIAL VARIATION HAS BEEN DONE.
!           PRESET=NO
!  END    - TERMINATES EACH PARAMETER LIST.
!
!
!  WRITTEN AND COPYRIGHTED BY:
!  PAUL HENKART, SCRIPPS INSTITUTION OF OCEANOGRAPHY, JULY 1983
!  ALL RIGHTS ARE RESERVED BY THE AUTHOR.  PERMISSION TO COPY OR REPRODUCE THIS
!  SUBROUTINE, BY COMPUTER OR OTHER MEANS, MAY BE OBTAINED ONLY FROM THE AUTHOR.
!  mod 12 oct 89 to convert interval vtps to average vtps.
!  mod 10 mar 90 by pch to to the above on interval velocities ONLY!
!  mod 6 Aug 97 - Add check for OSI > 32 because OSI*1000 must fit in a
!       SEG-Y 16 bit header word.
!  mod 6 Oct 06 - Change sdepth & edepth preset to -1, so each they can change
!                 on every trace like delay and nsamps do
!  mod 6 Oct 06 - Change OSI preset to 1 (like as in 1m, or 1mil)
!
!  ARGUMENTS:
!  SCR    - A SCRATCH ARRAY AT LEAST 60 32BIT WORDS LONG.
!  LSCR   - THE SAME ARRAY BUT USED FOR LONG INTEGERS.  PRIME DOESN'T ALLOW
!           EQUIVALENCING OF ARGUMENTS.
!
!  DISC PARAMETER LIST ORDER:
!  WORD 1)  FNO      (ALL ENTRIES ARE 32 BITS LONG)
!       2   LNO
!       3)  OSI
!       4)  SDEPTH
!       5)  EDEPTH
!       6)  ADDWB
!       7)  LPRINT
!       8)  NVTPS  - THE NUMBER OF ENTRIES IN THE VTP GIVEN BY THE USER
!       9)  ODELAY - THE DELAY OF THE OUTPUT
!       10)  VTYPE
!       11-60) VTP  - ALWAYS 50 LONG
!
!
!
      PARAMETER (NPARS=10)                                              !* THE NUMBER OF USER PARAMETERS
      PARAMETER (MAXVTP=50)                                             ! THE MAXIMUM NUMBER OF VTPS THAT T2DEX CAN HANDLE
      PARAMETER (NWRDS=MAXVTP+NPARS)                                    ! THE NUMBER OF WORDS IN EACH DISK OPARAMETER LIST
      PARAMETER (MULTIV=10)                                             ! THE INDEX OF THE FIRST MULTIVALUED PARAMETER
      CHARACTER*6 NAMES(NPARS)
      CHARACTER*1 TYPES(NPARS)
      DIMENSION LENGTH(NPARS)
      CHARACTER*80 TOKEN
      DIMENSION VALS(NPARS),LVALS(NPARS)
      COMMON /EDITS/ IERROR,IWARN,IRUN,NOW,ICOMPT
      COMMON /T2DCOM/ MUNIT,NLISTS
      COMMON /READT/ILUN,NUMHDR,NUMDAT,IUNHDR,IREELM,INTRCS,IFMT,NSKIP,
     *   SECS,LRENUM,ISRCF,IDTYPE
      DIMENSION SCR(111),LSCR(111)
      INTEGER FNO,SDEPTH,EDEPTH,VTYPE
!
!
      EQUIVALENCE (FNO,LVALS(1)),
     2            (LNO,LVALS(2)),
     3            (ADDWB,LVALS(3)),
     4            (LPRINT,LVALS(4)),
     5            (OSI,VALS(5)),
     6            (SDEPTH,LVALS(6)),
     7            (EDEPTH,LVALS(7)),
     9            (VTYPE,LVALS(9)),
     *            (VTP,VALS(10))
      DATA NAMES/'FNO   ',
     2           'LNO   ',
     3           'ADDWB ',
     4           'LPRINT',
     5           'OSI   ',
     6           'SDEPTH',
     7           'EDEPTH',
     8           'ODELAY',
     9           'VTYPE ',
     *           'VTP   '/
      DATA LENGTH/3,3,5,6,3,6,6,6,5,3/
      DATA TYPES/'L','L','A','L','F','L','L','F','A','F'/
      DATA NVTPS/0/, NVTPS1/-1/
!****
!****      SET THE PRESETS
!****
      FNO=1
      LNO=9999999
      VTP=-1.
      OSI=1.                                                            ! 5. METERS BETWEEN SAMPLES
      SDEPTH=-1.
      EDEPTH=-1
      VTYPE=2                                                           ! 1= INTERVAL, 2= AVERAGE
      LPRINT=0
      IADDWB=0
      LLNO = 0
      NLISTS=0
      NS=0
!****
!****   GET A PARAMETER FILE
!****
      CALL GETFIL(1,MUNIT,TOKEN,ISTAT)
!****
!****   THE CURRENT COMMAND LINE IN THE SYSTEM BUFFER MAY HAVE THE PARAMETERS.
!****   GET A PARAMETER LIST FROM THE USER.
!****
      NTOKES=1
  100 CONTINUE
      CALL GETOKE(TOKEN,NCHARS)                                         ! GET A TOKEN FROM THE USER PARAMETER LINE
      CALL UPCASE(TOKEN,NCHARS)                                         ! CONVERT THE TOKEN TO UPPERCASE
      IF(NCHARS.GT.0) GO TO 150
      IF(NOW.EQ.1) PRINT 140
  140 FORMAT(' <  ENTER PARAMETERS  >')
      CALL RDLINE                                                       ! GET ANOTHER USER PARAMETER LINE
      NTOKES=0
      GO TO 100
  150 CONTINUE
      NTOKES=NTOKES+1
      DO 190 I=1,NPARS                                                  ! SEE IF IT IS A PARAMETER NAME
      LEN=LENGTH(I)                                                     ! GET THE LEGAL PARAMETER NAME LENGTH
      IPARAM=I                                                          ! SAVE THE INDEX
      IF(TOKEN(1:NCHARS).EQ.NAMES(I)(1:LEN).AND.NCHARS.EQ.LEN) GO TO 200
  190 CONTINUE                                                          ! STILL LOOKING FOR THE NAME
      IF(TOKEN(1:NCHARS).EQ.'END'.AND.NCHARS.EQ.3) GO TO 1000           ! END OF PARAM LIST?
      IF(NS.NE.0) GO TO 230
      PRINT 191, TOKEN(1:NCHARS)
  191 FORMAT(' ***  ERROR  *** T2D DOES NOT HAVE A PARAMETER ',
     *  'NAMED ',A10)
      IERROR=IERROR+1
      GO TO 100
!****
!****    FOUND THE PARAMETER NAME, NOW FIND THE VALUE
!****
  200 CONTINUE
      NS=0
      NPARAM=IPARAM
  210 CONTINUE                                                          !  NOW FIND THE VALUE
      CALL GETOKE(TOKEN,NCHARS)
      CALL UPCASE(TOKEN,NCHARS)
      NTOKES=NTOKES+1
      IF(NCHARS.GT.0) GO TO 230                                         ! END OF LINE?
      IF(NOW.EQ.1) PRINT 140                                            ! THIS ALLOWS A PARAMETER TO BE ON A DIFFERENT LINE FROM THE NAME
      CALL RDLINE                                                       ! GET ANOTHER LINE
      NTOKES=0
      GO TO 210
  230 CONTINUE
      IF(TYPES(NPARAM).NE.'A') GO TO 240
      IF(NAMES(NPARAM).EQ.'ADDWB'.AND.TOKEN(1:NCHARS).EQ.'YES')
     *    IADDWB=1
      IF(NAMES(NPARAM).EQ.'VTYPE'.AND.TOKEN(1:NCHARS).EQ.'AVE') VTYPE=2
      GOTO 100
  240 CONTINUE
      CALL DCODE(TOKEN,NCHARS,AREAL,ISTAT)                              ! TRY AND DECODE IT
      IF(ISTAT.EQ.2) GO TO 420                                          ! =2 MEANS IT IS A NUMERIC
      IERROR=IERROR+1                                                   ! DCODE PRINTED AN ERROR
      GO TO 100
  420 IF(TYPES(NPARAM).EQ.'L') GO TO 500
      IF(NPARAM.LT.MULTIV) GO TO 490                                    !  IS IT A MULTIVALUED PARAMETER
      NS=NS+1                                                           !  THE TOKEN WAS A MULTI-VALUED PARAMETER
      nvtps = ns
      ITEMP=MULTIV
      SCR(NS+ITEMP)=AREAL
      GO TO 100
  490 VALS(NPARAM)=AREAL                                                !  FLOATING POINT VALUES
      GO TO 100
  500 CONTINUE                                                          !  32 BIT INTEGER VALUES
      LVALS(NPARAM)=AREAL
      GO TO 100
!****
!****   FINISHED A LIST, NOW DO THE ERROR AND VALIDITY CHECKS
!****
 1000 CONTINUE                                                          ! MAKE SURE ALL SHOT & RP NUMBERS INCREASE
      IF(LNO.EQ.9999999) LNO=FNO                                        ! DEFAULT LNO TO FNO
      IF(FNO.GT.LLNO) GO TO 1020                                        !  IS FNO LARGER THAN THE LAST LNO
      PRINT 1010
 1010 FORMAT(' ***  ERROR  ***  SHOT AND RP NUMBERS MUST INCREASE.')
      IERROR=IERROR+1
 1020 IF(LNO.GE.FNO) GO TO 1030                                         ! DO THEY INCREASE IN THIS LIST
      PRINT 1010
      IERROR=IERROR+1
 1030 IF( nvtps .EQ. 0 ) THEN
         PRINT *,' ***  ERROR  ***  VTP NOT GIVEN.'
         IERROR=IERROR+1
      ENDIF
! 1070 IF( SDEPTH .GE.-1.) GO TO 1080
!      PRINT 1075
! 1075 FORMAT(' ***  ERROR  ***  SDEPTH MUST BE POSITIVE.')
!      IERROR=IERROR+1
! 1080 IF(EDEPTH.GT.0.) GO TO 1090
!      PRINT 1085
! 1085 FORMAT(' ***  ERROR  ***  EDEPTH IS REQUIRED.')
!      IERROR=IERROR+1
! 1090 IF((EDEPTH-SDEPTH)/OSI.GT.100.) GO TO 1100
!      PRINT 1095
! 1095 FORMAT(' ***  WARNING  ***  UNUSUALLY SHORT OUTPUT REQUESTED.')
!      IWARN=IWARN+1
 1100 CONTINUE
      IF(OSI.GT.0.) GO TO 1110
      PRINT 1105,OSI
 1105 FORMAT(' ***  ERROR  ***  INVALID VALUE OF ',F10.5,' FOR OSI.')
      IERROR=IERROR+1
 1110 CONTINUE
      IF( osi .GT. 32 ) THEN
          PRINT *,' ***  ERROR  ***  OSI must be less than 32.'
          ierror = ierror + 1
      ENDIF
      DO 1200 II=1,NVTPS,2                                              !  CHECK THE VTP FOR ERRORS
      I=II+NPARS
      SCR(I)=SCR(I)/2.                                                  ! ****  CONVERT TWO WAY TRAVE TIMES TO ONE WAY!!!!!
      IF(SCR(I).GT.0) GO TO 1160
      PRINT 1150,SCR(I)
 1150 FORMAT(' ***  ERROR  ***  ILLEGAL VTP VELOCITY OF ',F10.4)
      IERROR=IERROR+1
 1160 J=I+1
      IF(SCR(J).GE.0.AND.SCR(J).LT.20.) GO TO 1180
      PRINT 1170,SCR(J)
 1170 FORMAT(' ***  ERROR  ***  ILLEGAL VTP TIME OF ',F10.4)
      IERROR=IERROR+1
 1180 IF(II.EQ.1) GO TO 1200
      IF(SCR(J).GT.SCR(J-2)) GO TO 1200
      PRINT 1190,SCR(J)
 1190 FORMAT(' ***  ERROR  ***  THE VTP TIME OF ',F10.4,' DECREASED.')
      IERROR=IERROR+1
 1200 CONTINUE
      NVTPS1=NVTPS
      LLNO=LNO
!****
!****      WRITE THE PARAMETER LIST TO DISC
!****
      IF(NVTPS.LE.MAXVTP) GO TO 1360
      ITEMP=MAXVTP/2
      PRINT 1350,ITEMP
 1350 FORMAT(' ***  ERROR  ***  T2D CAN HANDLE ONLY ',I3,' VTPS.')
      IERROR=IERROR+1
 1360 CONTINUE
      LSCR(1)=FNO
      LSCR(2)=LNO
      SCR(3)=OSI
      LSCR(4)=SDEPTH
      LSCR(5)=EDEPTH
      LSCR(6)=IADDWB
      LSCR(7)=LPRINT
      LSCR(8)=NVTPS
      LSCR(10)=VTYPE
      ITEMP=NPARS+1
      ITEMP1=NPARS+NVTPS
      IF(IAND(LPRINT,1).EQ.1)  PRINT *,(LSCR(I),I=1,2),
     *   SCR(3),(LSCR(J),J=4,8),LSCR(10),(SCR(J),J=ITEMP,ITEMP1)
      CALL WRDISC(MUNIT,SCR,NWRDS)
      NLISTS=NLISTS+1
      LLNO=LNO
      LNO=9999999                                                       ! DEFAULT THE DEFAULTS
      NS=0                                                              ! SET THE NUMBER OF MULTI-VALUED PARAMETER ENTRIES BACK TO ZER0
      VTP=-1.
 2020 CALL GETOKE(TOKEN,NCHARS)                                         ! GET THE NEXT TOKEN
      CALL UPCASE(TOKEN,NCHARS)
      NTOKES=NTOKES+1
      IF(NCHARS.GT.0) GO TO 2030                                        ! WAS IT THE END OF A LINE?
      IF(NOW.EQ.1) PRINT 140
      CALL RDLINE                                                       ! GET ANOTHER LINE
      NTOKES=0
      GO TO 2020
 2030 IF(TOKEN(1:NCHARS).NE.'END'.OR.NCHARS.NE.3) GO TO 150
      RETURN                                                            !  FINISHED ALL OF THE PARAMETERS!!!
      END

ext: f


text:
	SUBROUTINE HCLEAR  ( iret )
C************************************************************************
C* HCLEAR - GIF								*
C* 									*
C* This subroutine clears the current device.  On a direct access 	*
C* device, HCLEAR erases the screen.  On a continuous paper plotter, 	*
C* HCLEAR will advance  to the next page.  On a single page plotter, 	*
C* HCLEAR will unload the paper or place the plotter into pause mode so *
C* another sheet can be loaded.						*
C* 									*
C* HCLEAR  ( IRET )							*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* M. desJardins/GSFC	12/90						*
C* M. desJardins/NMC	 4/91	Check that something has been plotted	*
C* D. Austin		 5/96	Modified for GIF driver			*
C* T. Lee/GSC		 7/00	Renamed gdr_clear to wclear		*
C************************************************************************
	INCLUDE		'ERROR.PRM'
	INCLUDE		'driver.cmn'
C------------------------------------------------------------------------
	iret = NORMAL
C
	IF  ( ( .not. opnfil ) .or. ( .not. gfplot ) )  RETURN
C
C*	Clear the current device.
C
	CALL WCLEAR ( iret )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
c     ****************************************************************
c     *                                                              *
c     *                      subroutine setelb                       *
c     *                                                              *
c     *                       written by : bh                        *
c     *                                                              *
c     *                   last modified : 08/11/2017 rhd             *
c     *                                                              *
c     *     this subroutine sets the element library to include all  *
c     *     of the elements to be implemented by the program.        *
c     *                                                              *
c     ****************************************************************
c
c
c
      subroutine setelb( nlibel, elelib, outmap, mxlbel, out )
c
      implicit none
c
      integer :: nlibel, outmap(mxlbel,*), mxlbel, out
      character(len=8) :: elelib(mxlbel)
c
      integer :: i
c
      nlibel = 19
      if( nlibel > mxlbel ) then
         write(out,9000)
         call die_abort
      end if
c
      elelib(1:mxlbel) = ' '
c
c                       element type one: q3disop
c
      elelib(1) = 'q3disop '
      do i = 1, 30
        outmap(1,i) = i
      end do
c
c                       element type two: l3disop
c
      elelib(2) = 'l3disop '
      do i = 1, 30
        outmap(2,i) = i
      end do
c
c                       element type 3: ts12isop
c
      elelib(3) = 'ts12isop'
      do i = 1, 30
        outmap(3,i) = i
      end do
c
c                       element type 4: ts15isop
c
      elelib(4) = 'ts15isop'
      do i = 1, 30
        outmap(4,i) = i
      end do
c
c                       element type 5: ts9isop
c
      elelib(5) = 'ts9isop'
      do i = 1, 30
        outmap(5,i) = i
      end do
c
c                       element type 6: tet10
c
      elelib(6) = 'tet10'
      do i = 1, 30
        outmap(6,i) = i
      end do
c
c                       element type 7: wedge15
c
      elelib(7) = 'wedge15'
      do i = 1, 30
        outmap(7,i) = i
      end do
c
c                       element type 8: tri6
c
      elelib(8) = 'tri6'
      do i = 1, 30
        outmap(8,i) = i
      end do
c
c                       element type 9: quad8
c
      elelib(9) = 'quad8'
      do i = 1, 30
        outmap(9,i) = i
      end do
c
c                       element type 10: axiquad8
c
      elelib(10) = 'axiquad8'
      do i = 1, 30
        outmap(10,i) = i
      end do
c
c                       element type 11: axitri6
c
      elelib(11) = 'axitri6'
      do i = 1, 30
        outmap(11,i) = i
      end do
c
c                       element type 12: 'inter_8'
c
      elelib(12) = 'inter_8'
      do i = 1,30
        outmap(12,i) = i
      end do
c
c                       element type 13: tet4
c
      elelib(13) = 'tet4'
      do i = 1, 30
        outmap(13,i) = i
      end do
c
c                       element type 14: 'trint6'
c
      elelib(14) = 'trint6'
      do i = 1,30
        outmap(14,i) = i
      end do
c
c                       element type 15: 'trint12'
c
      elelib(15) = 'trint12'
      do i = 1,30
        outmap(15,i) = i
      end do
c
c                       element type 16: not a real element.
c                                        used to provide integration
c                                        points for 4-node quad
c
      elelib(16) = 'not_real'
      do i = 1,30
        outmap(16,i) = i
      end do
c
c                       element type 17: not a real element.
c                                        used to provide integration
c                                        points for 3 node line.
c                                        used in domain integral work
c
      elelib(17) = 'line3'
      do i = 1,30
        outmap(17,i) = i
      end do
c
c
c                       element type 18: 'bar2'
c
      elelib(18) = 'bar2'
      do i = 1,30
        outmap(18,i) = i
      end do
c
c
c                       element type 19: 'link2'
c
      elelib(19) = 'link2'
      do i = 1,30
        outmap(19,i) = i
      end do
c
      return
c
 9000 format(/,'>>>> FATAL ERROR: setelb. incorrect (nlibel,mxlebl)',
     & ' values',
     &       /,'                  job terminated',//)
      end


text:
!-----------------------------------------------------------------------
!
! MODULE: sweep_module
!> @brief
!> This module controls the mesh sweep scheduling. It directs the flow
!> of KBA stages for different octants, groups, spatial work chunks.
!
!-----------------------------------------------------------------------

MODULE sweep_module

  USE global_module, ONLY: i_knd, zero, one

  USE data_module, ONLY: ng

  USE geom_module, ONLY: ndiag

  USE control_module, ONLY: inrdone, swp_typ, ncor, corner_sch

  USE octsweep_module, ONLY: octsweep

  USE solvar_module, ONLY: flkx, flky, flkz, fmin, fmax

  USE plib_module, ONLY: nthreads, waitinit

  USE thrd_comm_module, ONLY: assign_thrd_set

  IMPLICIT NONE

  PRIVATE

  PUBLIC :: sweep

  SAVE

  CONTAINS

  SUBROUTINE sweep ( t, do_grp, ng_per_thrd, nnstd_used, grp_act )

!-----------------------------------------------------------------------
!
! Driver for the mesh sweeps. Manages the loops over octant pairs.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: t

    INTEGER(i_knd), INTENT(INOUT) :: ng_per_thrd, nnstd_used

    INTEGER(i_knd), DIMENSION(ng), INTENT(INOUT) :: do_grp

    INTEGER(i_knd), DIMENSION(ng,nthreads), INTENT(INOUT) :: grp_act

!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: corner, jd, kd, g, n

    INTEGER(i_knd), DIMENSION(2) :: reqs
!_______________________________________________________________________
!
!   Assign the work to threads. Main level threads always applied to
!   energy groups. Apply nested threads additionally to groups if
!   swp_typ is 0. Apply nested threads to mini-KBA if swp_typ is 1.
!_______________________________________________________________________

  !$OMP MASTER

    do_grp = 1
    WHERE ( inrdone ) do_grp = 0

    CALL assign_thrd_set ( do_grp, ng, ng_per_thrd, ndiag, nnstd_used, &
      grp_act )

  !$OMP END MASTER
  !$OMP BARRIER
!_______________________________________________________________________
!
!   Each thread initializes the reqs (send request) array for
!   asynchronous sends.
!_______________________________________________________________________

    CALL waitinit ( reqs, SIZE( reqs ) )
!_______________________________________________________________________
!
!   Clean up and initialize some values.
!_______________________________________________________________________

    clean_loop: DO n = 1, ng_per_thrd

      g = grp_act(n,t)
      IF ( g == 0 ) EXIT clean_loop

      fmin(g) = HUGE( one )
      fmax(g) = zero

      flkx(:,:,:,g) = zero
      flky(:,:,:,g) = zero
      flkz(:,:,:,g) = zero

    END DO clean_loop
!_______________________________________________________________________
!
!   Loop over octant pairs; set the starting corner, i.e., the direction
!   in y and z. Spawn the nested threads if used.
!_______________________________________________________________________

  !$OMP PARALLEL NUM_THREADS(nnstd_used) IF(nnstd_used>1)              &
  !$OMP& DEFAULT(SHARED) PRIVATE(corner,jd,kd,g) PROC_BIND(CLOSE)
    corner_loop: DO corner = 1, ncor

      jd = corner_sch(1,corner)
      kd = corner_sch(2,corner)
!_______________________________________________________________________
!
!     Loop over the groups assigned to each thread.
!_______________________________________________________________________

      grp_loop: DO n = 1, ng_per_thrd

        g = grp_act(n,t)
!_______________________________________________________________________
!
!       Sweep all the chunks of an octant pair (+/- x-dir).
!_______________________________________________________________________

        CALL octsweep ( g, jd, kd, t, reqs, SIZE( reqs ) )
!_______________________________________________________________________
!
!       End the loops. Destroy the task set.
!_______________________________________________________________________

      END DO grp_loop

    END DO corner_loop
  !$OMP END PARALLEL
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE sweep

END MODULE sweep_module

max_stars_repo_name: tomdeakin/SNAP

max_issues_repo_name: tomdeakin/SNAP

max_forks_repo_name: tomdeakin/SNAP


text:
 ! ==============================================================
 ! Copyright © 2020 Intel Corporation
 !
 ! SPDX-License-Identifier: MIT
 ! =============================================================
 !
 ! [DESCRIPTION]
 ! This program computes the integral (area under the curve) of a user-supplied
 ! function over an interval in a stepwise fashion. The interval is split into
 ! segments, and at each segment position the area of a rectangle is computed
 ! whose height is the value of sine at that point and the width is the segment
 ! width.  The areas of the rectangles are then summed.
 !
 ! The process is repeated with smaller and smaller width rectangles, more
 ! closely approximating the true value.
 !
 ! The source for this program also demonstrates recommended Fortran
 ! coding practices.
 !
 ! Compile the sample several times using different optimization options.
 !
 ! Read the Intel(R) Fortran Compiler Documentation for more information about these options.
 !
 ! Some of these automatic optimizations use features and options
 ! that can restrict program execution to specific architectures.
 !
 ! [COMPILE]
 ! Use the one of the following compiler options:
 !
 ! Windows*: /O1, /O2, /O3
 !
 ! Linux* and macOS*: -O1, -O2, -O3
 !

program int_sin
implicit none

! Create a value DP that is the "kind" number of a double precision value
! We will use this value in our declarations and constants.
integer, parameter :: DP = kind(0.0D0)

! Declare a named constant for pi, specifying the kind type
real(DP), parameter :: pi = 3.141592653589793238_DP

! Declare interval begin and end
real(DP), parameter :: interval_begin = 0.0_DP
real(DP), parameter :: interval_end   = 2.0_DP * pi

real(DP) :: step, sum, x_i
integer :: N, i, j
real clock_start, clock_finish

write (*,'(A)') "  "
write (*,'(A)') "    Number of    | Computed Integral |"
write (*,'(A)') " Interior Points |                   |"
call cpu_time (clock_start)

do j=2,26
  write (*,'(A)') "--------------------------------------"
  N = 2**j
  ! Compute stepsize for N-1 internal rectangles
  step = (interval_end - interval_begin) / real(N,DP);

  ! Approximate 1/2 area in first rectangle: f(x0) * (step/2)
  sum = INTEG_FUNC(interval_begin) * (step / 2.0_DP)

  do i=1,N-1
    x_i = real(i,DP) * step
    ! Apply midpoint rule:
    ! Given length = f(x), compute the area of the
    ! rectangle of width step
    sum = sum + (INTEG_FUNC(x_i) * step)
    end do

  ! Add approximate area in last rectangle for f(xN) * (step/2)
  sum = sum + (INTEG_FUNC(interval_end) * (step / 2.0_DP))

  write (*,'(T5,I10,T18,"|",2X,1P,E14.7,T38,"|")') N, sum
  end do

call cpu_time(clock_finish)
write (*,'(A)') "--------------------------------------"
write (*,'(A)') "  "
write (*,*) "CPU Time = ",(clock_finish - clock_start), " seconds"

contains

! Function to integrate
real(DP) function INTEG_FUNC (x)
  real(DP), intent(IN) :: x

  INTEG_FUNC = abs(sin(x))
  return
end function INTEG_FUNC

end program int_sin

max_forks_repo_name: junxnone/oneAPI-samples


text:
      program inter01

C     Check that insufficient calls do not cause core dump in
C     transformers and preconditions: it is now trapped by the
C     effects analysis

      call foo(i)

      print *, i

      end

      subroutine foo(i, j)
      i = 1
      j = 2
      end