text:
CM
C->>> -------------------------------------------> ems_se_i_ct_vr_rg <<<
c     Set the ranges for the integer control variables.
c
      subroutine ems_se_i_ct_vr_rg
      implicit none
      include 'EMSV.INC'
      include 'EMSMMGR.INC'
      include 'EMSPM.INC'
      include 'EMSFI.INC'
      include 'EMSMEM.INC'
CM      IF (emsol_xa .EQ. 1) THEN
C?      include 'CRASH.INC'
CM      ENDIF
      include 'ICTVR.INC'

      call ems_cp_i_a(n_ems_i_ct_vr, 0, i_ct_vr_lb, 0)
      call ems_cp_i_a(n_ems_i_ct_vr, 0, i_ct_vr_ub, 0)
      i_ct_vr_lb(ix_msg_fq) = 1
      i_ct_vr_ub(ix_msg_fq) = i_inf
      i_ct_vr_ub(ix_ems_wr_cn) = 99
      i_ct_vr_ub(ix_usr_mx_n_u) = 999
      i_ct_vr_ub(ix_n_si_it) = i_inf
      i_ct_vr_ub(ix_mx_n_si_it) = i_inf
      i_ct_vr_ub(ix_rsmi_msg_msk) = 31
      i_ct_vr_lb(ix_mx_n_r) = -mx_mx_ml_a_dim
      i_ct_vr_ub(ix_mx_n_r) =  mx_mx_ml_a_dim
      i_ct_vr_lb(ix_mx_n_c) = -mx_mx_ml_a_dim
      i_ct_vr_ub(ix_mx_n_c) =  mx_mx_ml_a_dim
      i_ct_vr_lb(ix_ml_nm_n_ch) = 8
      i_ct_vr_ub(ix_ml_nm_n_ch) = ml_nm_mx_n_ch
      i_ct_vr_ub(ix_dvx_mode) = pc_alg_sed
      i_ct_vr_lb(ix_it_usr_xit_fq) = 1
      i_ct_vr_ub(ix_it_usr_xit_fq) = i_inf
      i_ct_vr_ub(ix_n_el_in_aux_blk) = i_inf
      i_ct_vr_ub(ix_n_aux_blk) = i_inf
      i_ct_vr_lb(ix_n_pg_li) = 10
      i_ct_vr_ub(ix_n_pg_li) = i_inf
      i_ct_vr_ub(ix_n_r) = mx_n_r
      i_ct_vr_ub(ix_n_c) = mx_n_c
      i_ct_vr_ub(ix_n_pr_ifs) = n_r
      i_ct_vr_ub(ix_n_du_ifs) = n_c
      i_ct_vr_ub(ix_wr_ml_da_msk) = 1023
      i_ct_vr_ub(ix_wr_ml_st_msk) = 15 + 16 + 32
      i_ct_vr_lb(ix_mx_n_el_in_aux_blk) = 1
      i_ct_vr_ub(ix_mx_n_el_in_aux_blk) = i_inf
      i_ct_vr_lb(ix_wr_li_ln) = 60
      i_ct_vr_ub(ix_wr_li_ln) = i_inf
      i_ct_vr_ub(ix_mx_nd_n) = i_inf
      i_ct_vr_ub(ix_mx_n_nd) = i_inf
      i_ct_vr_ub(ix_mx_n_i_fs_sol) = i_inf
      i_ct_vr_ub(ix_n_i_fs_sol) = i_inf
      i_ct_vr_ub(ix_n_i_vr) = i_inf
      i_ct_vr_ub(ix_n_sos_se) = i_inf
      i_ct_vr_ub(ix_n_i_vr_no_i) = i_inf
      i_ct_vr_lb(ix_prob_st) = -1
      i_ct_vr_ub(ix_prob_st) =  6
      i_ct_vr_ub(ix_wr_rg_da_msk) = 1023
      i_ct_vr_lb(ix_mx_n_i_vr) = -i_inf
      i_ct_vr_ub(ix_mx_n_i_vr) = i_inf
      i_ct_vr_lb(ix_mx_n_sos_se) = -i_inf
      i_ct_vr_ub(ix_mx_n_sos_se) = i_inf
      i_ct_vr_ub(ix_i_alg) = 3
      i_ct_vr_lb(ix_mx_n_i_da_en) = -i_inf
      i_ct_vr_ub(ix_mx_n_i_da_en) = i_inf
c
CM      IF (emsol_xa .EQ. 1) THEN
C?      i_ct_vr_ub(ix_crsh_msk) = crsh_msk_ze_fs_bt +
C?     &     crsh_msk_ze_pr_act_bt
C?      i_ct_vr_lb(ix_crsh_tl_dl_pri) = -i_inf
C?      i_ct_vr_ub(ix_crsh_tl_dl_pri) = i_inf
C?      i_ct_vr_lb(ix_crsh_fn_dn) = -i_inf
C?      i_ct_vr_ub(ix_crsh_fn_dn) = i_inf
CM      ENDIF

      i_ct_vr_lb(ix_ems_er_cn) = -1
      i_ct_vr_ub(ix_ems_er_cn) = 99
      i_ct_vr_ub(ix_ts_parsmi) = i_inf

      i_ct_vr_lb(ix_ems_fi_wr_msk) = -1
      i_ct_vr_ub(ix_ems_fi_wr_msk) = 2
      i_ct_vr_ub(ix_scl_mode) = 1
      i_ct_vr_ub(ix_scl_bd) = i_inf
      i_ct_vr_lb(ix_mx_n_a_el) = -i_inf
      i_ct_vr_ub(ix_mx_n_a_el) = i_inf
      i_ct_vr_ub(ix_n_a_el) = i_inf
      i_ct_vr_ub(ix_ml_blk_st_msk) = bt16-1
      i_ct_vr_ub(ix_ml_da_st_msk) = bt16-1
      i_ct_vr_ub(ix_du_sol_mode) = mx_du_sol_mode
      i_ct_vr_ub(ix_ml_nm_n_rl) = ml_nm_mx_n_ch/8

      i_ct_vr_ub(ix_sslv_en_msk) = sslv_en_no_bd_cg + sslv_en_no_co_cg
      i_ct_vr_ub(ix_sslv_xit_msk) =
     &     sslv_xit_no_inv +
     &     sslv_xit_no_reset_non_bc_pr_act +
     &     sslv_xit_no_reset_bc_pr_act +
     &     sslv_xit_no_reset_non_bc_du_act
      i_ct_vr_ub(ix_n_bs) = i_inf
      i_ct_vr_ub(ix_n_reset) = i_inf

      i_ct_vr_ub(ix_u_pc) = 1
      i_ct_vr_ub(ix_r_pc) = 1
      i_ct_vr_ub(ix_eta_fi_mode_msk) =
     &     eta_fi_fwd_p_poss_bt +
     &     eta_fi_r_eta_poss_bt +
     &     eta_fi_bwd_p_poss_bt +
     &     eta_fi_fwd_p_y_bt +
     &     eta_fi_r_eta_y_bt +
     &     eta_fi_bwd_p_y_bt
      i_ct_vr_ub(ix_eta_fi_da_st_msk) =
     &     eta_fi_da_st_fwd_p + eta_fi_da_st_r_eta + eta_fi_da_st_bwd_p
      i_ct_vr_ub(ix_pc_alg) = 3
      i_ct_vr_ub(ix_n_dvx_it) = i_inf
      i_ct_vr_ub(ix_mn_n_dvx_it) = i_inf
      i_ct_vr_ub(ix_n_dvx_fwk) = i_inf
      i_ct_vr_lb(ix_r_du_act_sgn) = -1
      i_ct_vr_ub(ix_r_du_act_sgn) =  1

      i_ct_vr_lb(ix_inv_alg_msk) = inv_alg_tom
      i_ct_vr_ub(ix_inv_alg_msk) = inv_alg_tom + inv_alg_perm +
     &     inv_alg_mwz + inv_alg_sus
      i_ct_vr_ub(ix_inv_msg_msk) = 7
      i_ct_vr_ub(ix_inv_log_msk) = 7
      i_ct_vr_ub(ix_inv_mx_n_eta) = i_inf
      i_ct_vr_ub(ix_inv_mx_n_eta_el) = i_inf
      i_ct_vr_ub(ix_n_inv_sing) = i_inf
      i_ct_vr_ub(ix_u_bs_cg) = 1
      i_ct_vr_ub(ix_rq_inv) = 1

      i_ct_vr_lb(ix_sto_ftran_ix_mode) = -1
      i_ct_vr_ub(ix_sto_ftran_ix_mode) =  1
      i_ct_vr_lb(ix_sto_btran_ix_mode) = -1
      i_ct_vr_ub(ix_sto_btran_ix_mode) =  1
      i_ct_vr_lb(ix_sto_tbu_r_ix_mode) = -1
      i_ct_vr_ub(ix_sto_tbu_r_ix_mode) =  1
      i_ct_vr_lb(ix_tbu_r_loop_mode) =   -1
      i_ct_vr_ub(ix_tbu_r_loop_mode) =    1
      i_ct_vr_lb(ix_sto_ftran_ix) = -1
      i_ct_vr_ub(ix_sto_ftran_ix) =  1
      i_ct_vr_lb(ix_sto_btran_ix) = -1
      i_ct_vr_ub(ix_sto_btran_ix) =  1
      i_ct_vr_lb(ix_sto_tbu_r_ix) = -1
      i_ct_vr_ub(ix_sto_tbu_r_ix) =  1
      i_ct_vr_ub(ix_bwd_tran_dse_rhs_n_r) = i_inf
      i_ct_vr_ub(ix_fwd_tran_dse_rhs_n_r) = i_inf
      i_ct_vr_ub(ix_dse_tbu_r_n_c) = i_inf

      i_ct_vr_ub(ix_usr_cz_c_msk) = cz_c_bk_bd_bt
      i_ct_vr_ub(ix_usr_cz_r_msk) =
     &     cz_r_l1_bt + cz_r_refine_bt + cz_r_growth_mode
      i_ct_vr_lb(ix_l1_cz_r_mx_n_cdd) = 1
      i_ct_vr_ub(ix_l1_cz_r_mx_n_cdd) = i_inf
      i_ct_vr_ub(ix_n_pc_vr) = i_inf
      i_ct_vr_ub(ix_n_pc_el) = i_inf
      i_ct_vr_ub(ix_n_non_bc_pr_ifs) = i_inf
      i_ct_vr_ub(ix_rq_re_pc) = i_inf
      i_ct_vr_ub(ix_cz_c_msk) = i_ct_vr_ub(ix_usr_cz_c_msk)
      i_ct_vr_ub(ix_cz_r_msk) = i_ct_vr_ub(ix_usr_cz_r_msk)

      i_ct_vr_ub(ix_u_bs_mode) = 2
      i_ct_vr_ub(ix_sto_u_eta_v_mode) = 2
      i_ct_vr_ub(ix_u_bs) = 2
      i_ct_vr_ub(ix_sto_u_eta_v) = 1
      i_ct_vr_ub(ix_sto_u_eta_se_ty) = n_eta_se_ty
      i_ct_vr_ub(ix_mx_n_u) = 999
      i_ct_vr_ub(ix_n_u) = i_inf

      i_ct_vr_ub(ix_lo_sto_mode) = 1
      i_ct_vr_ub(ix_asm_msk) = 15
      i_ct_vr_ub(ix_ck_msk) = 255
      i_ct_vr_ub(ix_it_o_sw_fm_wr_df) = i_inf
      i_ct_vr_ub(ix_it_o_sw_fm_ck_df) = i_inf
      i_ct_vr_ub(ix_wr_lp_da) = 1023
      i_ct_vr_ub(ix_wr_bs_fq) = i_inf

      i_ct_vr_lb(ix_g_tt_da) = -i_inf
      i_ct_vr_ub(ix_g_tt_da) = bt15-1
      i_ct_vr_ub(ix_g_it_tt_fq) = i_inf
      i_ct_vr_ub(ix_g_it_tt_ivl) = i_inf

      i_ct_vr_ub(ix_mx_n_bc_pr_act_rfn_it) = i_inf
      i_ct_vr_ub(ix_mx_n_reset_pi_rfn_it) = i_inf
      i_ct_vr_ub(ix_mx_n_pv_c_rfn_it) = i_inf
      i_ct_vr_ub(ix_mx_n_full_pi_rfn_it) = i_inf
      i_ct_vr_ub(ix_mx_n_u_pi_rfn_it) = i_inf

      i_ct_vr_ub(ix_eta_fi_n_grp) = eta_fi_mx_n_eta_grp
      i_ct_vr_ub(ix_eta_fi_n_inv_grp) = eta_fi_mx_n_eta_grp
      i_ct_vr_ub(ix_eta_fi_n_se) = i_inf
      i_ct_vr_ub(ix_eta_fi_n_inv_se) = i_inf
      i_ct_vr_ub(ix_eta_fi_n_eta) = i_inf
      i_ct_vr_ub(ix_eta_fi_n_inv_eta) = i_inf
      i_ct_vr_ub(ix_eta_fi_n_v) = i_inf
      i_ct_vr_ub(ix_eta_fi_n_inv_v) = i_inf
      i_ct_vr_ub(ix_eta_fi_n_ix) = i_inf
      i_ct_vr_ub(ix_eta_fi_n_inv_ix) = i_inf

      i_ct_vr_lb(ix_mx_n_pwl_vr) = -i_inf
      i_ct_vr_ub(ix_mx_n_pwl_vr) = i_inf
      i_ct_vr_lb(ix_mx_n_pwl_vr_da_en) = -i_inf
      i_ct_vr_ub(ix_mx_n_pwl_vr_da_en) = i_inf
      i_ct_vr_ub(ix_repl_non_std_vr_mode_msk) =
     &     repl_non_std_vr_bp_bt + repl_non_std_vr_pwl_bt
      i_ct_vr_ub(ix_n_bp_vr) = i_inf
      i_ct_vr_ub(ix_n_pwl_vr) = i_inf
      i_ct_vr_ub(ix_n_pwl_vr_da_en) = i_inf
c
      i_ct_vr_ub(ix_cu_is_n) = i_inf
      i_ct_vr_ub(ix_cu_ml_n) = i_inf
      i_ct_vr_ub(ix_ml_blk_mv_k) = i_inf
      i_ct_vr_ub(ix_ds_n_en_m1) = i_inf
      i_ct_vr_ub(ix_is_n_en_m1) = i_inf
      i_ct_vr_ub(ix_sv_ml_ct_vr) = 1

      i_ct_vr_ub(ix_n_lo_c_eta) = i_inf
      i_ct_vr_ub(ix_n_lo_r_eta) = i_inf
      i_ct_vr_ub(ix_n_lo_r_eta_el) = i_inf
      i_ct_vr_ub(ix_n_up_r_eta) = i_inf
      i_ct_vr_ub(ix_n_up_r_eta_el) = i_inf
      i_ct_vr_ub(ix_n_sus_eta_srch_op_cf) = i_inf
      i_ct_vr_ub(ix_sus_eta_srch_al_t_ls_ix_lm) = i_inf
      i_ct_vr_ub(ix_sus_eta_srch_ls_t_buk_ix_lm) = i_inf
      i_ct_vr_ub(ix_sus_eta_srch_dbl_n_buk_ix_lm) = i_inf
      return
      end
C->>> -------------------------------------------> ems_iz_i_ct_vr_df <<<
c     Initialise the settable logical and default values for the
c     integer control variables.
c
      subroutine ems_iz_i_ct_vr_df
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'ICTVR.INC'

      call ems_cp_i_a(n_ems_i_ct_vr, 0, i_ct_vr_df, 0)
      i_ct_vr_df(ix_msg_fq) = 999999
      i_ct_vr_df(ix_ems_wr_cn) = 6
      i_ct_vr_df(ix_usr_mx_n_u) = 25
      i_ct_vr_df(ix_n_si_it) = 0
      i_ct_vr_df(ix_mx_n_si_it) = 999999
      i_ct_vr_df(ix_rsmi_msg_msk) = 16
      i_ct_vr_df(ix_mx_n_r) = 0
      i_ct_vr_df(ix_mx_n_c) = 0
      i_ct_vr_df(ix_ml_nm_n_ch) = 8
      i_ct_vr_df(ix_dvx_mode) = pc_alg_approx_dvx
      i_ct_vr_df(ix_it_usr_xit_fq) = 1
      i_ct_vr_df(ix_n_el_in_aux_blk) = 0
      i_ct_vr_df(ix_n_aux_blk) = 0
      i_ct_vr_df(ix_n_pg_li) = 999999
      i_ct_vr_df(ix_n_r) = 0
      i_ct_vr_df(ix_n_c) = 0
      i_ct_vr_df(ix_n_pr_ifs) = 0
      i_ct_vr_df(ix_n_du_ifs) = 0
      i_ct_vr_df(ix_wr_ml_da_msk) = 3
      i_ct_vr_df(ix_wr_ml_st_msk) = 0
      i_ct_vr_df(ix_mx_n_el_in_aux_blk) = 1000
      i_ct_vr_df(ix_wr_li_ln) = 80
      i_ct_vr_df(ix_mx_n_nd) = 9999999
      i_ct_vr_df(ix_mx_n_i_fs_sol) = 9999999
      i_ct_vr_df(ix_prob_st) = -1
      i_ct_vr_df(ix_wr_rg_da_msk) = 0
      i_ct_vr_df(ix_i_alg) = 1

CM      IF (emsol_xa .EQ. 1) THEN
C?      i_ct_vr_df(ix_crsh_msk) = 0
C?      i_ct_vr_df(ix_crsh_tl_dl_pri) = -10
C?      i_ct_vr_df(ix_crsh_fn_dn) = 10
CM      ENDIF

CM      IF (emsol_xa .EQ. 1) THEN
C?      i_ct_vr_df(ix_ems_fi_wr_msk) = 0
CM      ELSE
      i_ct_vr_df(ix_ems_fi_wr_msk) = 1
CM      ENDIF

      i_ct_vr_df(ix_scl_mode) = 0
      i_ct_vr_df(ix_scl_bd) = 20
      i_ct_vr_df(ix_mx_n_a_el) = -1000
      i_ct_vr_df(ix_mx_n_a_el) = 0
      i_ct_vr_df(ix_ml_nm_n_rl) = 1

      i_ct_vr_df(ix_u_pc) = 1
      i_ct_vr_df(ix_r_pc) = 1
      i_ct_vr_df(ix_eta_fi_mode_msk) =
     &     eta_fi_fwd_p_poss_bt +
     &     eta_fi_r_eta_poss_bt +
     &     eta_fi_bwd_p_poss_bt
      i_ct_vr_df(ix_mn_n_dvx_it) = 25
      i_ct_vr_df(ix_r_du_act_sgn) = 1

      i_ct_vr_df(ix_inv_alg_msk) = inv_alg_tom
     &     + inv_alg_sus
      i_ct_vr_df(ix_rq_inv) = 1

      i_ct_vr_df(ix_sto_ftran_ix_mode) = sto_ix_poss
      i_ct_vr_df(ix_sto_btran_ix_mode) = sto_ix_poss
      i_ct_vr_df(ix_sto_tbu_r_ix_mode) = sto_ix_poss
      i_ct_vr_df(ix_tbu_r_loop_mode) = tbu_r_loop_poss
      i_ct_vr_df(ix_sto_ftran_ix) = i_ct_vr_df(ix_sto_ftran_ix_mode)
      i_ct_vr_df(ix_sto_btran_ix) = i_ct_vr_df(ix_sto_btran_ix_mode)
      i_ct_vr_df(ix_sto_tbu_r_ix) = i_ct_vr_df(ix_sto_tbu_r_ix_mode)

      i_ct_vr_df(ix_usr_cz_c_msk) = 0
      i_ct_vr_df(ix_usr_cz_r_msk) = cz_r_refine_bt + cz_r_growth_inv
      i_ct_vr_df(ix_l1_cz_r_mx_n_cdd) = df_l1_cz_r_mx_n_cdd
      i_ct_vr_df(ix_cz_c_msk) = i_ct_vr_df(ix_usr_cz_c_msk)
      i_ct_vr_df(ix_cz_r_msk) = i_ct_vr_df(ix_usr_cz_r_msk)

c      i_ct_vr_df(ix_u_bs_mode) = u_bs_pf
      i_ct_vr_df(ix_u_bs_mode) = u_bs_pf_r_cp
      i_ct_vr_df(ix_sto_u_eta_v_mode) = sto_pk_eta_v
      i_ct_vr_df(ix_sto_u_eta_v) = i_ct_vr_df(ix_sto_u_eta_v_mode)
      i_ct_vr_df(ix_u_bs) = u_bs_pf
      i_ct_vr_df(ix_sto_u_eta_se_ty) = pk_c_eta_se_ty

      i_ct_vr_df(ix_g_tt_da) = 0

      i_ct_vr_df(ix_mx_n_bc_pr_act_rfn_it) = 3
      i_ct_vr_df(ix_mx_n_reset_pi_rfn_it) = 3

      i_ct_vr_df(ix_mx_n_pwl_vr_da_en) = 0
c
c     0 => Replace nothing
c     1 => Replace BP
c     2 => Replace PWL
c
c      i_ct_vr_df(ix_repl_non_std_vr_mode_msk) = 0
      i_ct_vr_df(ix_repl_non_std_vr_mode_msk) = repl_non_std_vr_pwl_bt

      i_ct_vr_df(ix_cu_is_n) = 1
      i_ct_vr_df(ix_cu_ml_n) = 1
      i_ct_vr_df(ix_sv_ml_ct_vr) = 1

      i_ct_vr_df(ix_n_sus_eta_srch_op_cf) = 1

      call ems_cp_lg_a(n_ems_i_ct_vr, .false., i_ct_vr_lg, 0)
      i_ct_vr_lg(ix_msg_fq) = .true.
      i_ct_vr_lg(ix_ems_wr_cn) = .true.
      i_ct_vr_lg(ix_usr_mx_n_u) = .true.
      i_ct_vr_lg(ix_mx_n_si_it) = .true.
      i_ct_vr_lg(ix_rsmi_msg_msk) = .true.
      i_ct_vr_lg(ix_mx_n_r) = .true.
      i_ct_vr_lg(ix_mx_n_c) = .true.
      i_ct_vr_lg(ix_ml_nm_n_ch) = .true.
      i_ct_vr_lg(ix_dvx_mode) = .true.
      i_ct_vr_lg(ix_it_usr_xit_fq) = .true.
      i_ct_vr_lg(ix_n_pg_li) = .true.
      i_ct_vr_lg(ix_wr_ml_da_msk) = .true.
      i_ct_vr_lg(ix_wr_ml_st_msk) = .true.
c      i_ct_vr_lg(ix_mx_n_el_in_aux_blk) = .true.
      i_ct_vr_lg(ix_wr_li_ln) = .true.
      i_ct_vr_lg(ix_wr_rg_da_msk) = .true.
      i_ct_vr_lg(ix_mx_n_nd) = .true.
      i_ct_vr_lg(ix_mx_n_i_fs_sol) = .true.
      i_ct_vr_lg(ix_mx_n_i_vr) = .true.
      i_ct_vr_lg(ix_mx_n_sos_se) = .true.
      i_ct_vr_lg(ix_mx_n_i_da_en) = .true.
c
c     Made user-accessible on 21/11/97
c
      i_ct_vr_lg(ix_mx_n_a_el) = .true.
      i_ct_vr_lg(ix_ems_fi_wr_msk) = .true.
      i_ct_vr_lg(ix_repl_non_std_vr_mode_msk) = .true.

CM      IF (emsol_dev .EQ. 1) THEN
C?      i_ct_vr_lg(ix_crsh_msk) = .true.
C?      i_ct_vr_lg(ix_crsh_tl_dl_pri) = .true.
C?      i_ct_vr_lg(ix_crsh_fn_dn) = .true.
C?
C?      i_ct_vr_lg(ix_ems_er_cn) = .true.
C?      i_ct_vr_lg(ix_ts_parsmi) = .true.
C?
C?      i_ct_vr_lg(ix_scl_mode) = .true.
C?c      i_ct_vr_lg(ix_scl_bd) = .true.
C?
C?      i_ct_vr_lg(ix_sslv_en_msk) = .true.
C?      i_ct_vr_lg(ix_sslv_xit_msk) = .true.
C?
C?      i_ct_vr_lg(ix_u_pc) = .true.
C?      i_ct_vr_lg(ix_mn_n_dvx_it) = .true.
C?
C?      i_ct_vr_lg(ix_inv_msg_msk) = .true.
C?      i_ct_vr_lg(ix_inv_log_msk) = .true.
C?
C?      i_ct_vr_lg(ix_r_pc) = .true.
C?      i_ct_vr_lg(ix_eta_fi_mode_msk) = .true.
C?      i_ct_vr_lg(ix_inv_alg_msk) = .true.
C?      i_ct_vr_lg(ix_sto_btran_ix_mode) = .true.
C?      i_ct_vr_lg(ix_sto_ftran_ix_mode) = .true.
C?      i_ct_vr_lg(ix_sto_tbu_r_ix_mode) = .true.
C?      i_ct_vr_lg(ix_tbu_r_loop_mode) = .true.
C?      i_ct_vr_lg(ix_usr_cz_r_msk) = .true.
C?      i_ct_vr_lg(ix_mx_n_bc_pr_act_rfn_it) = .true.
C?      i_ct_vr_lg(ix_mx_n_reset_pi_rfn_it) = .true.
C?
C?      i_ct_vr_lg(ix_usr_cz_c_msk) = .true.
C?      i_ct_vr_lg(ix_usr_cz_r_msk) = .true.
C?      i_ct_vr_lg(ix_l1_cz_r_mx_n_cdd) = .true.
C?
C?      i_ct_vr_lg(ix_u_bs_mode) = .true.
C?c      i_ct_vr_lg(ix_sto_u_eta_v_mode) = .true.
C?
C?      i_ct_vr_lg(ix_lo_sto_mode) = .true.
C?      i_ct_vr_lg(ix_asm_msk) = .true.
C?      i_ct_vr_lg(ix_ck_msk) = .true.
C?      i_ct_vr_lg(ix_it_o_sw_fm_wr_df) = .true.
C?      i_ct_vr_lg(ix_it_o_sw_fm_ck_df) = .true.
C?      i_ct_vr_lg(ix_wr_lp_da) = .true.
C?      i_ct_vr_lg(ix_wr_bs_fq) = .true.
C?
C?      i_ct_vr_lg(ix_g_tt_da) = .true.
C?      i_ct_vr_lg(ix_g_it_tt_fq) = .true.
C?      i_ct_vr_lg(ix_g_it_tt_ivl) = .true.
C?
C?      i_ct_vr_lg(ix_mx_n_pv_c_rfn_it) = .true.
C?      i_ct_vr_lg(ix_mx_n_full_pi_rfn_it) = .true.
C?      i_ct_vr_lg(ix_mx_n_u_pi_rfn_it) = .true.
C?
C?      i_ct_vr_lg(ix_mx_n_pwl_vr) = .true.
C?      i_ct_vr_lg(ix_mx_n_pwl_vr_da_en) = .true.
CM      ENDIF
      i_ct_vr_lg(ix_n_sus_eta_srch_op_cf) = .true.
      i_ct_vr_lg(ix_sus_eta_srch_al_t_ls_ix_lm) = .true.
      i_ct_vr_lg(ix_sus_eta_srch_ls_t_buk_ix_lm) = .true.
      i_ct_vr_lg(ix_sus_eta_srch_dbl_n_buk_ix_lm) = .true.
      return
      end

C->>> ----------------------------------------------> ems_wr_i_ct_vr <<<
      subroutine ems_wr_i_ct_vr
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'ICTVR.INC'
      include 'EMSMSG.INC'
      integer ix
      do 10, ix = 1, n_ems_i_ct_vr
         if (i_ct_vr_lg(ix)) then
            if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9000)
     &           ix, i_ct_vr_lb(ix), i_ct_vr_ub(ix),
     &           i_ct_vr_df(ix), i_ct_vr(ix)
            call ems_msg_wr_li(info_msg_n)
         else
            if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9010)
     &           ix, i_ct_vr_lb(ix), i_ct_vr_ub(ix),
     &           i_ct_vr_df(ix), i_ct_vr(ix)
            call ems_msg_wr_li(info_msg_n)
         end if
         if (i_ct_vr(ix) .lt. i_ct_vr_lb(ix) .or.
     &        i_ct_vr(ix) .gt. i_ct_vr_ub(ix)) then
            if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9100)ix
            call ems_msg_wr_li(er_msg_n)
         end if
 10   continue
      return
 9000 format('i_ct_vr ', i3, '     is settable: range [', i9, ', ',
     &     i9, ']; default ', i9, '; value ', i9)
 9010 format('i_ct_vr ', i3, ' is not settable: range [', i9, ', ',
     &     i9, ']; default ', i9, '; value ', i9)
 9100 format('i_ct_vr ', i3, ' violates its range ')
      end

C->>> ------------------------------------------> ems_se_rl_ct_vr_rg <<<
c     Set the ranges for the real control variables.
c
      subroutine ems_se_rl_ct_vr_rg
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'RLCTVR.INC'

      call ems_cp_rl_a(n_ems_rl_ct_vr, zero, rl_ct_vr_lb, 0)
      call ems_cp_rl_a(n_ems_rl_ct_vr, zero, rl_ct_vr_ub, 0)
      rl_ct_vr_lb(ix_usr_tl_pr_ifs) = 1d-12
      rl_ct_vr_ub(ix_usr_tl_pr_ifs) = 1d-1
      rl_ct_vr_lb(ix_usr_tl_du_ifs) = 1d-12
      rl_ct_vr_ub(ix_usr_tl_du_ifs) = 1d-1
      rl_ct_vr_lb(ix_mx_mn) = -1d0
      rl_ct_vr_ub(ix_mx_mn) =  1d0
      rl_ct_vr_lb(ix_pr_wt) = zero
      rl_ct_vr_ub(ix_pr_wt) = 1d+10
      rl_ct_vr_lb(ix_bb_ub) = -1d20
      rl_ct_vr_ub(ix_bb_ub) = undn
      rl_ct_vr_lb(ix_ob_fn_v) = -undn
      rl_ct_vr_ub(ix_ob_fn_v) =  undn
      rl_ct_vr_lb(ix_su_pr_ifs) = -undn
      rl_ct_vr_ub(ix_su_pr_ifs) =  undn
      rl_ct_vr_lb(ix_su_du_ifs) = -undn
      rl_ct_vr_ub(ix_su_du_ifs) =  undn
      rl_ct_vr_lb(ix_mps_ze) = -undn
      rl_ct_vr_ub(ix_mps_ze) =  undn
      rl_ct_vr_ub(ix_dga_scl) = undn
      rl_ct_vr_lb(ix_bst_i_fs_ob_fn_v) = -undn
      rl_ct_vr_ub(ix_bst_i_fs_ob_fn_v) = undn
      rl_ct_vr_ub(ix_i_ifs_wt) = undn
      rl_ct_vr_lb(ix_tl_i_fs_ob_fn_v_imp) = -undn
      rl_ct_vr_ub(ix_tl_i_fs_ob_fn_v_imp) = undn
      rl_ct_vr_lb(ix_i_fs_ob_fn_v_tgt) = -undn
      rl_ct_vr_ub(ix_i_fs_ob_fn_v_tgt) = undn
      rl_ct_vr_lb(ix_tl_i_ifs) = 1d-12
      rl_ct_vr_ub(ix_tl_i_ifs) = 1d-1
      rl_ct_vr_lb(ix_bst_ob_fn_v) = -undn
      rl_ct_vr_ub(ix_bst_ob_fn_v) = undn
      rl_ct_vr_lb(ix_bst_est_ob_fn_v) = -undn
      rl_ct_vr_ub(ix_bst_est_ob_fn_v) = undn

      rl_ct_vr_lb(ix_wk_tl_pr_ifs) = 1d-12
      rl_ct_vr_ub(ix_wk_tl_pr_ifs) = 1d-1
c
c     ?? Temporary change to allow experiments
c
      rl_ct_vr_lb(ix_wk_tl_pr_ifs) = 0d0
      rl_ct_vr_ub(ix_wk_tl_pr_ifs) = inf

      rl_ct_vr_lb(ix_wk_tl_du_ifs) = 1d-12
      rl_ct_vr_ub(ix_wk_tl_du_ifs) = 1d-1
      rl_ct_vr_ub(ix_wk_xp_tau) = 1d-4
      rl_ct_vr_lb(ix_tl_pr_ifs) = 1d-12
      rl_ct_vr_ub(ix_tl_pr_ifs) = 1d-1
      rl_ct_vr_lb(ix_tl_du_ifs) = 1d-12
      rl_ct_vr_ub(ix_tl_du_ifs) = 1d-1
      rl_ct_vr_ub(ix_xp_tau) = 1d-4
      rl_ct_vr_lb(ix_tl_bwd_tran_dse_rhs) = zero
      rl_ct_vr_ub(ix_tl_bwd_tran_dse_rhs) = one
      rl_ct_vr_lb(ix_tl_fwd_tran_dse_rhs) = zero
      rl_ct_vr_ub(ix_tl_fwd_tran_dse_rhs) = one
      rl_ct_vr_lb(ix_tl_dse_tbu_r) = zero
      rl_ct_vr_ub(ix_tl_dse_tbu_r) = one

      rl_ct_vr_lb(ix_usr_tl_crsh_abs_pv_v) = 1d-16
      rl_ct_vr_ub(ix_usr_tl_crsh_abs_pv_v) = inf
      rl_ct_vr_lb(ix_usr_tl_crsh_rlv_pv_v) = 1d-16
      rl_ct_vr_ub(ix_usr_tl_crsh_rlv_pv_v) = one
      rl_ct_vr_lb(ix_tl_iz_bs_tran_er) = zero
      rl_ct_vr_ub(ix_tl_iz_bs_tran_er) = inf
      rl_ct_vr_lb(ix_tl_mx_iz_pr_act) = zero
      rl_ct_vr_ub(ix_tl_mx_iz_pr_act) = inf

      rl_ct_vr_ub(ix_tl_du_act_er) = undn
      rl_ct_vr_ub(ix_tl_pv_er) = undn
      rl_ct_vr_ub(ix_tl_ed_wt_er) = undn
      rl_ct_vr_ub(ix_tl_tran_er) = undn
      rl_ct_vr_lb(ix_tl_cz_r_growth) = one
      rl_ct_vr_ub(ix_tl_cz_r_growth) = undn
      rl_ct_vr_ub(ix_fac_long_prd) = one
      rl_ct_vr_ub(ix_fac_short_prd) = one

      rl_ct_vr_ub(ix_bwd_tran_ze) = undn
      rl_ct_vr_ub(ix_pc_ze) = undn
      rl_ct_vr_ub(ix_fwd_tran_ze) = undn
      rl_ct_vr_ub(ix_u_ed_wt_ze) = undn
      rl_ct_vr_ub(ix_pk_pv_c_ze) = undn
      rl_ct_vr_ub(ix_scl_v_ze) = undn

      rl_ct_vr_lb(ix_av_eta_dse) = -undn
      rl_ct_vr_ub(ix_av_eta_dse) =  undn
      rl_ct_vr_lb(ix_av_ftran_sol_dse) = -undn
      rl_ct_vr_ub(ix_av_ftran_sol_dse) =  undn
      rl_ct_vr_lb(ix_av_btran_sol_dse) = -undn
      rl_ct_vr_ub(ix_av_btran_sol_dse) =  undn
      rl_ct_vr_lb(ix_av_tbu_r_dse) = -undn
      rl_ct_vr_ub(ix_av_tbu_r_dse) =  undn
      rl_ct_vr_ub(ix_tl_dvx_wt) =  undn
      rl_ct_vr_ub(ix_nw_dvx_fwk_fq) = undn

      rl_ct_vr_lb(ix_inv_pv_tl) =     1d-16
      rl_ct_vr_ub(ix_inv_pv_tl) =     inf
      rl_ct_vr_lb(ix_inv_wr_eta_tl) = 1d-19
      rl_ct_vr_ub(ix_inv_wr_eta_tl) = inf
      rl_ct_vr_lb(ix_inv_c_rlv_tl) =  1d-6
      rl_ct_vr_ub(ix_inv_c_rlv_tl) =  1d0
      rl_ct_vr_lb(ix_inv_unit_v_tl) = 0d0
      rl_ct_vr_ub(ix_inv_unit_v_tl) = 1d0
      rl_ct_vr_lb(ix_inv_l_pv_rsdu_er_tl) = 0d0
      rl_ct_vr_ub(ix_inv_l_pv_rsdu_er_tl) = 1d0

      rl_ct_vr_lb(ix_usr_cz_r_pv_tl) = 0d0
      rl_ct_vr_ub(ix_usr_cz_r_pv_tl) = 1d0
      rl_ct_vr_lb(ix_wk_cz_r_pv_tl) =  rl_ct_vr_lb(ix_usr_cz_r_pv_tl)
      rl_ct_vr_ub(ix_wk_cz_r_pv_tl) =  rl_ct_vr_ub(ix_usr_cz_r_pv_tl)
      rl_ct_vr_lb(ix_cz_r_pv_tl) =     rl_ct_vr_lb(ix_usr_cz_r_pv_tl)
      rl_ct_vr_ub(ix_cz_r_pv_tl) =     rl_ct_vr_ub(ix_usr_cz_r_pv_tl)

      rl_ct_vr_ub(ix_tl_bc_pr_act_it_rfn) = undn
      rl_ct_vr_ub(ix_tl_reset_pi_it_rfn) = undn
      rl_ct_vr_ub(ix_tl_pv_c_it_rfn) = undn
      rl_ct_vr_ub(ix_tl_full_pi_it_rfn) = undn
      rl_ct_vr_ub(ix_tl_u_pi_it_rfn) = undn
      rl_ct_vr_ub(ix_bc_pr_act_it_rfn_tran_ze) = undn
      rl_ct_vr_ub(ix_reset_pi_it_rfn_tran_ze) = undn
      rl_ct_vr_ub(ix_pv_c_it_rfn_tran_ze) = undn
      rl_ct_vr_ub(ix_full_pi_it_rfn_tran_ze) = undn
      rl_ct_vr_ub(ix_u_pi_it_rfn_tran_ze) = undn

      rl_ct_vr_lb(ix_ob_fn_cs) = -undn
      rl_ct_vr_ub(ix_ob_fn_cs) = undn
      rl_ct_vr_ub(ix_su_non_bc_pr_ifs) = undn

      return
      end

C->>> ------------------------------------------> ems_iz_rl_ct_vr_df <<<
c     Initialise the settable logical and default values for the
c     real control variables.
c
      subroutine ems_iz_rl_ct_vr_df
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'ICTVR.INC'
      include 'RLCTVR.INC'

      call ems_cp_rl_a(n_ems_rl_ct_vr, zero, rl_ct_vr_df, 0)
      rl_ct_vr_df(ix_usr_tl_pr_ifs) = 1d-6
      rl_ct_vr_df(ix_usr_tl_du_ifs) = 1d-6
      rl_ct_vr_df(ix_mx_mn) =         one
      rl_ct_vr_df(ix_pr_wt) =         zero
      rl_ct_vr_df(ix_bb_ub) =         inf
      rl_ct_vr_df(ix_dga_scl) =       one
      rl_ct_vr_df(ix_mps_ze) =        1d-12
      rl_ct_vr_df(ix_bst_i_fs_ob_fn_v) = inf
      rl_ct_vr_df(ix_i_ifs_wt) =      one
      rl_ct_vr_df(ix_tl_i_fs_ob_fn_v_imp) = 1d-5
      rl_ct_vr_df(ix_tl_i_ifs) = 1d-6

      rl_ct_vr_df(ix_wk_tl_pr_ifs) =   1d-6
      rl_ct_vr_df(ix_wk_tl_du_ifs) =   1d-6
      rl_ct_vr_df(ix_wk_xp_tau) =      0d0
      rl_ct_vr_df(ix_tl_pr_ifs) =      rl_ct_vr_df(ix_wk_tl_pr_ifs)
      rl_ct_vr_df(ix_tl_du_ifs) =      rl_ct_vr_df(ix_wk_tl_du_ifs)
      rl_ct_vr_df(ix_xp_tau) =         rl_ct_vr_df(ix_wk_xp_tau)
      rl_ct_vr_df(ix_tl_bwd_tran_dse_rhs) = 4d-1
      rl_ct_vr_df(ix_tl_fwd_tran_dse_rhs) = 1d-1
      rl_ct_vr_df(ix_tl_dse_tbu_r) =        1d-1

      rl_ct_vr_df(ix_usr_tl_crsh_abs_pv_v) =      1d-4
      rl_ct_vr_df(ix_usr_tl_crsh_rlv_pv_v) =      1d-2
      rl_ct_vr_df(ix_tl_iz_bs_tran_er) =          1d-10
      rl_ct_vr_df(ix_tl_mx_iz_pr_act) = 1d+10

      rl_ct_vr_df(ix_tl_du_act_er) =     1d-6
      rl_ct_vr_df(ix_tl_pv_er) =         1d-6
      rl_ct_vr_df(ix_tl_ed_wt_er) =      1d-2
      rl_ct_vr_df(ix_tl_tran_er) =       1d-6
      rl_ct_vr_df(ix_tl_cz_r_growth) =   1d8
      rl_ct_vr_df(ix_fac_long_prd) =     5d-2
      rl_ct_vr_df(ix_fac_short_prd) =    2d-1

      rl_ct_vr_df(ix_bwd_tran_ze) =         1d-12
      rl_ct_vr_df(ix_pc_ze) =            1d-12
      rl_ct_vr_df(ix_fwd_tran_ze) =         1d-12
      rl_ct_vr_df(ix_u_ed_wt_ze) =       1d-12
      rl_ct_vr_df(ix_pk_pv_c_ze) =       1d-12
      rl_ct_vr_ub(ix_scl_v_ze) =         0d0

      rl_ct_vr_df(ix_av_eta_dse) =       1d0
      rl_ct_vr_df(ix_av_ftran_sol_dse) =
     &     rl_ct_vr_df(ix_tl_fwd_tran_dse_rhs)
      rl_ct_vr_df(ix_av_btran_sol_dse) =
     &     rl_ct_vr_df(ix_tl_fwd_tran_dse_rhs)
      rl_ct_vr_df(ix_av_tbu_r_dse) =
     &     rl_ct_vr_df(ix_tl_fwd_tran_dse_rhs)
      rl_ct_vr_df(ix_tl_dvx_wt) =        3d0
      rl_ct_vr_df(ix_nw_dvx_fwk_fq) =    1d-2

      rl_ct_vr_df(ix_inv_pv_tl) =     1d-6
      rl_ct_vr_df(ix_inv_wr_eta_tl) = 1d-9
      rl_ct_vr_df(ix_inv_c_rlv_tl) =  1d-2
      rl_ct_vr_df(ix_inv_unit_v_tl) = 1d-6
      rl_ct_vr_df(ix_inv_l_pv_rsdu_er_tl) = 1d-6

      rl_ct_vr_df(ix_usr_cz_r_pv_tl) = 0d0
      rl_ct_vr_df(ix_wk_cz_r_pv_tl) =  rl_ct_vr_df(ix_usr_cz_r_pv_tl)
      rl_ct_vr_df(ix_cz_r_pv_tl) =     rl_ct_vr_df(ix_usr_cz_r_pv_tl)

      rl_ct_vr_df(ix_tl_bc_pr_act_it_rfn) = rl_ct_vr_df(ix_tl_pr_ifs)
      rl_ct_vr_df(ix_tl_reset_pi_it_rfn) =  rl_ct_vr_df(ix_tl_du_ifs)
      rl_ct_vr_df(ix_tl_pv_c_it_rfn) =      rl_ct_vr_df(ix_tl_pr_ifs)
      rl_ct_vr_df(ix_tl_full_pi_it_rfn) =   rl_ct_vr_df(ix_tl_du_ifs)
      rl_ct_vr_df(ix_tl_u_pi_it_rfn) =      rl_ct_vr_df(ix_tl_du_ifs)
      rl_ct_vr_df(ix_bc_pr_act_it_rfn_tran_ze) = zero
      rl_ct_vr_df(ix_reset_pi_it_rfn_tran_ze) =  zero
      rl_ct_vr_df(ix_pv_c_it_rfn_tran_ze) =      zero
      rl_ct_vr_df(ix_full_pi_it_rfn_tran_ze) =   zero
      rl_ct_vr_df(ix_u_pi_it_rfn_tran_ze) =      zero

      call ems_cp_lg_a(n_ems_rl_ct_vr, .false., rl_ct_vr_lg, 0)
      rl_ct_vr_lg(ix_usr_tl_pr_ifs) = .true.
      rl_ct_vr_lg(ix_usr_tl_du_ifs) = .true.
      rl_ct_vr_lg(ix_mx_mn) = .true.
      rl_ct_vr_lg(ix_pr_wt) = .true.
      rl_ct_vr_lg(ix_bb_ub) = .true.
      rl_ct_vr_lg(ix_dga_scl) = .true.
      rl_ct_vr_lg(ix_mps_ze) = .true.
      rl_ct_vr_lg(ix_i_ifs_wt) = .true.
      rl_ct_vr_lg(ix_tl_i_fs_ob_fn_v_imp) = .true.
      rl_ct_vr_lg(ix_i_fs_ob_fn_v_tgt) = .true.
      rl_ct_vr_lg(ix_tl_i_ifs) = .true.

      rl_ct_vr_lg(ix_tl_iz_bs_tran_er) =  .true.

      rl_ct_vr_lg(ix_tl_bc_pr_act_it_rfn) = .true.
      rl_ct_vr_lg(ix_tl_reset_pi_it_rfn) = .true.
      rl_ct_vr_lg(ix_tl_pv_c_it_rfn) = .true.
      rl_ct_vr_lg(ix_tl_full_pi_it_rfn) = .true.
      rl_ct_vr_lg(ix_tl_u_pi_it_rfn) = .true.
      rl_ct_vr_lg(ix_bc_pr_act_it_rfn_tran_ze) = .true.
      rl_ct_vr_lg(ix_reset_pi_it_rfn_tran_ze) = .true.
      rl_ct_vr_lg(ix_pv_c_it_rfn_tran_ze) = .true.
      rl_ct_vr_lg(ix_full_pi_it_rfn_tran_ze) = .true.
      rl_ct_vr_lg(ix_u_pi_it_rfn_tran_ze) = .true.

CM      IF (emsol_dev .EQ. 1) THEN
C?      rl_ct_vr_lg(ix_wk_tl_pr_ifs) = .true.
C?      rl_ct_vr_lg(ix_wk_tl_du_ifs) = .true.
C?      rl_ct_vr_lg(ix_wk_xp_tau) = .true.
C?      rl_ct_vr_lg(ix_tl_bwd_tran_dse_rhs) = .true.
C?      rl_ct_vr_lg(ix_tl_fwd_tran_dse_rhs) = .true.
C?      rl_ct_vr_lg(ix_tl_dse_tbu_r) = .true.
C?
C?      rl_ct_vr_lg(ix_usr_tl_crsh_abs_pv_v) = .true.
C?      rl_ct_vr_lg(ix_usr_tl_crsh_rlv_pv_v) = .true.
C?      rl_ct_vr_lg(ix_tl_mx_iz_pr_act) =  .true.
C?
C?      rl_ct_vr_lg(ix_tl_du_act_er) = .true.
C?      rl_ct_vr_lg(ix_tl_pv_er) = .true.
C?      rl_ct_vr_lg(ix_tl_ed_wt_er) = .true.
C?      rl_ct_vr_lg(ix_tl_tran_er) = .true.
C?      rl_ct_vr_lg(ix_tl_cz_r_growth) = .true.
C?
C?      rl_ct_vr_lg(ix_bwd_tran_ze) = .true.
C?      rl_ct_vr_lg(ix_pc_ze) = .true.
C?      rl_ct_vr_lg(ix_fwd_tran_ze) = .true.
C?      rl_ct_vr_lg(ix_u_ed_wt_ze) = .true.
C?      rl_ct_vr_lg(ix_pk_pv_c_ze) = .true.
C?      rl_ct_vr_lg(ix_scl_v_ze) = .true.
C?
C?      rl_ct_vr_lg(ix_tl_dvx_wt) =  .true.
C?      rl_ct_vr_lg(ix_nw_dvx_fwk_fq) = .true.
C?
C?      rl_ct_vr_lg(ix_inv_pv_tl) =  .true.
C?      rl_ct_vr_lg(ix_inv_wr_eta_tl) =  .true.
C?      rl_ct_vr_lg(ix_inv_c_rlv_tl) =  .true.
C?      rl_ct_vr_lg(ix_inv_unit_v_tl) =  .true.
C?      rl_ct_vr_lg(ix_inv_l_pv_rsdu_er_tl) =  .true.
C?
C?      rl_ct_vr_lg(ix_usr_cz_r_pv_tl) =  .true.
C?      rl_ct_vr_lg(ix_wk_cz_r_pv_tl) =  .true.
C?
CM      ENDIF
      return
      end

C->>> -----------------------------------------> ems_iz_osl_ct_vr_nm <<<
c     Set up an array of OSL names for integer and real control
c     variables.
c
      subroutine ems_iz_osl_ct_vr_nm
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'ICTVR.INC'
      include 'RLCTVR.INC'
c
c     Set the integer control variable names
c
      osl_i_ct_vr_nm( 1) = 'ILOGFREQ'
      osl_i_ct_vr_nm( 2) = 'IPRINTUNIT'
      osl_i_ct_vr_nm( 3) = 'IMAXFACTOR'
      osl_i_ct_vr_nm( 4) = 'IITERNUM'
      osl_i_ct_vr_nm( 5) = 'IMAXITER'
      osl_i_ct_vr_nm( 6) = 'ILOGLEVEL'
      osl_i_ct_vr_nm( 7) = 'IONEOBJ'
      osl_i_ct_vr_nm( 8) = 'IQPARNUMITER'
      osl_i_ct_vr_nm( 9) = 'IMAXROWS'
      osl_i_ct_vr_nm(10) = 'IMAXCOLS'
      osl_i_ct_vr_nm(11) = 'INUMCHAR'
      osl_i_ct_vr_nm(12) = 'ISTOPMASK'
      osl_i_ct_vr_nm(13) = 'IMAXITERB'
      osl_i_ct_vr_nm(14) = 'IADJACTYPE'
      osl_i_ct_vr_nm(15) = 'IFORMNTYPE'
      osl_i_ct_vr_nm(16) = 'IDENSECOL'
      osl_i_ct_vr_nm(17) = 'IDEVEXMODE'
      osl_i_ct_vr_nm(18) = 'INULLCHECK'
      osl_i_ct_vr_nm(19) = 'IDROPROWCT'
      osl_i_ct_vr_nm(20) = 'IITERUFREQ'
      osl_i_ct_vr_nm(21) = 'IPOSSBASIS'
      osl_i_ct_vr_nm(22) = 'IMAXPROJNS'
      osl_i_ct_vr_nm(23) = 'INUMELS'
      osl_i_ct_vr_nm(24) = 'INUMBLOCKS'
      osl_i_ct_vr_nm(25) = 'IMSGPOS'
      osl_i_ct_vr_nm(26) = 'IPAGELINES'
      osl_i_ct_vr_nm(27) = 'INUMROWS'
      osl_i_ct_vr_nm(28) = 'INUMCOLS'
      osl_i_ct_vr_nm(29) = 'INUMPINF'
      osl_i_ct_vr_nm(30) = 'INUMDINF'
      osl_i_ct_vr_nm(31) = 'IMODELMASK'
      osl_i_ct_vr_nm(32) = 'IPRTINFOMASK'
      osl_i_ct_vr_nm(33) = 'ISOLMASK'
      osl_i_ct_vr_nm(33) = 'IPRTMTRXMASK'
      osl_i_ct_vr_nm(34) = 'IEXTRABLK'
      osl_i_ct_vr_nm(35) = 'IQPARMAXITER'
      osl_i_ct_vr_nm(36) = 'ILINELEN'
      osl_i_ct_vr_nm(37) = 'INUMNODES'
      osl_i_ct_vr_nm(38) = 'IINTMASK'
      osl_i_ct_vr_nm(39) = 'IFASTITS'
      osl_i_ct_vr_nm(40) = 'IMAXNODES'
      osl_i_ct_vr_nm(41) = 'IMAXSOLS'
      osl_i_ct_vr_nm(42) = 'INUMSOLS'
      osl_i_ct_vr_nm(43) = 'INUMINTS'
      osl_i_ct_vr_nm(44) = 'INUMSETS'
      osl_i_ct_vr_nm(45) = 'INUMUNSAT'
      osl_i_ct_vr_nm(46) = 'IVECTOR'
      osl_i_ct_vr_nm(47) = 'IPROBSTAT'
      osl_i_ct_vr_nm(48) = 'IMAJORITS'
      osl_i_ct_vr_nm(49) = 'IPRINTSENS'
      osl_i_ct_vr_nm(50) = 'IMAXINTS'
      osl_i_ct_vr_nm(51) = 'IMAXSETS'
      osl_i_ct_vr_nm(52) = 'ISTRATEGY'
      osl_i_ct_vr_nm(53) = 'IMAXINTINFO'
      osl_i_ct_vr_nm(54) = 'ITHRESHOLD'
      osl_i_ct_vr_nm(55) = 'IHEURPASS'
      osl_i_ct_vr_nm(56) = 'ISUPERTOL'
      osl_i_ct_vr_nm(57) = 'IROWORD'
      osl_i_ct_vr_nm(58) = 'IORDUNIT'
      osl_i_ct_vr_nm(59) = 'IMIPLENGTH'
      osl_i_ct_vr_nm(60) = 'IITERBNUM'
      osl_i_ct_vr_nm(61) = 'IPRICETYPE'
      osl_i_ct_vr_nm(62) = 'INUMCPU'
      osl_i_ct_vr_nm(63) = 'IPROBSTAT2'
      osl_i_ct_vr_nm(64) = 'ISMDLTYPEMASK'
      osl_i_ct_vr_nm(65) = 'ITOTALCPU'
      osl_i_ct_vr_nm(66) = 'IWHICHCPU'
      osl_i_ct_vr_nm(67) = 'ISTRIPES'
      osl_i_ct_vr_nm(69) = 'INUMELQ'
      osl_i_ct_vr_nm(70) = 'IORDTHRSH'
      osl_i_ct_vr_nm(71) = 'ILPDCFLAG'
c
c     Set the real control variable names
c
      osl_rl_ct_vr_nm( 1) = 'RTOLPINF'
      osl_rl_ct_vr_nm( 2) = 'RTOLDINF'
      osl_rl_ct_vr_nm( 3) = 'RMAXMIN'
      osl_rl_ct_vr_nm( 4) = 'RMUFACTOR'
      osl_rl_ct_vr_nm( 5) = 'RMULIMIT'
      osl_rl_ct_vr_nm( 6) = 'RRGFACTOR'
      osl_rl_ct_vr_nm( 7) = 'RRGLIMIT'
      osl_rl_ct_vr_nm( 8) = 'RFIXVAR1'
      osl_rl_ct_vr_nm( 9) = 'RFIXVAR2'
      osl_rl_ct_vr_nm(10) = 'RCHOLABSTOL'
      osl_rl_ct_vr_nm(11) = 'RCHOLRELTOL'
      osl_rl_ct_vr_nm(12) = 'RMULINFAC'
      osl_rl_ct_vr_nm(13) = 'RPROJTOL'
      osl_rl_ct_vr_nm(14) = 'RPWEIGHT'
      osl_rl_ct_vr_nm(15) = 'RCHANGEWEIGHT'
      osl_rl_ct_vr_nm(16) = 'RBBCUTOFF'
      osl_rl_ct_vr_nm(17) = 'RDWEIGHT'
      osl_rl_ct_vr_nm(18) = 'ROBJVALUE'
      osl_rl_ct_vr_nm(19) = 'RSUMPINF'
      osl_rl_ct_vr_nm(20) = 'RSUMDINF'
      osl_rl_ct_vr_nm(21) = 'RTOLMPS'
      osl_rl_ct_vr_nm(22) = 'RDEGSCALE'
      osl_rl_ct_vr_nm(23) = 'RBESTSOL'
      osl_rl_ct_vr_nm(24) = 'RIWEIGHT'
      osl_rl_ct_vr_nm(25) = 'RIMPROVE'
      osl_rl_ct_vr_nm(26) = 'RTARGET'
      osl_rl_ct_vr_nm(27) = 'RTOLINT'
      osl_rl_ct_vr_nm(28) = 'RBESTPOSS'
      osl_rl_ct_vr_nm(29) = 'RBESTEST'
      osl_rl_ct_vr_nm(30) = 'RSTEPMULT'
      osl_rl_ct_vr_nm(31) = 'RMUINIT'
      osl_rl_ct_vr_nm(32) = 'RDENSETHR'
      osl_rl_ct_vr_nm(33) = 'ROBJWEIGHT'
      osl_rl_ct_vr_nm(34) = 'RLAMBDAVAL'
      osl_rl_ct_vr_nm(35) = 'RDCCUTOFF'
      osl_rl_ct_vr_nm(36) = 'RDOBJVAL'
      osl_rl_ct_vr_nm(37) = 'RSLAMBDA'
      osl_rl_ct_vr_nm(38) = 'RSLAMBDALIM'
      osl_rl_ct_vr_nm(39) = 'RSLAMBDADELTA'
      osl_rl_ct_vr_nm(40) = 'RTHRESHOLD'
      osl_rl_ct_vr_nm(41) = 'RPDGAPTOL'
      osl_rl_ct_vr_nm(42) = 'RPDSTEPMULT'
      osl_rl_ct_vr_nm(43) = 'RPERTDIAG'
      osl_rl_ct_vr_nm(44) = 'RNETSAMP'
      osl_rl_ct_vr_nm(45) = 'RPRINTCPU'
      return
      end
C->>> -----------------------------------------> ems_iz_ems_ct_vr_nm <<<
c     Set up an array of EMSOL names for integer and real control
c     variables.
c
      subroutine ems_iz_ems_ct_vr_nm
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'ICTVR.INC'
      include 'RLCTVR.INC'
c
c     Set the integer control variable names
c
      ems_i_ct_vr_nm(ix_msg_fq) = 'msg_fq'
      ems_i_ct_vr_nm(ix_ems_wr_cn) = 'ems_wr_cn'
      ems_i_ct_vr_nm(ix_usr_mx_n_u) = 'usr_mx_n_u'
      ems_i_ct_vr_nm(ix_n_si_it) = 'n_si_it'
      ems_i_ct_vr_nm(ix_mx_n_si_it) = 'mx_n_si_it'
      ems_i_ct_vr_nm(ix_rsmi_msg_msk) = 'rsmi_msg_msk'
      ems_i_ct_vr_nm(ix_mx_n_r) = 'mx_n_r'
      ems_i_ct_vr_nm(ix_mx_n_c) = 'mx_n_c'
      ems_i_ct_vr_nm(ix_ml_nm_n_ch) = 'ml_nm_n_ch'
      ems_i_ct_vr_nm(ix_dvx_mode) = 'dvx_mode'
      ems_i_ct_vr_nm(ix_it_usr_xit_fq) = 'it_usr_xit_fq'
      ems_i_ct_vr_nm(ix_n_el_in_aux_blk) = 'n_el_in_aux_blk'
      ems_i_ct_vr_nm(ix_n_aux_blk) = 'n_aux_blk'
      ems_i_ct_vr_nm(ix_n_pg_li) = 'n_pg_li'
      ems_i_ct_vr_nm(ix_n_r) = 'n_r'
      ems_i_ct_vr_nm(ix_n_c) = 'n_c'
      ems_i_ct_vr_nm(ix_n_pr_ifs) = 'n_pr_ifs'
      ems_i_ct_vr_nm(ix_n_du_ifs) = 'n_du_ifs'
      ems_i_ct_vr_nm(ix_wr_ml_da_msk) = 'wr_ml_da_msk'
      ems_i_ct_vr_nm(ix_wr_ml_st_msk) = 'wr_ml_st_msk'
      ems_i_ct_vr_nm(ix_mx_n_el_in_aux_blk) = 'mx_n_el_in_aux_blk'
      ems_i_ct_vr_nm(ix_wr_li_ln) = 'wr_li_ln'
      ems_i_ct_vr_nm(ix_mx_nd_n) = 'mx_nd_n'
      ems_i_ct_vr_nm(ix_i_lp_msg_msk) = 'i_lp_msg_msk'
      ems_i_ct_vr_nm(ix_mx_n_nd) = 'mx_n_nd'
      ems_i_ct_vr_nm(ix_mx_n_i_fs_sol) = 'mx_n_i_fs_sol'
      ems_i_ct_vr_nm(ix_n_i_fs_sol) = 'n_i_fs_sol'
      ems_i_ct_vr_nm(ix_n_i_vr) = 'n_i_vr'
      ems_i_ct_vr_nm(ix_n_sos_se) = 'n_sos_se'
      ems_i_ct_vr_nm(ix_n_i_vr_no_i) = 'n_i_vr_no_i'
      ems_i_ct_vr_nm(ix_prob_st) = 'prob_st'
      ems_i_ct_vr_nm(ix_wr_rg_da_msk) = 'wr_rg_da_msk'
      ems_i_ct_vr_nm(ix_mx_n_i_vr) = 'mx_n_i_vr'
      ems_i_ct_vr_nm(ix_mx_n_sos_se) = 'mx_n_sos_se'
      ems_i_ct_vr_nm(ix_i_alg) = 'i_alg'
      ems_i_ct_vr_nm(ix_mx_n_i_da_en) = 'mx_n_i_da_en'
      ems_i_ct_vr_nm(ix_crsh_msk) = 'crsh_msk'
      ems_i_ct_vr_nm(ix_crsh_tl_dl_pri) = 'crsh_tl_dl_pri'
      ems_i_ct_vr_nm(ix_crsh_fn_dn) = 'crsh_fn_dn'
      ems_i_ct_vr_nm(ix_ems_er_cn) = 'ems_er_cn'
      ems_i_ct_vr_nm(ix_ts_parsmi) = 'ts_parsmi'
      ems_i_ct_vr_nm(ix_ems_fi_wr_msk) = 'ems_fi_wr_msk'
      ems_i_ct_vr_nm(ix_scl_mode) = 'scl_mode'
      ems_i_ct_vr_nm(ix_scl_bd) = 'scl_bd'
      ems_i_ct_vr_nm(ix_mx_n_a_el) = 'mx_n_a_el'
      ems_i_ct_vr_nm(ix_n_a_el) = 'n_a_el'
      ems_i_ct_vr_nm(ix_ml_blk_st_msk) = 'ml_blk_st_msk'
      ems_i_ct_vr_nm(ix_ml_da_st_msk) = 'ml_da_st_msk'
      ems_i_ct_vr_nm(ix_ml_da_no_cg_msk) = 'ml_da_no_cg_msk'
      ems_i_ct_vr_nm(ix_du_sol_mode) = 'du_sol_mode'
      ems_i_ct_vr_nm(ix_ml_nm_n_rl) = 'ml_nm_n_rl'
      ems_i_ct_vr_nm(ix_sslv_en_msk) = 'sslv_en_msk'
      ems_i_ct_vr_nm(ix_sslv_xit_msk) = 'sslv_xit_msk'
      ems_i_ct_vr_nm(ix_n_bs) = 'n_bs'
      ems_i_ct_vr_nm(ix_n_reset) = 'n_reset'
      ems_i_ct_vr_nm(ix_rq_reset) = 'rq_reset'
      ems_i_ct_vr_nm(ix_u_pc) = 'u_pc'
      ems_i_ct_vr_nm(ix_r_pc) = 'r_pc'
      ems_i_ct_vr_nm(ix_eta_fi_mode_msk) = 'eta_fi_mode_msk'
      ems_i_ct_vr_nm(ix_eta_fi_da_st_msk) = 'eta_fi_da_st_msk'
      ems_i_ct_vr_nm(ix_mn_n_dvx_it) = 'mn_n_dvx_it'
      ems_i_ct_vr_nm(ix_pc_alg) = 'pc_alg'
      ems_i_ct_vr_nm(ix_n_dvx_it) = 'n_dvx_it'
      ems_i_ct_vr_nm(ix_n_dvx_fwk) = 'n_dvx_fwk'
      ems_i_ct_vr_nm(ix_r_du_act_sgn) = 'r_du_act_sgn'
      ems_i_ct_vr_nm(ix_inv_alg_msk) = 'inv_alg_msk'
      ems_i_ct_vr_nm(ix_inv_msg_msk) = 'inv_msg_msk'
      ems_i_ct_vr_nm(ix_inv_log_msk) = 'inv_log_msk'
      ems_i_ct_vr_nm(ix_inv_mx_n_eta) = 'inv_mx_n_eta'
      ems_i_ct_vr_nm(ix_inv_mx_n_eta_el) = 'inv_mx_n_eta_el'
      ems_i_ct_vr_nm(ix_n_inv_sing) = 'n_inv_sing'
      ems_i_ct_vr_nm(ix_u_bs_cg) = 'u_bs_cg'
      ems_i_ct_vr_nm(ix_rq_inv) = 'rq_inv'
      ems_i_ct_vr_nm(ix_sto_ftran_ix_mode) = 'sto_ftran_ix_mode'
      ems_i_ct_vr_nm(ix_sto_btran_ix_mode) = 'sto_btran_ix_mode'
      ems_i_ct_vr_nm(ix_sto_tbu_r_ix_mode) = 'sto_tbu_r_ix_mode'
      ems_i_ct_vr_nm(ix_tbu_r_loop_mode) = 'tbu_r_loop_mode'
      ems_i_ct_vr_nm(ix_sto_ftran_ix) = 'sto_ftran_ix'
      ems_i_ct_vr_nm(ix_sto_btran_ix) = 'sto_btran_ix'
      ems_i_ct_vr_nm(ix_sto_tbu_r_ix) = 'sto_tbu_r_ix'
      ems_i_ct_vr_nm(ix_bwd_tran_dse_rhs_n_r) = 'bwd_tran_dse_rhs_n_r'
      ems_i_ct_vr_nm(ix_fwd_tran_dse_rhs_n_r) = 'fwd_tran_dse_rhs_n_r'
      ems_i_ct_vr_nm(ix_dse_tbu_r_n_c) = 'dse_tbu_r_n_c'
      ems_i_ct_vr_nm(ix_usr_cz_c_msk) = 'usr_cz_c_msk'
      ems_i_ct_vr_nm(ix_usr_cz_r_msk) = 'usr_cz_r_msk'
      ems_i_ct_vr_nm(ix_l1_cz_r_mx_n_cdd) = 'l1_cz_r_mx_n_cdd'
      ems_i_ct_vr_nm(ix_n_pc_vr) = 'n_pc_vr'
      ems_i_ct_vr_nm(ix_n_pc_el) = 'n_pc_el'
      ems_i_ct_vr_nm(ix_n_non_bc_pr_ifs) = 'n_non_bc_pr_ifs'
      ems_i_ct_vr_nm(ix_rq_re_pc) = 'rq_re_pc'
      ems_i_ct_vr_nm(ix_cz_c_msk) = 'cz_c_msk'
      ems_i_ct_vr_nm(ix_cz_r_msk) = 'cz_r_msk'
      ems_i_ct_vr_nm(ix_u_bs_mode) = 'u_bs_mode'
      ems_i_ct_vr_nm(ix_sto_u_eta_v_mode) = 'sto_u_eta_v_mode'
      ems_i_ct_vr_nm(ix_u_bs) = 'u_bs'
      ems_i_ct_vr_nm(ix_sto_u_eta_v) = 'sto_u_eta_v'
      ems_i_ct_vr_nm(ix_sto_u_eta_se_ty) = 'sto_u_eta_se_ty'
      ems_i_ct_vr_nm(ix_mx_n_u) = 'mx_n_u'
      ems_i_ct_vr_nm(ix_n_u) = 'n_u'
      ems_i_ct_vr_nm(ix_lo_sto_mode) = 'lo_sto_mode'
      ems_i_ct_vr_nm(ix_asm_msk) = 'asm_msk'
      ems_i_ct_vr_nm(ix_ck_msk) = 'ck_msk'
      ems_i_ct_vr_nm(ix_it_o_sw_fm_wr_df) = 'it_o_sw_fm_wr_df'
      ems_i_ct_vr_nm(ix_it_o_sw_fm_ck_df) = 'it_o_sw_fm_ck_df'
      ems_i_ct_vr_nm(ix_wr_lp_da) = 'wr_lp_da'
      ems_i_ct_vr_nm(ix_wr_bs_fq) = 'wr_bs_fq'
      ems_i_ct_vr_nm(ix_g_tt_da) = 'g_tt_da'
      ems_i_ct_vr_nm(ix_g_it_tt_fq) = 'g_it_tt_fq'
      ems_i_ct_vr_nm(ix_g_it_tt_ivl) = 'g_it_tt_ivl'
      ems_i_ct_vr_nm(ix_mx_n_bc_pr_act_rfn_it) = 'mx_n_bc_pr_act_rfn_it'
      ems_i_ct_vr_nm(ix_mx_n_reset_pi_rfn_it) = 'mx_n_reset_pi_rfn_it'
      ems_i_ct_vr_nm(ix_mx_n_pv_c_rfn_it) = 'mx_n_pv_c_rfn_it'
      ems_i_ct_vr_nm(ix_mx_n_full_pi_rfn_it) = 'mx_n_full_pi_rfn_it'
      ems_i_ct_vr_nm(ix_mx_n_u_pi_rfn_it) = 'mx_n_u_pi_rfn_it'
      ems_i_ct_vr_nm(ix_eta_fi_n_grp) = 'eta_fi_n_grp'
      ems_i_ct_vr_nm(ix_eta_fi_n_inv_grp) = 'eta_fi_n_inv_grp'
      ems_i_ct_vr_nm(ix_eta_fi_n_se) = 'eta_fi_n_se'
      ems_i_ct_vr_nm(ix_eta_fi_n_inv_se) = 'eta_fi_n_inv_se'
      ems_i_ct_vr_nm(ix_eta_fi_n_eta) = 'eta_fi_n_eta'
      ems_i_ct_vr_nm(ix_eta_fi_n_inv_eta) = 'eta_fi_n_inv_eta'
      ems_i_ct_vr_nm(ix_eta_fi_n_v) = 'eta_fi_n_v'
      ems_i_ct_vr_nm(ix_eta_fi_n_inv_v) = 'eta_fi_n_inv_v'
      ems_i_ct_vr_nm(ix_eta_fi_n_ix) = 'eta_fi_n_ix'
      ems_i_ct_vr_nm(ix_eta_fi_n_inv_ix) = 'eta_fi_n_inv_ix'
      ems_i_ct_vr_nm(ix_mx_n_pwl_vr) = 'mx_n_pwl_vr'
      ems_i_ct_vr_nm(ix_mx_n_pwl_vr_da_en) = 'mx_n_pwl_vr_da_en'
      ems_i_ct_vr_nm(ix_repl_non_std_vr_mode_msk) =
     &     'repl_non_std_vr_mode_msk'
      ems_i_ct_vr_nm(ix_repl_non_std_vr_st_msk) =
     &     'repl_non_std_vr_st_msk'
      ems_i_ct_vr_nm(ix_n_bp_vr) = 'n_bp_vr'
      ems_i_ct_vr_nm(ix_n_pwl_vr) = 'n_pwl_vr'
      ems_i_ct_vr_nm(ix_n_pwl_vr_da_en) = 'n_pwl_vr_da_en'
      ems_i_ct_vr_nm(ix_cu_is_n) = 'cu_is_n'
      ems_i_ct_vr_nm(ix_cu_ml_n) = 'cu_ml_n'
      ems_i_ct_vr_nm(ix_ml_blk_mv_k) = 'ml_blk_mv_k'
      ems_i_ct_vr_nm(ix_ds_n_en_m1) = 'ds_n_en_m1'
      ems_i_ct_vr_nm(ix_is_n_en_m1) = 'is_n_en_m1'
      ems_i_ct_vr_nm(ix_ns_n_en_m1) = 'ns_n_en_m1'
      ems_i_ct_vr_nm(ix_sv_ml_ct_vr) = 'sv_ml_ct_vr'
      ems_i_ct_vr_nm(ix_n_lo_c_eta) = 'n_lo_c_eta'
      ems_i_ct_vr_nm(ix_n_lo_r_eta) = 'n_lo_r_eta'
      ems_i_ct_vr_nm(ix_n_lo_r_eta_el) = 'n_lo_r_eta_el'
      ems_i_ct_vr_nm(ix_n_up_r_eta) = 'n_up_r_eta'
      ems_i_ct_vr_nm(ix_n_up_r_eta_el) = 'n_up_r_eta_el'
      ems_i_ct_vr_nm(ix_n_sus_eta_srch_op_cf) = 'n_sus_eta_srch_op_cf'
      ems_i_ct_vr_nm(ix_sus_eta_srch_al_t_ls_ix_lm) =
     &     'sus_eta_srch_al_t_ls_ix_lm'
      ems_i_ct_vr_nm(ix_sus_eta_srch_ls_t_buk_ix_lm) =
     &     'sus_eta_srch_ls_t_buk_ix_lm'
      ems_i_ct_vr_nm(ix_sus_eta_srch_dbl_n_buk_ix_lm) =
     &     'sus_eta_srch_dbl_n_buk_ix_lm'
c
c     Set the real control variable names
c
      ems_rl_ct_vr_nm(ix_usr_tl_pr_ifs) = 'usr_tl_pr_ifs'
      ems_rl_ct_vr_nm(ix_usr_tl_du_ifs) = 'usr_tl_du_ifs'
      ems_rl_ct_vr_nm(ix_mx_mn) = 'mx_mn'
      ems_rl_ct_vr_nm(ix_pr_wt) = 'pr_wt'
      ems_rl_ct_vr_nm(ix_bb_ub) = 'bb_ub'
      ems_rl_ct_vr_nm(ix_ob_fn_v) = 'ob_fn_v'
      ems_rl_ct_vr_nm(ix_su_pr_ifs) = 'su_pr_ifs'
      ems_rl_ct_vr_nm(ix_su_du_ifs) = 'su_du_ifs'
      ems_rl_ct_vr_nm(ix_mps_ze) = 'mps_ze'
      ems_rl_ct_vr_nm(ix_dga_scl) = 'dga_scl'
      ems_rl_ct_vr_nm(ix_bst_i_fs_ob_fn_v) = 'bst_i_fs_ob_fn_v'
      ems_rl_ct_vr_nm(ix_i_ifs_wt) = 'i_ifs_wt'
      ems_rl_ct_vr_nm(ix_tl_i_fs_ob_fn_v_imp) = 'tl_i_fs_ob_fn_v_imp'
      ems_rl_ct_vr_nm(ix_i_fs_ob_fn_v_tgt) = 'i_fs_ob_fn_v_tgt'
      ems_rl_ct_vr_nm(ix_tl_i_ifs) = 'tl_i_ifs'
      ems_rl_ct_vr_nm(ix_bst_ob_fn_v) = 'bst_ob_fn_v'
      ems_rl_ct_vr_nm(ix_bst_est_ob_fn_v) = 'bst_est_ob_fn_v'
      ems_rl_ct_vr_nm(ix_wk_tl_pr_ifs) = 'wk_tl_pr_ifs'
      ems_rl_ct_vr_nm(ix_wk_tl_du_ifs) = 'wk_tl_du_ifs'
      ems_rl_ct_vr_nm(ix_wk_xp_tau) = 'wk_xp_tau'
      ems_rl_ct_vr_nm(ix_tl_pr_ifs) = 'tl_pr_ifs'
      ems_rl_ct_vr_nm(ix_tl_du_ifs) = 'tl_du_ifs'
      ems_rl_ct_vr_nm(ix_xp_tau) = 'xp_tau'
      ems_rl_ct_vr_nm(ix_tl_bwd_tran_dse_rhs) = 'tl_bwd_tran_dse_rhs'
      ems_rl_ct_vr_nm(ix_tl_fwd_tran_dse_rhs) = 'tl_fwd_tran_dse_rhs'
      ems_rl_ct_vr_nm(ix_sus_eta_srch_co_mu) = 'sus_eta_srch_co_mu'
      ems_rl_ct_vr_nm(ix_tl_dse_tbu_r) = 'tl_dse_tbu_r'
      ems_rl_ct_vr_nm(ix_usr_tl_crsh_abs_pv_v) = 'usr_tl_crsh_abs_pv_v'
      ems_rl_ct_vr_nm(ix_usr_tl_crsh_rlv_pv_v) = 'usr_tl_crsh_rlv_pv_v'
      ems_rl_ct_vr_nm(ix_tl_iz_bs_tran_er) = 'tl_iz_bs_tran_er'
      ems_rl_ct_vr_nm(ix_tl_mx_iz_pr_act) = 'tl_mx_iz_pr_act'
      ems_rl_ct_vr_nm(ix_tl_du_act_er) = 'tl_du_act_er'
      ems_rl_ct_vr_nm(ix_tl_pv_er) = 'tl_pv_er'
      ems_rl_ct_vr_nm(ix_tl_ed_wt_er) = 'tl_ed_wt_er'
      ems_rl_ct_vr_nm(ix_tl_tran_er) = 'tl_tran_er'
      ems_rl_ct_vr_nm(ix_tl_cz_r_growth) = 'tl_cz_r_growth'
      ems_rl_ct_vr_nm(ix_fac_long_prd) = 'fac_long_prd'
      ems_rl_ct_vr_nm(ix_fac_short_prd) = 'fac_short_prd'
      ems_rl_ct_vr_nm(ix_bwd_tran_ze) = 'bwd_tran_ze'
      ems_rl_ct_vr_nm(ix_pc_ze) = 'pc_ze'
      ems_rl_ct_vr_nm(ix_fwd_tran_ze) = 'fwd_tran_ze'
      ems_rl_ct_vr_nm(ix_u_ed_wt_ze) = 'u_ed_wt_ze'
      ems_rl_ct_vr_nm(ix_pk_pv_c_ze) = 'pk_pv_c_ze'
      ems_rl_ct_vr_nm(ix_av_eta_dse) = 'av_eta_dse'
      ems_rl_ct_vr_nm(ix_av_ftran_sol_dse) = 'av_ftran_sol_dse'
      ems_rl_ct_vr_nm(ix_av_btran_sol_dse) = 'av_btran_sol_dse'
      ems_rl_ct_vr_nm(ix_av_tbu_r_dse) = 'av_tbu_r_dse'
      ems_rl_ct_vr_nm(ix_scl_v_ze) = 'scl_v_ze'
      ems_rl_ct_vr_nm(ix_tl_dvx_wt) = 'tl_dvx_wt'
      ems_rl_ct_vr_nm(ix_nw_dvx_fwk_fq) = 'nw_dvx_fwk_fq'
      ems_rl_ct_vr_nm(ix_inv_pv_tl) = 'inv_pv_tl'
      ems_rl_ct_vr_nm(ix_inv_wr_eta_tl) = 'inv_wr_eta_tl'
      ems_rl_ct_vr_nm(ix_inv_c_rlv_tl) = 'inv_c_rlv_tl'
      ems_rl_ct_vr_nm(ix_inv_unit_v_tl) = 'inv_unit_v_tl'
      ems_rl_ct_vr_nm(ix_inv_l_pv_rsdu_er_tl) = 'inv_l_pv_rsdu_er_tl'
      ems_rl_ct_vr_nm(ix_usr_cz_r_pv_tl) = 'usr_cz_r_pv_tl'
      ems_rl_ct_vr_nm(ix_wk_cz_r_pv_tl) = 'wk_cz_r_pv_tl'
      ems_rl_ct_vr_nm(ix_cz_r_pv_tl) = 'cz_r_pv_tl'
      ems_rl_ct_vr_nm(ix_tl_bc_pr_act_it_rfn) = 'tl_bc_pr_act_it_rfn'
      ems_rl_ct_vr_nm(ix_tl_reset_pi_it_rfn) = 'tl_reset_pi_it_rfn'
      ems_rl_ct_vr_nm(ix_tl_pv_c_it_rfn) = 'tl_pv_c_it_rfn'
      ems_rl_ct_vr_nm(ix_tl_full_pi_it_rfn) = 'tl_full_pi_it_rfn'
      ems_rl_ct_vr_nm(ix_tl_u_pi_it_rfn) = 'tl_u_pi_it_rfn'
      ems_rl_ct_vr_nm(ix_bc_pr_act_it_rfn_tran_ze) =
     &     'bc_pr_act_it_rfn_tran_ze'
      ems_rl_ct_vr_nm(ix_reset_pi_it_rfn_tran_ze) =
     &     'reset_pi_it_rfn_tran_ze'
      ems_rl_ct_vr_nm(ix_pv_c_it_rfn_tran_ze) =
     &     'pv_c_it_rfn_tran_ze'
      ems_rl_ct_vr_nm(ix_full_pi_it_rfn_tran_ze) =
     &     'full_pi_it_rfn_tran_ze'
      ems_rl_ct_vr_nm(ix_u_pi_it_rfn_tran_ze) = 'u_pi_it_rfn_tran_ze'
      ems_rl_ct_vr_nm(ix_rl_ml_nm) = 'rl_ml_nm'
      ems_rl_ct_vr_nm(ix_ob_fn_cs) = 'ob_fn_cs'
      ems_rl_ct_vr_nm(ix_su_non_bc_pr_ifs) = 'su_non_bc_pr_ifs'
      return
      end

C->>> ---------------------------------------------> ems_wr_rl_ct_vr <<<
      subroutine ems_wr_rl_ct_vr
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'RLCTVR.INC'
      include 'EMSMSG.INC'
      integer ix
      do 10, ix = 1, n_ems_rl_ct_vr
         if (rl_ct_vr_lg(ix)) then
            if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9000)
     &           ix, rl_ct_vr_lb(ix), rl_ct_vr_ub(ix),
     &           rl_ct_vr_df(ix), rl_ct_vr(ix)
            call ems_msg_wr_li(info_msg_n)
         else
            if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9010)
     &           ix, rl_ct_vr_lb(ix), rl_ct_vr_ub(ix),
     &           rl_ct_vr_df(ix), rl_ct_vr(ix)
            call ems_msg_wr_li(info_msg_n)
         end if
         if (rl_ct_vr(ix) .lt. rl_ct_vr_lb(ix) .or.
     &        rl_ct_vr(ix) .gt. rl_ct_vr_ub(ix)) then
            if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9100)ix
            call ems_msg_wr_li(er_msg_n)
         end if
 10   continue
      return
 9000 format('rl_ct_vr ', i3, '     is settable: range [', g11.4, ', ',
     &     g11.4, ']; default ', g11.4, '; value ', g11.4)
 9010 format('rl_ct_vr ', i3, ' is not settable: range [', g11.4, ', ',
     &     g11.4, ']; default ', g11.4, '; value ', g11.4)
 9100 format('rl_ct_vr ', i3, ' violates its range ')
      end

C->>> ---------------------------------------------> ems_iz_ch_ct_vr <<<
c     Initialise the character control variables.
c
      subroutine ems_iz_ch_ct_vr
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'CHCTVR.INC'
      ch_ml_nm =  ch_un_nm
      ch_ob_nm =  ch_un_nm
      ch_rhs_nm = ch_un_nm
      ch_rg_nm =  ch_un_nm
      ch_bd_nm =  ch_un_nm
      ch_bs_nm =  ch_un_nm
      return
      end

C->>> ---------------------------------> ems_cp_ct_vr_t_or_fm_bs_blk <<<
c     Copies the control variables
c     to   the base block (if t=.true. ) or
c     from the base block (if t=.false.).
c
      subroutine ems_cp_ct_vr_t_or_fm_bs_blk(t, p_bs_blk, ds, is)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'EMSMMGR.INC'
      include 'EMSMEM.INC'
      include 'CHCTVR.INC'
      include 'ICTVR.INC'
      include 'RLCTVR.INC'
      include 'EMSMSG.INC'
      logical t
      double precision ds(0:*)
      integer is(0:*), p_bs_blk
      integer ems_mem_mgr_nx_i_n_en_p
      integer ln, og_p_bs_blk, p, ch_n

      og_p_bs_blk = p_bs_blk
      p_bs_blk = ems_mem_mgr_nx_i_n_en_p(p_bs_blk, rl_wo_z)
      p = p_bs_blk/rl_wo_z
      ln = n_ems_rl_ct_vr*rl_wo_z
      if (t) then
         call ems_cp_rl_a1(n_ems_rl_ct_vr, rl_ct_vr, ds(p))
      else
         call ems_cp_rl_a1(n_ems_rl_ct_vr, ds(p), rl_ct_vr)
      end if
      p_bs_blk = p_bs_blk + ln
c
      p_bs_blk = ems_mem_mgr_nx_i_n_en_p(p_bs_blk, ch_wo_z)
      p = p_bs_blk/i_wo_z
      ln = n_ch_ct_vr*ch_wo_z
      if (t) then
         do 10, ch_n = 1, n_ch_ct_vr
            call ems_ch8_t_i(ch_ct_vr(ch_n), is(p))
            p = p + ch_wo_z
 10      continue
      else
         do 20, ch_n = 1, n_ch_ct_vr
            call ems_i_t_ch8(is(p), ch_ct_vr(ch_n))
            p = p + ch_wo_z
 20      continue
      end if
      p_bs_blk = p_bs_blk + ln
c
      p_bs_blk = ems_mem_mgr_nx_i_n_en_p(p_bs_blk, i_wo_z)
      p = p_bs_blk/i_wo_z
      ln = n_ems_i_ct_vr*i_wo_z
      if (t) then
         call ems_cp_i_a1(n_ems_i_ct_vr, i_ct_vr, is(p))
      else
         call ems_cp_i_a1(n_ems_i_ct_vr, is(p), i_ct_vr)
      end if
      p_bs_blk = p_bs_blk + ln
c
      if (p_bs_blk-og_p_bs_blk .gt. ct_vr_n_wo) go to 8102
      if (.not. t) then
         ems_msg_wr_cn = ems_wr_cn
         ems_msg_er_cn = ems_er_cn
         ems_msg_wr_li_ln = wr_li_ln
         ems_msg_n_pg_li = n_pg_li
      endif
 7000 continue
      return
 8102 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9820)
     &     p_bs_blk-og_p_bs_blk, ct_vr_n_wo
      call ems_msg_wr_li(serious_msg_n)
CM      IF (emsol_deb .EQ. 1) THEN
C?      call ems_dump
CM      ENDIF
      go to 7000
 9820 format(i9, ' words required for ct_vr rather than ', i9)
      end

max_stars_repo_name: ERGO-Code/EMSOL

max_issues_repo_name: ERGO-Code/EMSOL

max_forks_repo_name: ERGO-Code/EMSOL


!! ALL MODULES

!! SUBROUTINES FOR WAVE FUNCTION PROPAGATION including core prop_loop subroutine
!! prop_loop calls cmf_prop_wfn which calls cmf_prop_avector and propspfs
!! for regular forward time propagation (improvedrelaxflag=0)

#include "Definitions.INC"

module prop_loop_sub_mod
contains

!! ***************************************************************************************** !!
!!                        EXACT (VMF) PROPAGATION ROUTINE
!! ***************************************************************************************** !!

  subroutine prop_wfn(tin, tout)
    use parameters
    use xxxmod
    use mpimod
    use orbdermod
    use configmod
    use basissubmod
    use mpisubmod
    use odexmod
    use getstuffmod
    use spfsubmod
    implicit none
    real*8,intent(in) :: tout, tin
    real*8 :: mytime,nullreal, nulldouble,gbsstepsize
    integer ::  iflag,zzz, nullint, idid, itime, jtime, time=0,  time2=0 , numiters=0,rkworkdim,&
         rkiworkdim,psilength
    real*8, allocatable :: rkwork(:)
    integer, allocatable :: rkiwork(:)
    DATATYPE,allocatable :: psivec(:)

    call myclock(itime)

    if (numfrozen.gt.0) then
       print *, "program frozen in all_derivs for VMF.";     stop
    endif
    mytime=tin

    zzz=2
#ifdef REALGO
    zzz=1
#endif

    psilength=totspfdim+tot_adim*mcscfnum
    allocate(psivec(psilength))
    psivec(tot_adim*mcscfnum+1:tot_adim*mcscfnum+totspfdim)=yyy%cmfspfs(:,0)
    if (tot_adim.gt.0) then
       psivec(1:tot_adim*mcscfnum)=RESHAPE(yyy%cmfavec(:,:,0),(/tot_adim*mcscfnum/))
    endif

    if (nprocs.gt.1.and.(parorbsplit.ne.0.or.par_consplit.ne.0)) then
       OFLWR "error, for parallel options can't do VMF."
       WRFL  "set parorbsplit=0 and par_consplit=0 in namelist &parinp"; CFLST
    endif

    if (intopt==0) then
       rkworkdim=(3+6*psilength)*zzz; rkiworkdim=100
       allocate(rkwork(rkworkdim),rkiwork(rkiworkdim))
       rkiwork(:)=0; rkwork=0
       rkwork(psilength*zzz+1)=0.00001d0
       iflag=1
       call rkf45(all_derivs, psilength*zzz, psivec(:), mytime, tout, &
            relerr, abserr, iflag, rkwork, rkiwork)
       if (iflag/=2) then
          OFLWR; WRFL "RK45F ERR   Iflag = ", iflag, abserr, rkiwork(1); CFLST
       endif
       numiters=rkiwork(1)
       deallocate(rkwork,rkiwork)
    else if (intopt.eq.1) then
       rkworkdim=20*psilength*zzz; rkiworkdim=20*psilength*zzz
       allocate(rkwork(rkworkdim),rkiwork(rkiworkdim))
       rkwork=0; rkiwork=0
       call odexq(psilength*zzz,gbs_derivs,mytime,psivec(:),tout,gbsstepsize,relerr,abserr,&
            0,dummysub,0,rkwork,rkworkdim,rkiwork,rkiworkdim,idid)
       if (.not.(idid.eq.1)) then
          OFLWR; WRFL "ERR ODEX", idid; CFLST
       endif
       numiters=rkiwork(17)
       deallocate(rkwork,rkiwork)
    else
       OFLWR "INTOPT NOT REC", intopt; CFLST
    endif
!!    OFLWR "       Number of evaluations: ", numiters;  CFL

    yyy%cmfspfs(:,0)=psivec(tot_adim*mcscfnum+1:tot_adim*mcscfnum+totspfdim)
    if (tot_adim.gt.0) then
       yyy%cmfavec(:,:,0)=RESHAPE(psivec(1:tot_adim*mcscfnum),(/tot_adim,mcscfnum/))
    endif
    deallocate(psivec)

    call myclock(jtime);  time=time+jtime-itime;  itime=jtime

    call basis_project(www,numr,yyy%cmfavec(:,:,0))

    call apply_spf_constraints(yyy%cmfspfs(:,0))

    call spf_orthogit(yyy%cmfspfs(:,0), nulldouble)

!! prevent drift.. need timings
    if (parorbsplit.ne.3) then
       call mympibcast(yyy%cmfspfs(:,0),1,totspfdim)
    endif
    if (par_consplit.eq.0) then
       call mympibcast(yyy%cmfavec(:,:,0),1,tot_adim*mcscfnum)
    endif

    call get_stuff(tout)

    call myclock(jtime);  time2=time2+jtime-itime

! rkiwork!
!  if ((myrank.eq.1).and.(notiming.eq.0)) then
!     open(853, file=timingdir(1:getlen(timingdir))//"/vmf_prop.time.dat", &
!             status="unknown", position="append")
!     write(853,'(A3,F12.3,100I15)') "T=",  tout, time/1000, time2/1000,rkiwork(1);
!     close(853)
!  endif

  contains
    subroutine dummysub()
    end subroutine dummysub

  end subroutine prop_wfn

  subroutine propspfs(inspfs,outspfs,tin, tout,inlinearflag,numiters)
    use parameters
    implicit none
    DATATYPE,intent(in) :: inspfs(spfsize,nspf)
    DATATYPE,intent(out) :: outspfs(spfsize,nspf)
    DATATYPE,allocatable :: tempspfs2(:,:)
    integer :: inlinearflag,numiters,jj,k
    real*8 :: tout, tin,timea,timeb

    outspfs(:,:)=inspfs(:,:)
    numiters=0

    allocate(tempspfs2(spfsize,nspf));   tempspfs2(:,:)=0

    do k=1,littlesteps
       timea=tin+(tout-tin)*(k-1)/littlesteps;     timeb=tin+(tout-tin)*k/littlesteps

       call propspfs0(outspfs,tempspfs2,timea,timeb,inlinearflag,jj)
       numiters=numiters+jj
       outspfs(:,:)=tempspfs2(:,:)
    enddo

    deallocate(tempspfs2)

  end subroutine propspfs

  subroutine propspfs0(inspfs,outspfs,tin, tout,inlinearflag,numiters)
    use parameters
    use linearmod
    use orbdermod
    use odexmod
    use expospfpropmod
    use verletmod
    use spfsubmod
    implicit none
    DATATYPE,intent(in) :: inspfs(spfsize,nspf)
    DATATYPE,intent(out) :: outspfs(spfsize,nspf)
    integer ::  numiters,iflag,zzz, nullint, idid, inlinearflag,rkworkdim,rkiworkdim
    real*8 :: tout, tin,nullreal, time1, time2, nulldouble,gbsstepsize
    real*8, allocatable :: rkwork(:)
    integer, allocatable :: rkiwork(:)

    time1=tin;     time2=tout; numiters=0

    zzz=2
#ifdef REALGO
    zzz=1
#endif

    effective_cmf_linearflag=inlinearflag     !!LINEARMOD

    outspfs(:,:)=inspfs(:,:)

    if (intopt==0) then
       rkworkdim=(3+6*totspfdim)*zzz;      rkiworkdim=100
       allocate(rkwork(rkworkdim),rkiwork(rkiworkdim))
       iflag=1    ; rkwork=0;  rkiwork=0
       call rkf45(spf_linear_derivs, totspfdim*zzz, outspfs, time1, time2, &
            relerr, abserr, iflag, rkwork, rkiwork)
       if (iflag/=2) then
          OFLWR "RK45F ERR   Iflag = ", iflag, abserr, relerr; CFLST
       endif
       numiters=rkiwork(1)
       deallocate(rkwork,rkiwork)
    else if (intopt==1) then
       gbsstepsize=1.d-6
       rkworkdim=20*totspfdim*zzz; rkiworkdim=20*totspfdim*zzz
       allocate(rkwork(rkworkdim),rkiwork(rkiworkdim))
       rkwork=0; rkiwork=0
       call odexq(totspfdim*zzz,gbs_linear_derivs,time1,outspfs,time2,gbsstepsize,relerr,abserr,&
            0,dummysub,0,rkwork,rkworkdim,rkiwork,rkiworkdim,idid)
       numiters=rkiwork(17)
       if (.not.(idid.eq.1)) then
          OFLWR "ERR ODEX", idid; CFLST
       endif
       deallocate(rkwork,rkiwork)
    else if (intopt.eq.3) then
       call expospfprop(time1,time2,outspfs,numiters)
    else if (intopt.eq.4) then
       call verlet(outspfs,time1,time2,numiters)
    else
       OFLWR "Intopt not recognized: ", intopt; CFLST
    endif  !! intopt

    call apply_spf_constraints(outspfs)

    call spf_orthogit(outspfs, nulldouble)

  contains
    subroutine dummysub()
    end subroutine dummysub

  end subroutine propspfs0

!! ***************************************************************************************** !!
!!                        MEAN FIELD (CMF) PROPAGATION ROUTINE
!! ***************************************************************************************** !!

  subroutine cmf_prop_wfn(tin, tout)
    use parameters
    use linearmod
    use xxxmod
    use configmod
    use mpimod
    use mpisubmod
    use derivativemod !! conpropspfs
    use quadavecmod
    use quadspfsmod
    use configstuffmod
    use repnatmod
    use denutilmod
    use getstuffmod
    use dfconsubmod
    use fockrepsubmod
    use meansubmod
    implicit none
    real*8,intent(in) :: tout, tin
    integer,save :: times(0:20)=0
    integer, save :: xxcount=0
    integer :: numaiters,numiters,linearflag,imc,printflag=1,getlen,qq,&
         myiostat,ii,itime,jtime
    real*8 :: time1, time2
    DATAECS :: myvalues(mcscfnum)

    firsttime=tin;   lasttime=tout

    numiters=0;   numaiters=0

    xxcount=xxcount+1

    call myclock(itime)

!!     MOVE SAVED DATA (MEAN FIELDS ETC.) ONE TIME SLOT BACK.

!! should clean this up into derived type

    yyy%invdenmat(:,:,1) = yyy%invdenmat(:,:,0)
    yyy%denmat(:,:,1) = yyy%denmat(:,:,0)
    yyy%reducedinvrsq(:,:,1)=yyy%reducedinvrsq(:,:,0)
    yyy%reducedinvr(:,:,1)=yyy%reducedinvr(:,:,0)
    yyy%reducedr(:,:,1)=yyy%reducedr(:,:,0)
    if (nonuc_checkflag/=1) then
       yyy%reducedproderiv(:,:,1)=yyy%reducedproderiv(:,:,0)
    endif
    yyy%reducedpot(:,:,:,1) = yyy%reducedpot(:,:,:,0)
    yyy%reducedpottally(:,:,:,:,1) = yyy%reducedpottally(:,:,:,:,0)
    if (drivingflag.ne.0.and.tot_adim.gt.0) then
       yyy%drivingavectorsxx(:,:,:,1)=yyy%drivingavectorsxx(:,:,:,0)
       yyy%drivingavectorsyy(:,:,:,1)=yyy%drivingavectorsyy(:,:,:,0)
       yyy%drivingavectorszz(:,:,:,1)=yyy%drivingavectorszz(:,:,:,0)
       yyy%drivingorbsxx(:,:,1)=yyy%drivingorbsxx(:,:,0)
       yyy%drivingorbsyy(:,:,1)=yyy%drivingorbsyy(:,:,0)
       yyy%drivingorbszz(:,:,1)=yyy%drivingorbszz(:,:,0)
    endif

    if (numfrozen.ne.0) then
       yyy%frozenexchinvr(:,:,1) = yyy%frozenexchinvr(:,:,0)
       yyy%frozenexchmat(:,:,1) = yyy%frozenexchmat(:,:,0)
    endif

    if (use_fockmatrix) then
       yyy%fockmatrix(:,:,1) = yyy%fockmatrix(:,:,0)
    endif

    call assign_cptr(yyy%cptr(1),yyy%cptr(0),DATAONE)

    if (sparseopt.ne.0) then
       call assign_sptr(yyysptr(1),yyysptr(0),DATAONE,www)
    endif
    if (use_dfwalktype) then
       if (sparseopt.ne.0) then
          if (shuffle_dfwalktype) then
             call assign_sptr(yyysfdptr(1),yyysfdptr(0),DATAONE,fdww)
          else
             call assign_sptr(yyysdfptr(1),yyysdfptr(0),DATAONE,dfww)
          endif
       endif
    endif
    yyy%cmfspfs(:,1)= yyy%cmfspfs(:,0)
    if (tot_adim.gt.0) then
       yyy%cmfavec(:,:,1)= yyy%cmfavec(:,:,0)
    endif

    call myclock(jtime);  times(6)=times(6)+jtime-itime

    if ( (threshflag.ne.0.and.improvedrelaxflag.le.0).and.&
         ( ( (improvedquadflag.eq.1.or.improvedquadflag.eq.3).and.tin.ge.aquadstarttime ).or.&
         ( (improvedquadflag.ge.2).and.tin.ge.quadstarttime ) ) ) then
       OFLWR " ** SETTING IMPROVEDRELAXFLAG=1 ** "; CFL
       improvedrelaxflag=max(1,(-1)*improvedrelaxflag)
    endif

    if (improvedrelaxflag.gt.0) then

       call myclock(itime)

       if (spf_flag.ne.0) then
          if (improvedquadflag.gt.1.and.tin.ge.quadstarttime) then
             call quadspfs(yyy%cmfspfs(:,0), qq)
             numiters=numiters+qq
          else
             time1=tin;        time2=tout
             call propspfs(yyy%cmfspfs(:,1), yyy%cmfspfs(:,0), time1,time2,0,qq)
             numiters=numiters+qq
          endif
          call myclock(jtime);     times(4)=times(4)+jtime-itime;     itime=jtime
!! prevent drift
          if (parorbsplit.ne.3) then
             call mympibcast(yyy%cmfspfs(:,0),1,totspfdim)
          endif
          call myclock(jtime);     times(9)=times(9)+jtime-itime;   itime=jtime
       endif

       if (improvednatflag.ne.0) then
          call replace_withnat(printflag)
       elseif (improvedfockflag.ne.0) then
          call replace_withfock(printflag)
       endif
       call myclock(jtime);     times(7)=times(7)+jtime-itime;   itime=jtime

       call all_matel()
       call myclock(jtime);     times(1)=times(1)+jtime-itime;   itime=jtime

       if (avector_flag.ne.0) then
          if ((improvedquadflag.eq.1.or.improvedquadflag.eq.3).and.tin.ge.aquadstarttime) then
             call quadavector(yyy%cmfavec(:,:,0),qq)
             numaiters=numaiters+qq
          else
             if (followflag.ne.0) then
                call myconfigeig(yyy%cptr(0),yyy%cmfavec(:,:,0),&
                     myvalues,mcscfnum,eigprintflag,2,0d0,max(0,improvedrelaxflag-1))
             else
                call myconfigeig(yyy%cptr(0),yyy%cmfavec(:,:,0),&
                     myvalues,mcscfnum,eigprintflag,1,0d0,max(0,improvedrelaxflag-1))
             endif
          endif
          call myclock(jtime);     times(5)=times(5)+jtime-itime;    itime=jtime
!! prevent drift
          if (par_consplit.eq.0) then
             call mympibcast(yyy%cmfavec(:,:,0),1,tot_adim*mcscfnum)
          endif
          call myclock(jtime);     times(9)=times(9)+jtime-itime;    itime=jtime
       endif

       call get_allden()
       call myclock(jtime);        times(2)=times(2)+jtime-itime;     itime=jtime

       if (use_fockmatrix) then
          call get_fockmatrix()
       endif
       call get_reducedpot()
       if (numfrozen.gt.0) then
          call get_frexchange()
       endif
       call myclock(jtime);        times(3)=times(3)+jtime-itime;       itime=jtime

       if (constraintflag.ne.0) then    !! probably need denmat right
          call get_constraint(tout);
          call myclock(jtime);        times(7)=times(7)+jtime-itime
       endif

    else  !! IMPROVEDRELAX

!! ******* CMF PREDICTOR ***** !!
!!            and              !!
!! ******* LMF CORRECTOR ***** !!

       linearflag=0

       call myclock(itime)

!! PRE-PROP

       if (prepropflag.ne.0) then
          if (constraintflag.ne.0) then
             time1=tin;        time2=tout
             call conpropspfs(yyy%cmfspfs(:,1),yyy%cmfspfs(:,0), time1,time2)
             call myclock(jtime);     times(4)=times(4)+jtime-itime;   itime=jtime
!! prevent drift
             if (parorbsplit.ne.3) then
                call mympibcast(yyy%cmfspfs(:,0),1,totspfdim)
                call myclock(jtime);     times(9)=times(9)+jtime-itime;   itime=jtime
             endif
             call all_matel()
             call myclock(jtime);  times(1)=times(1)+jtime-itime;   itime=jtime
          endif

          if(avector_flag.ne.0) then
             call cmf_prop_avector_and_stuff()

             call myclock(itime)
             call get_allden()
             call myclock(jtime);  times(2)=times(2)+jtime-itime;   itime=jtime
          endif

          call end_stuff()

          linearflag=1

       endif  !! prepropflag

!! MAIN STEPS

       if (step_flag.lt.0) then
          OFLWR "what? step_flag", step_flag; CFLST
       endif

       do ii=1,step_flag
          call propspfs_and_stuff()
          call cmf_prop_avector_and_stuff()
          call get_stuff0(tout,times)
          linearflag=1
       enddo

!! POST-PROP

       if (postpropflag.eq.1.and.spf_flag.ne.0) then

          call propspfs_and_stuff()

          call myclock(itime)
          call all_matel()
          call myclock(jtime);  times(1)=times(1)+jtime-itime

          call end_stuff()

       endif

       if (postpropflag.eq.2.and.avector_flag.ne.0) then

          call cmf_prop_avector_and_stuff()

          call myclock(itime)
          call get_allden()
          call myclock(jtime);  times(2)=times(2)+jtime-itime

          call end_stuff()

       endif  !! postpropflag

    endif  !!improvedrelax

!!!!! ******** !!!!!!!
!!!!! ******** !!!!!!!
!!!!! ******** !!!!!!!

    if (myrank.eq.1.and.notiming.le.1) then
       call output_denmat(1,tout)
    endif

    if ((myrank.eq.1).and.(notiming.le.1)) then
       if (xxcount==1) then
          open(853, file=timingdir(1:getlen(timingdir))//"/cmf_prop.time.dat", &
               status="unknown",iostat=myiostat)
          call checkiostat(myiostat," opening cmf_prop_time.dat")
          write(853,'(A15,100A11)',iostat=myiostat) &
               "Time", &      !!
               "matel", &     !! (1)
               "denmat", &    !! (2)
               "reduced",&    !! (3)
               "spfprop",&    !! (4)
               "aprop", &     !! (5)
               "advance", &   !! (6)
               "constrain", & !! (7)
               "driving",&    !! (8)
               "MPI",&        !! (9)
               "#SDERIVS",&   !!
               "#ADERIVS"     !!
          call checkiostat(myiostat," writing cmf_prop_time.dat")
          close(853)
       endif
       open(853, file=timingdir(1:getlen(timingdir))//"/cmf_prop.time.dat", &
            status="unknown", position="append",iostat=myiostat)
       call checkiostat(myiostat," opening cmf_prop_time.dat")
       write(853,'(A3,F12.3,T16, 100I11)',iostat=myiostat)  "T=", tout, &
            times(1:9)/1000, numiters,numaiters
       call checkiostat(myiostat," writing cmf_prop_time.dat")
       close(853)
    endif

  contains
    subroutine propspfs_and_stuff()
      implicit none
      if (spf_flag.ne.0) then
         call myclock(itime)
         time1=tin;        time2=tout
         call propspfs(yyy%cmfspfs(:,1),yyy%cmfspfs(:,0), time1,time2,linearflag,qq)
         numiters=numiters+qq
         call myclock(jtime);     times(4)=times(4)+jtime-itime;   itime=jtime
!! prevent drift
         if (parorbsplit.ne.3) then
            call mympibcast(yyy%cmfspfs(:,0),1,totspfdim)
            call myclock(jtime);     times(9)=times(9)+jtime-itime
         endif
      endif
    end subroutine propspfs_and_stuff

    subroutine cmf_prop_avector_and_stuff()
      implicit none
      if(avector_flag.ne.0) then
         call myclock(itime)
         do imc=1,mcscfnum
            call cmf_prop_avector(yyy%cmfavec(:,imc,1),  &
                 yyy%cmfavec(:,imc,0), linearflag,tin,tout,imc,qq)
            numaiters=numaiters+qq
         enddo
         call myclock(jtime);     times(5)=times(5)+jtime-itime;  itime=jtime
!! prevent drift
         if (par_consplit.eq.0) then
            call mympibcast(yyy%cmfavec(:,:,0),1,tot_adim*mcscfnum)
            call myclock(jtime);     times(9)=times(9)+jtime-itime
         endif
      endif
    end subroutine cmf_prop_avector_and_stuff

    subroutine end_stuff()
      implicit none
      if (constraintflag.ne.0) then
         call myclock(itime)
         call get_constraint(tout)
         call myclock(jtime); times(7)=times(7)+jtime-itime
      endif

      if (drivingflag.ne.0) then
         call myclock(itime)
         call drivingtrans(tout)
         call myclock(jtime); times(8)=times(8)+jtime-itime
      endif

      call myclock(itime)
      if (use_fockmatrix) then
         call get_fockmatrix()
      endif
      call get_reducedpot()
      if (numfrozen.gt.0) then
         call get_frexchange()
      endif
      call myclock(jtime);     times(3)=times(3)+jtime-itime
    end subroutine end_stuff

  end subroutine cmf_prop_wfn

!!!   A-VECTOR PROPAGATION ROUTINE FOR CMF

!! now just cmf and lmf.

  subroutine cmf_prop_avector(avectorin,avectorout,linearflag,time1,time2,imc,numiters)
    use parameters
    use mpisubmod
    implicit none
    DATATYPE,intent(in) :: avectorin(tot_adim)
    DATATYPE,intent(out) :: avectorout(tot_adim)
    integer,intent(in) :: linearflag,imc
    integer,intent(out) :: numiters
    real*8,intent(in) :: time1,time2
    DATATYPE :: csum
    DATATYPE,allocatable :: tempvector(:)
    integer :: k,qq
    real*8 :: timea,timeb

    if (avector_flag.eq.0) then
       if (tot_adim.gt.0) then
          avectorout(:)=avectorin(:)
       endif
       return
    endif

    allocate(tempvector(tot_adim))
    if (tot_adim.gt.0) then
       tempvector(:)=avectorin(:)
    endif
    numiters=0
    do k=1,littlesteps
       timea=time1+(time2-time1)*(k-1)/littlesteps
       timeb=time1+(time2-time1)*k/littlesteps

       call cmf_prop_avector0(tempvector,avectorout,linearflag,timea,timeb,imc,qq)
       if (tot_adim.gt.0) then
          tempvector(:)=avectorout(:)
       endif
       numiters=numiters+qq
    enddo

    if (threshflag.ne.0) then
       csum=0d0
       if (tot_adim.gt.0) then
          csum=dot(avectorout,avectorout,tot_adim)
       endif
       if (par_consplit.ne.0) then
          call mympireduceone(csum)
       endif
       if (tot_adim.gt.0) then
          avectorout(:)=avectorout(:)/sqrt(csum)
       endif
    endif

    deallocate(tempvector)

  contains
    subroutine cmf_prop_avector0(avectorin,avectorout,linearflag,time1,time2,imc,numiters)
      use linearmod
      use parameters
      use mpimod
      use xxxmod
      use configmod
      use configpropmod
      use pulsesubmod
      use configstuffmod
      implicit none
      DATATYPE,intent(in) :: avectorin(tot_adim)
      DATATYPE,intent(out) :: avectorout(tot_adim)
      integer,intent(in) :: imc,linearflag
      integer,intent(out) :: numiters
      real*8, intent(in) :: time1,time2
      DATATYPE :: sum1,sum0,pots(3)=0d0
      DATATYPE :: nullvector1(numr),nullvector2(numr)
      integer :: itime,jtime,getlen,ii,iflag,myiostat
      real*8 :: thisstep,midtime,rsum
      integer, save :: times(2)=0, icalled=0

      nullvector1(:)=0;nullvector2(:)=0

      call myclock(itime)

      if (sparseopt.eq.0) then
         call zero_cptr(workconfigpointer)
      else
         call zero_sptr(worksparsepointer,www)
         if (use_dfwalktype) then
            if (shuffle_dfwalktype) then
               call zero_sptr(workfdsparsepointer,fdww)
            else
               call zero_sptr(workdfsparsepointer,dfww)
            endif
         endif
      endif

      thisstep=time2-time1;

!! 062714
      midtime=(time2+time1)/2d0

      iflag=0
      if (drivingflag.ne.0) then
         call vectdpot(midtime,velflag,pots,imc)
         rsum=0d0
         do ii=1,3
            rsum=rsum+abs(pots(ii))**2
         enddo
         if (rsum.eq.0d0) then
            iflag=0
         else
            iflag=1
         endif
      endif

      if (linearflag.eq.0) then

         if (sparseopt.eq.0) then
            call assign_cptr(workconfigpointer,yyy%cptr(0),thisstep*DATAONE)
         else
            call assign_sptr(worksparsepointer,yyysptr(0),thisstep*DATAONE,www)
            if (use_dfwalktype) then
               if (shuffle_dfwalktype) then
                  call assign_sptr(workfdsparsepointer,yyysfdptr(0),thisstep*DATAONE,fdww)
               else
                  call assign_sptr(workdfsparsepointer,yyysdfptr(0),thisstep*DATAONE,dfww)
               endif
            endif
         endif

         if (drivingflag.ne.0.and.tot_adim.gt.0) then
            workdrivingavec(:,:)=0d0

            if (iflag.eq.1) then
               workdrivingavec(:,:)=( yyy%drivingavectorsxx(:,:,imc,0)*pots(1) + &
                    yyy%drivingavectorsyy(:,:,imc,0)*pots(2) + &
                    yyy%drivingavectorszz(:,:,imc,0)*pots(3) )*thisstep
            endif
         endif
      else

!! MAJOR BUG 062714     midtime=(time2+time1)/2d0

         sum0=(midtime-firsttime)/(lasttime-firsttime)
         sum1=(lasttime-midtime)/(lasttime-firsttime)

         if (drivingflag.ne.0.and.tot_adim.gt.0) then
            workdrivingavec(:,:)=0d0
            if (iflag.eq.1) then
               workdrivingavec(:,:) = thisstep*( &
   (sum1*yyy%drivingavectorsxx(:,:,imc,1)+ sum0*yyy%drivingavectorsxx(:,:,imc,0))*pots(1) + &
   (sum1*yyy%drivingavectorsyy(:,:,imc,1)+ sum0*yyy%drivingavectorsyy(:,:,imc,0))*pots(2) + &
   (sum1*yyy%drivingavectorszz(:,:,imc,1)+ sum0*yyy%drivingavectorszz(:,:,imc,0))*pots(3) )
            endif
         endif

         if (sparseopt.eq.0) then
            call add_cptr(yyy%cptr(1),yyy%cptr(0),workconfigpointer,sum1*thisstep,sum0*thisstep)
         else
            call add_sptr(yyysptr(1),yyysptr(0),worksparsepointer,sum1*thisstep,sum0*thisstep,www)
            if (use_dfwalktype) then
               if (shuffle_dfwalktype) then
                  call add_sptr(yyysfdptr(1),yyysfdptr(0),workfdsparsepointer,sum1*thisstep,&
                       sum0*thisstep,fdww)
               else
                  call add_sptr(yyysdfptr(1),yyysdfptr(0),workdfsparsepointer,sum1*thisstep,&
                       sum0*thisstep,dfww)
               endif
            endif
         endif
      endif
      call myclock(jtime); times(1)=times(1)+jtime-itime; itime=jtime

      if (tot_adim.gt.0) then
         call myconfigprop(avectorin(:),avectorout(:),midtime,imc,numiters)
      else
         call myconfigprop(nullvector1(:),nullvector2(:),midtime,imc,numiters)
      endif
      call myclock(jtime); times(2)=times(2)+jtime-itime;

      if (notiming.eq.0.and.myrank.eq.1) then

         if (icalled.eq.0) then
            open(11766, file=timingdir(1:getlen(timingdir))//"/aprop.time.dat", &
                 status="unknown",iostat=myiostat)
            call checkiostat(myiostat," writing aprop.time.dat")
            write(11766,'(100A11)',iostat=myiostat)  "init", "prop"
            call checkiostat(myiostat," writing aprop.time.dat")
            close(11766)
         endif
         open(11766, file=timingdir(1:getlen(timingdir))//"/aprop.time.dat", &
              status="unknown", position="append",iostat=myiostat)
         call checkiostat(myiostat," writing aprop.time.dat")
         write(11766,'(100I11)',iostat=myiostat)  times(1:2)
         call checkiostat(myiostat," writing aprop.time.dat")
         close(11766)
      endif

      icalled=1

    end subroutine cmf_prop_avector0

  end subroutine cmf_prop_avector

end module prop_loop_sub_mod

module proploopmod
contains

subroutine prop_loop( starttime)
  use parameters
  use mpimod
  use xxxmod
  use configmod
  use derivativemod
  use sparsemultmod
  use basissubmod
  use savenormmod
  use mpisubmod
  use prop_loop_sub_mod
  use orbdermod
  use getstuffmod
  use loadstuffmod
  use spfsubmod
  implicit none
  integer ::  jj,flag, itime, jtime, times(20)=0, qq,imc,getlen,myiostat
  DATAECS :: thisenergy(mcscfnum), lastenergy(mcscfnum) ,thisenergyavg,&
       lastenergyavg,startenergy(mcscfnum)
  CNORMTYPE :: norms(mcscfnum)
  real*8 :: thistime, starttime, thattime,error=1d10,avecerror=1d10
  DATATYPE :: sum2,sum,drivingoverlap(mcscfnum)
  DATATYPE, allocatable :: avectorp(:),outspfs(:)

  thistime=starttime;  flag=0;    call zero_mpi_times()

  allocate(avectorp(tot_adim),outspfs(totspfdim))
  outspfs=0
  if (tot_adim.gt.0) then
     avectorp=0;
  endif

  lastenergy(:)=1.d+3;  lastenergyavg=1.d+3

  call myclock(itime)

  do imc=1,mcscfnum
     call basis_project(www,numr,yyy%cmfavec(:,imc,0))
  enddo

  if (normboflag.ne.0) then
     OFLWR "    Enforcing BO norms due to normboflag"; CFL
     call enforce_bonorms(mcscfnum,  yyy%cmfavec(:,:,0),savenorms(:,:))
  endif

  call get_stuff(0.0d0)

  do imc=1,mcscfnum

     call sparseconfigmult(www,yyy%cmfavec(:,imc,0),avectorp,&
          yyy%cptr(0),yyysptr(0),1,1,1,0,0d0,imc)

     sum=0;     sum2=0
     if (tot_adim.gt.0) then
        sum = dot(yyy%cmfavec(:,imc,0),yyy%cmfavec(:,imc,0),tot_adim)
        sum2=dot(yyy%cmfavec(:,imc,0),avectorp,tot_adim)
     endif
     if (par_consplit.ne.0) then
        call mympireduceone(sum);        call mympireduceone(sum2)
     endif
     OFLWR "IN PROP: VECTOR NORM ",sqrt(sum);CFL
     startenergy(imc)=sum2/sum
     OFLWR "         ENERGY ", startenergy(imc); CFL

!!$  moving to main
!!$     if (normboflag.ne.0) then
!!$        OFLWR "    ... will enforce BO norms due to normboflag"; CFL
!!$        call get_bonorms(1,yyy%cmfavec(:,imc,0),savenorms(:,imc))
!!$     endif

  enddo

  if (drivingflag.ne.0) then
     OFLWR "call drivinginit"; CFL
     call drivinginit(startenergy)
     call get_stuff(0.0d0)
     OFLWR "called drivinginit"; CFL
  endif

  if (debugflag.eq.956) then
     OFLWR "Stopping due to debugflag=956"; CFLST
  endif

  if ((myrank.eq.1).and.(notiming.le.1)) then
     call system("echo -n > "//timingdir(1:getlen(timingdir))//"/abstiming.dat")
     open(853, file=timingdir(1:getlen(timingdir))//"/Main.time.dat", &
          status="unknown",iostat=myiostat)
     call checkiostat(myiostat," opening Main.time.dat")
     write(853,'(T16,100A15)')  &
          "Prop ", &     !! (1)
          "Act ", &      !! (2)
          "After ", &    !! (3)
          "Init", &      !! (4)
          "Save",&       !! (5)
          "MPI",&        !! (6)
          "Non MPI"      !! (7)
     close(853)
  endif

  call myclock(jtime)  ;  times(4)=times(4)+jtime-itime

  jj=0
  do while (flag==0)      !! BEGIN MAIN LOOP
     jj=jj+1

     call myclock(itime)

     if (save_every.ne.0.and.mod(jj,save_every).eq.0) then
        call save_vector(yyy%cmfavec(:,:,0),yyy%cmfspfs(:,0),avectoroutfile,spfoutfile)
     endif

     if ((jj==numpropsteps) .and. (threshflag/=1)) then
        flag=1
     endif

     call checkstopfile()

     thattime=thistime+par_timestep

     call myclock(jtime)  ;  times(5)=times(5)+jtime-itime;    itime=jtime

                                 !! ************************************************************* !!
     call actionsub( thistime)   !! ACTIONS.  IF ACTION CHANGES PSI THEN CALL GET_STUFF AFTER ACTION.
                                 !! ************************************************************* !!

     call myclock(jtime)  ;     times(2)=times(2)+jtime-itime;     itime=jtime

!! norms not gotten yet, tempfix 08-2016 for odex
!     rsum=0d0
!     do imc=1,mcscfnum
!        rsum=rsum+norms(imc)**2
!     enddo
!     rsum=sqrt(rsum);     abserr=sqrt(rsum)*myrelerr
!! tempfix:
     abserr=myrelerr

     if ((cmf_flag==1)) then
        call cmf_prop_wfn(thistime, thattime)
     else
        call prop_wfn(thistime, thattime)
     endif

     do imc=1,mcscfnum
        call basis_project(www,numr,yyy%cmfavec(:,imc,0))
     enddo

     if (normboflag.ne.0) then
!!        OFLWR "    ... enforcing BO norms due to normboflag"; CFL
        call enforce_bonorms(mcscfnum,  yyy%cmfavec(:,:,0),savenorms(:,:))
     endif

     call myclock(jtime)  ;  times(1)=times(1)+jtime-itime;    itime=jtime

     do imc=1,mcscfnum

        call sparseconfigmult(www,yyy%cmfavec(:,imc,0),avectorp,&
             yyy%cptr(0),yyysptr(0),1,1,timedepexpect,0,thattime,imc)

        call basis_project(www,numr,avectorp)

        sum=0; sum2=0
        if (tot_adim.gt.0) then
           sum=dot(yyy%cmfavec(:,imc,0),avectorp(:),tot_adim)
           sum2=dot(yyy%cmfavec(:,imc,0),yyy%cmfavec(:,imc,0),tot_adim)
        endif
        if (par_consplit.ne.0) then
           call mympireduceone(sum); call mympireduceone(sum2)
        endif

        thisenergy(imc) = sum/sum2   !! ok conversion
        norms(imc)=sqrt(sum2)   !! ok conversion

        OFL
#ifdef ECSFLAG
#ifdef CNORMFLAG
        write(mpifileptr,'(A3,F16.5, 2(A10, 2E18.10),E18.10)') "T= ",thattime,&
             "Energy: ", thisenergy(imc), "Norm: ", abs(norms(imc)),norms(imc)
#else
        write(mpifileptr,'(A3,F16.5, 2(A10, 2E18.10))') "T= ",thattime, &
             "Energy: ", thisenergy(imc), "Norm: ", norms(imc)
#endif
#else
        write(mpifileptr,'(A3,F16.5, 2(A10, E18.10))') "T= ", thattime, &
             "Energy: ", thisenergy(imc), "Norm: ", norms(imc)
#endif
        CFL

        if (drivingflag.ne.0) then
          call getdrivingoverlap(drivingoverlap,mcscfnum)

          OFL
#ifdef ECSFLAG
           write(mpifileptr,'(A3,F16.5, 2(A10, 2E18.10))') "t= ",thattime, "DEnergy ", &
#else
           write(mpifileptr,'(A3,F16.5, 2(A10, E18.10))') "t= ", thattime, "DEnergy ", &
#endif

    (    (drivingenergies(imc)*drivingproportion**2 + CONJUGATE(drivingenergies(imc) * drivingoverlap(imc)) + &
    drivingenergies(imc)*drivingoverlap(imc) + thisenergy(imc)*norms(imc)**2)     )/  &
    (drivingproportion**2 + drivingoverlap(imc) + CONJUGATE(drivingoverlap(imc)) + norms(imc)**2), &
    "DNorm ", (sqrt &
    (drivingproportion**2 + drivingoverlap(imc) + CONJUGATE(drivingoverlap(imc)) + norms(imc)**2))

           CFL

!! ok so < Psi(t) | H | Psi(t) > =

!!   < Psi_0 | H | Psi_0 > +          drivingenergy * drivingproportion**2
!!   < Psi_0 | H | Psi' > +          = CONJUGATE(drivingenergy) * < Psi_0 | Psi' >  !! factor in psi0 in code * drivingproportion
!!   < Psi' | H | Psi_0 > +          = drivingenergy * < Psi_0 | Psi' >  !! factor in psi0 in code * drivingproportion
!!   < Psi' | H | Psi' >               in code

!! ok so < Psi(t) | Psi(t) > =

!!   < Psi_0 | Psi_0 > +          drivingproportion**2
!!   < Psi_0 |  Psi' > +          = < Psi_0 | Psi' >    !! factor in psi_0 in code * drivingproportion
!!   < Psi' | Psi_0 > +          = < Psi' | Psi_0 >     !!  ditto   * drivingproportion
!!   < Psi' | Psi' >               in code

        endif

     enddo  !! imc

     thisenergyavg=0
     do imc=1,mcscfnum
        thisenergyavg=thisenergyavg+thisenergy(imc)/mcscfnum
     enddo

     if ((myrank.eq.1).and.(notiming.le.1)) then
        call system("date >> "//timingdir(1:getlen(timingdir))//"/abstiming.dat")
        open(853, file=timingdir(1:getlen(timingdir))//"/Main.time.dat", &
             status="old", position="append",iostat=myiostat)
        call checkiostat(myiostat," opening Main.time.dat")
        write(853,'(F13.3,T16,100I15)',iostat=myiostat)  thistime, &
             times(1:5)/1000, mpitime/1000, nonmpitime/1000
        call checkiostat(myiostat," writing Main.time.dat")
        close(853)
     endif

     if (debugflag.eq.42) then
        OFLWR "Stopping due to debugflag=42"; CFLST
     endif

     thistime=thattime

     if (threshflag.eq.1) then

!! 06-16 why change to actreduced
!! not spf_linear_derivs?  need exchange
!!        call actreduced0(1,0d0,yyy%cmfspfs(:,0),yyy%cmfspfs(:,0),outspfs,0,1,1)

        call spf_linear_derivs0(0,1,0d0,yyy%cmfspfs(:,0),outspfs,1,1)

        call apply_spf_constraints(outspfs)

        avecerror=abs(thisenergyavg-lastenergyavg)/par_timestep

!!$ OOPS 04-24-15        error=sqrt(abs(hermdot(outspfs,outspfs,totspfdim)))
!                        if (parorbsplit.eq.3) then
!                           call mympirealreduceone(error)
!                        endif

        error=abs(hermdot(outspfs,outspfs,totspfdim))
        if (parorbsplit.eq.3) then
           call mympirealreduceone(error)
        endif
        error=sqrt(error)

        OFL; write(mpifileptr,'(A24,2E10.2,A10,2E10.2)') &
             " STOPTEST : ORBITALS ", error,stopthresh, " AVECTOR ",avecerror,astoptol;CFL

        if ( error.lt.stopthresh .or.spf_flag.eq.0) then
           if (avecerror.gt.astoptol) then
              OFL; write(mpifileptr,'(A67,2E12.5)') &
                   "   Orbitals Converged, Avector not necessarily converged "; CFL
           else
              flag=1
              OFLWR; WRFL "   ***  CONVERGED *** "; WRFL
              do qq=1,mcscfnum
                 write(mpifileptr,'(T5,10F20.12)') thisenergy(qq), thisenergy(qq)-lastenergy(qq)
              enddo
              write(mpifileptr, *) "   ***   ";              write(mpifileptr, *);CFL
           endif
        endif
        do qq=1,mcscfnum
           lastenergy(qq)=thisenergy(qq)
        enddo
        lastenergyavg=thisenergyavg

        par_timestep=par_timestep*timestepfac
        if (par_timestep.gt.max_timestep) then
           par_timestep=max_timestep
        endif

     endif

     call myclock(jtime)  ;        times(3)=times(3)+jtime-itime

  enddo    !! END MAIN LOOP

  norms=0
  if (tot_adim.gt.0) then
     do imc=1,mcscfnum
        norms(imc)=dot(yyy%cmfavec(:,imc,0),& !! ok conversion
             yyy%cmfavec(:,imc,0),tot_adim)   !! ok conversion
     enddo
  endif
  if (par_consplit.ne.0) then
#ifndef REALGO
#ifndef CNORMFLAG
     call mympirealreduce(norms,mcscfnum)
#else
     call mympireduce(norms,mcscfnum)
#endif
#else
     call mympireduce(norms,mcscfnum)
#endif
  endif
  norms(:)=sqrt(norms(:))

  if (threshflag.eq.1) then
     if (tot_adim.gt.0) then
        do imc=1,mcscfnum
           yyy%cmfavec(:,imc,0)=yyy%cmfavec(:,imc,0)/norms(imc)
        enddo
     endif
     norms(:)=1d0
  endif

  call actions_final()

  OFLWR "   ...done prop..."; CFL
  call mpibarrier()

  if (saveflag.ne.0) then
     OFLWR "  ...saving vector..."; CFL
     call mpibarrier()
     call save_vector(yyy%cmfavec(:,:,0),yyy%cmfspfs(:,0),avectoroutfile,spfoutfile)
  endif
  OFLWR "   ...end prop..."; CFL
  call mpibarrier()

  deallocate(avectorp,outspfs)

end subroutine prop_loop

end module proploopmod


!========================================================================
!
!                   S P E C F E M 2 D  Version 7 . 0
!                   --------------------------------
!
!     Main historical authors: Dimitri Komatitsch and Jeroen Tromp
!                        Princeton University, USA
!                and CNRS / University of Marseille, France
!                 (there are currently maNZ_IMAGE_color more authors!)
! (c) Princeton University and CNRS / University of Marseille, April 2014
!
! This software is a computer program whose purpose is to solve
! the two-dimensional viscoelastic anisotropic or poroelastic wave equation
! using a spectral-element method (SEM).
!
! This software is governed by the CeCILL license under French law and
! abiding by the rules of distribution of free software. You can use,
! modify and/or redistribute the software under the terms of the CeCILL
! license as circulated by CEA, CNRS and Inria at the following URL
! "http://www.cecill.info".
!
! As a counterpart to the access to the source code and rights to copy,
! modify and redistribute granted by the license, users are provided only
! with a limited warranty and the software's author, the holder of the
! economic rights, and the successive licensors have only limited
! liability.
!
! In this respect, the user's attention is drawn to the risks associated
! with loading, using, modifying and/or developing or reproducing the
! software by the user in light of its specific status of free software,
! that may mean that it is complicated to manipulate, and that also
! therefore means that it is reserved for developers and experienced
! professionals having in-depth computer knowledge. Users are therefore
! encouraged to load and test the software's suitability as regards their
! requirements in conditions enabling the security of their systems and/or
! data to be ensured and, more generally, to use and operate it in the
! same conditions as regards security.
!
! The full text of the license is available in file "LICENSE".
!
!========================================================================
  subroutine write_postscript_snapshot()

  use specfem_par, only: myrank,p_sv,it,imagetype_postscript, &
                         potential_acoustic,potential_gravitoacoustic, &
                         potential_gravito,displ_elastic,displs_poroelastic, &
                         potential_dot_acoustic,potential_dot_gravitoacoustic, &
                         potential_dot_gravito,veloc_elastic,velocs_poroelastic, &
                         potential_dot_dot_acoustic,potential_dot_dot_gravitoacoustic, &
                         potential_dot_dot_gravito,accel_elastic,accels_poroelastic

  implicit none
  include "constants.h"

  if( myrank == 0 ) then
    write(IOUT,*)
    write(IOUT,*) 'Writing PostScript vector plot for time step ',it
  endif

  if( imagetype_postscript == 1 .and. p_sv ) then

    if( myrank == 0 ) write(IOUT,*) 'drawing displacement vector as small arrows...'
    call compute_vector_whole_medium(potential_acoustic,potential_gravitoacoustic, &
                 potential_gravito,displ_elastic,displs_poroelastic)

    call plotpost()

  else if( imagetype_postscript == 2 .and. p_sv ) then

    if( myrank == 0 ) write(IOUT,*) 'drawing velocity vector as small arrows...'
    call compute_vector_whole_medium(potential_dot_acoustic,potential_dot_gravitoacoustic, &
                 potential_dot_gravito,veloc_elastic,velocs_poroelastic)

    call plotpost()

  else if( imagetype_postscript == 3 .and. p_sv ) then

    if( myrank == 0) write(IOUT,*) 'drawing acceleration vector as small arrows...'
    call compute_vector_whole_medium(potential_dot_dot_acoustic,potential_dot_dot_gravitoacoustic, &
                 potential_dot_dot_gravito,accel_elastic,accels_poroelastic)

    call plotpost()

  else if( .not. p_sv ) then
    call exit_MPI('cannot draw a SH scalar field as a vector plot, turn PostScript plots off')
  else
    call exit_MPI('wrong type for PostScript snapshots')
  endif

  if( myrank == 0 .and. imagetype_postscript /= 4 .and. p_sv ) write(IOUT,*) 'PostScript file written'

  end subroutine write_postscript_snapshot


text:
PROGRAM gs_comp

  use gs
  use gentests

  implicit none

  ! PARAMETERS
  integer, parameter :: MAXIT = 100000
  real(8), parameter :: TOL = 1.0D-5

  ! LOCAL ARRAYS
  real(8), allocatable :: M(:,:),b(:),x(:),rv(:),sol(:)
  integer, allocatable :: seed(:)

  ! LOCAL SCALARS
  real(8) :: r,mr,msdiff,mt,miter
  integer :: i,k,n,status,rsize,t1,t2,clock_rate,clock_max,iter, &
       trials,strials,nfailed,type

  write(*,*) 'Type dimension and number of trials...'
  read(*,*) n,trials
  write(*,*) '1- Diagonal Dominant 2- Stochastic Least Squares'
  write(*,*) '3- Symmetric Definite Positive'
  write(*,*) '4- Symmetric Semi-Definite Positive'
  write(*,*) '5- Symmetric Indefinite'
  read(*,*) type

  allocate(M(n,n),b(n),x(n),rv(n),sol(n),STAT=status)

  if (status .ne. 0) then
     write(*,*) "Memory problems"
     stop
  end if

  CALL RANDOM_SEED(SIZE=rsize)

  allocate(seed(rsize))
  do i = 1,rsize
     seed(i) = 123456789
  end do

  write(*,*) "Starting serial..."

  CALL RANDOM_SEED(PUT=seed)

  nfailed = 0
  mt = 0.0D0
  mr = 0.0D0
  miter = 0.0D0
  msdiff = 0.0D0

  do k = 1,trials

     if (type .eq. 1) then
        call genDDRandSystem(n,M,b,sol)
     else if (type .eq. 2) then
        call genRandStochSystem(n,M,b,sol)
     else if (type .eq. 3) then
        call genRandSDPSystem(n,M,b,sol,0.0D0,1.0D0 * n)
     else if (type .eq. 4) then
        call genRandSDPSystem(n,M,b,sol,4.0D-1,1.0D0 * n)
     else if (type .eq. 5) then
        call genRandSDPSystem(n,M,b,sol,2.0D-1,- 1.0D0 * n)
     else
        write(*,*) 'Invalid type.'
        exit
     end if

     do i = 1,n
        x(i) = 0.0D0
     end do

     CALL SYSTEM_CLOCK(t1,clock_rate,clock_max)

     CALL serialgs(n,M,b,x,r,iter,TOL,MAXIT)

     CALL SYSTEM_CLOCK(t2,clock_rate,clock_max)

     if (r .le. TOL) then
        miter = miter + iter
        mt = mt + (t2 - t1)
        mr = mr + r
        msdiff = msdiff + maxval(abs(x - sol))
     else
        nfailed = nfailed + 1
     end if

  end do

  strials = trials - nfailed

  write(*,9000) mr / (1.0D0 * strials),mt / (1.0D0 * clock_rate * strials), &
       INT(miter / (1.0D0 * strials)),msdiff / (1.0D0 * strials), &
       (1.0D+2 * nfailed) / trials

  write(*,*) "Starting parallel..."

  CALL RANDOM_SEED(PUT=seed)

  nfailed = 0
  mt = 0.0D0
  mr = 0.0D0
  miter = 0.0D0
  msdiff = 0.0D0

  do k = 1,trials

     if (type .eq. 1) then
        call genDDRandSystem(n,M,b,sol)
     else if (type .eq. 2) then
        call genRandStochSystem(n,M,b,sol)
     else if (type .eq. 3) then
        call genRandSDPSystem(n,M,b,sol,0.0D0,1.0D0 * n)
     else if (type .eq. 4) then
        call genRandSDPSystem(n,M,b,sol,4.0D-1,1.0D0 * n)
     else if (type .eq. 5) then
        call genRandSDPSystem(n,M,b,sol,2.0D-1,- 1.0D0 * n)
     else
        write(*,*) 'Invalid type.'
        exit
     end if

     do i = 1,n
        x(i) = 0.0D0
     end do

     CALL SYSTEM_CLOCK(t1,clock_rate,clock_max)

     CALL sharedgs(n,M,b,x,r,iter,TOL,MAXIT)

     CALL SYSTEM_CLOCK(t2,clock_rate,clock_max)

     if (r .le. TOL) then
        miter = miter + iter
        mt = mt + (t2 - t1)
        mr = mr + r
        msdiff = msdiff + maxval(abs(x - sol))
     else
        nfailed = nfailed + 1
     end if

  end do

  strials = trials - nfailed

  write(*,9000) mr / (1.0D0 * strials),mt / (1.0D0 * clock_rate * strials), &
       INT(miter / (1.0D0 * strials)),msdiff / (1.0D0 * strials), &
       (1.0D+2 * nfailed) / trials

!  write(*,2000) x
!  write(*,2000) sol

!2000 FORMAT(/,'Point:',/,3(1X,E20.10))
9000 FORMAT('RESIDUAL:',1X,D12.5,1X,'TIME (in secs):',1X,F10.4,&
          1X,'ITERATIONS:',1X,I8,1X,'DIFF',1X,D12.5,1X,&
          'FAILED:',1X,F6.2,'%')

END PROGRAM gs_comp

max_stars_repo_name: fsobral/ccs

max_issues_repo_name: fsobral/ccs

max_forks_repo_name: fsobral/ccs