program test_simsphere
  use, intrinsic :: ieee_arithmetic
  use simsphere_mod
  use mod_testing, only: assert,initialize_tests,report_tests
  implicit none

  ! mod_testing test setup

  logical, dimension(:), allocatable :: tests
  logical :: test_failed
  integer :: n, ntests

  ! end mod_testing test setup

  integer :: i

  logical, parameter :: splint_test = .true.
  logical, parameter :: spline_test = .true.
  logical, parameter :: transm_ftabs_test = .true.
  logical, parameter :: transm_ftabs_test2 = .true.
  logical, parameter :: transm_ftscatT_test = .true.
  logical, parameter :: transm_ftscatT_test2 = .true.
  logical, parameter :: transm_fbscatT_test = .true.
  logical, parameter :: transm_fbscatT_test2 = .true.
  logical, parameter :: advect_func_test = .true.
!  logical :: air_test
  logical, parameter :: cond_test = .true.
  logical, parameter :: stomfs_test = .true.
  logical, parameter :: stomrs_test_hi_temp = .true.
  logical, parameter :: stomrs_test_lo_temp = .true.
  logical, parameter :: stomrs_test_hi_psi = .true.
  logical, parameter :: stomrs_test_lo_psi = .true.
  logical, parameter :: stomc_test = .true.
  logical, parameter :: you_star_test = .true.
  logical, parameter :: r_ohms_test = .true.
  logical, parameter :: windf_test = .true.
  logical, parameter :: stab_test = .true.
  logical, parameter :: stabh_test = .true.
  logical, parameter :: fstabh_test = .true.
  logical, parameter :: fstabm_test = .true.
  logical, parameter :: restrn_test = .true.
  logical, parameter :: psgcal_test = .true.
  logical, parameter :: vegflx_test = .true.
  logical, parameter :: avr_test_init1 = .true.
!  logical, parameter :: avr_test_init_not1 = .true.
  logical, parameter :: ozone_test = .true.
  logical, parameter :: co2flx_test = .true.
  logical, parameter :: veghot_test = .true.
  logical, parameter :: hot_test = .true.
  logical, parameter :: albedo_test = .true.
  logical, parameter :: intpol_test = .true.
  logical, parameter :: prfile_test = .true.
  logical, parameter :: f_control_test = .false.

! Tolerance for comparisons
  real, parameter :: eq_tol = 1e-6

! splint_test variables
  integer, parameter :: splint_max_array = 50
  real, parameter :: splint_expected = 1.5
  real :: splint_arg1(splint_max_array), splint_arg2(splint_max_array)
  real(kind=4), allocatable :: splint_arg3(:)
  integer :: splint_arg4, splint_arg5
  real :: splint_output

! spline_test variables
  integer, parameter :: spline_max_array = 50
  integer :: spline_arg3
  real, parameter :: spline_expected(62) =                            &
    (/-8.84081032E-07,1.76816206E-06,-6.18856711E-06,2.29861053E-05,  &
      -8.57558480E-05,3.20037303E-04,-1.19439338E-03,4.45753615E-03,  &
      -1.66357514E-02,6.20854683E-02,-6.20854646E-02,1.66357514E-02,  &
      -4.45753615E-03,1.19439315E-03,-3.20036663E-04,8.57535633E-05,  &
      -2.29775978E-05,6.15682848E-06,-1.64971721E-06,4.42040374E-07,  &
      -1.18444362E-07,3.17370699E-08,-8.50392201E-09,2.27861907E-09,  &
      -6.10554107E-10,1.63597469E-10,-4.38358100E-11,1.17457693E-11,  &
      -3.14726977E-12,8.43309526E-13,-2.25968590E-13,6.05648480E-14,  &
      -1.62907932E-14,4.59832376E-15,-2.10250179E-15,3.81168384E-15,  &
      -1.31442336E-14,4.87652521E-14,-1.81916778E-13,6.78901867E-13,  &
      -2.53369049E-12,9.45585971E-12,-3.52897503E-11,1.31703135E-10,  &
      -4.91522767E-10,1.83438786E-09,-6.84602863E-09,2.55497259E-08,  &
      -9.53528740E-08,3.55861772E-07,-1.32809419E-06,4.95651511E-06,  &
      -1.84979654E-05,6.90353481E-05,-2.57643434E-04,9.61538346E-04,  &
      -3.58850998E-03,1.33925015E-02,-4.99814972E-02,0.186533481,     &
      -0.696152449,2.59807634/)
  real :: spline_arg1(spline_max_array), spline_arg2(spline_max_array)
  real :: spline_arg4, spline_arg5
  real, allocatable :: spline_output(:)

! transm_ftabsT variables
  real :: transm_ftabs_arg1
  real, parameter :: transm_ftabs_expected = 0.722923875
  real, parameter :: transm_ftabs_expected2 = 0.437107921

! transm_ftscatT variables
  real :: transm_ftscatT_arg1
  real, parameter :: transm_ftscatT_expected = 0.553545594
  real, parameter :: transm_ftscatT_expected2 = 0.173223004

! transm_fbscatT variables
  real :: transm_fbscatT_arg1
  real, parameter :: transm_fbscatT_expected = 0.307550341
  real, parameter :: transm_fbscatT_expected2 = 0.304252803

! netrad_init_test variables
!  real, parameter :: aepsi_expected = 0.804384410
!  real, parameter :: aepsi_cloud_flag_expected = 1.02347386

! advect_test variables
  real, parameter :: advect_test_expected = 9.41291146E-05

! air_test variables
!  real, parameter :: air_test_expected = 0.0

! cond_test variables
  real, parameter :: cond_test_expected = 1.46884334

! stomfs_test variables
  real, parameter :: stomfs_test_expected = 1.58197677

! stomrs_test_* variables
  real, parameter :: stomrs_test_hi_temp_expected = 5000.0
  real, parameter :: stomrs_test_lo_temp_expected = 5000.0
  real, parameter :: stomrs_test_hi_psi_expected = 3.16395354
  real, parameter :: stomrs_test_lo_psi_expected = 7.90988398

! stomc_test variables
  real, parameter :: stomc_test_expected = 4.74593019

! you_star_test variables
  real, parameter :: you_star_test_expected = 1.73396431E-02
  real :: you_star_arg1, you_star_arg2, you_star_arg3, you_star_arg4, you_star_arg5

! r_ohms_test variables
  real, parameter :: r_ohms_test_expected = 42.6767731
  real :: r_ohms_arg1, r_ohms_arg2, r_ohms_arg3, r_ohms_arg4, r_ohms_arg5

! windf_test variables
  real, parameter :: windf_test_expected = 8.98633766
  real :: windf_arg1, windf_arg2, windf_arg3, windf_arg4, windf_arg5

! stab_test variables
  real, parameter :: stab_test_expected = 0.707106769
  real :: stab_arg1, stab_arg2

! stabh_test variables
  real, parameter :: stabh_test_expected = 0.741619825

! fstabh_test variables
  real, parameter :: fstabh_test_expected = -0.810930133
  real :: fstabh_arg1, fstabh_arg2

! fstabm_test variables
 real, parameter :: fstabm_test_expected = -0.572894096

! restrn_test variables
  real, parameter :: restrn_test_expected = 0.590972006
  real :: restrn_arg1, restrn_arg2, restrn_arg3, restrn_arg4

! psgcal_test variables
  real, parameter :: psgcal_test_expected = -2.14821539E-05

! vegflx_test variables
  real, parameter :: vegflx_test_expected = 0.986029029
  real :: vegflx_arg1

! avr_test_* variables
  real, parameter :: avr_test_init1_expected = 300.5
!  real, parameter :: avr_test_init_not1_expected = 20.0
  real :: avr_arg1, avr_arg2
  integer :: avr_arg3

! co2flx_test variables
  real, parameter :: co2flx_ccan_expected = 1.77409637
  real, parameter :: co2flx_fco2_expected = 3.13184571

! veghot_test variables
  real, parameter :: veghot_heatv_expected = 1.5
  real :: veghot_arg1, veghot_arg2

! hot_test variables
  real, parameter :: hot_caseI_expected = 2.49819231
  real, parameter :: hot_caseII_expected = 1.62364423
  real, parameter :: hot_caseIII_expected = 0.999096155
  real :: hot_arg1, hot_arg2, hot_arg3, hot_arg4

! albedo_test variables
  real, parameter :: albedo_caseI_expected = 0.396643847
  real, parameter :: albedo_caseII_expected = 0.396643847
  real, parameter :: albedo_caseIII_expected = 0.396643847
  real, parameter :: albedo_caseIV_expected = 0.396643847
  real :: albedo_test_arg1

! intpol_test variables
  real, dimension(50), parameter :: intpol_ug_exp =  &
           (/9.18181801,8.36363602,7.54545498,6.72727299,5.90909100, &
             5.09090948,4.27272749,3.45454597,2.63636398,1.81818199, &
             1.00000000,0.181818962,-0.636363029,-1.45454502,-2.27272701, &
             -3.09090805,-3.90909004,-4.72727203,-5.54545403,-6.36363602, &
             0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
             0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
             0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
             0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
             0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
             0.00000000,0.00000000,0.00000000,0.00000000,0.00000000/)
  real, dimension(50), parameter :: intpol_vg_exp =  &
           (/9.18181801,8.36363602,7.54545498,6.72727299,5.90909100, &
             5.09090948,4.27272749,3.45454597,2.63636398,1.81818199, &
             1.00000000,0.181818962,-0.636363029,-1.45454502,-2.27272701, &
             -3.09090805,-3.90909004,-4.72727203,-5.54545403,-6.36363602, &
             0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
             0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
             0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
             0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
             0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
             0.00000000,0.00000000,0.00000000,0.00000000,0.00000000/)
  real, dimension(50), parameter :: intpol_ugd_exp =  &
            (/9.57142830,7.42857170,5.28571463,3.14285755,1.00000000, &
              1.00000000,1.00000000,1.00000000,1.00000000,1.00000000, &
              1.00000000,1.00000000,1.00000000,1.00000000,1.00000000, &
              1.00000000,1.00000000,1.00000000,1.00000000,1.00000000, &
              0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
              0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
              0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
              0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
              0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
              0.00000000,0.00000000,0.00000000,0.00000000,0.00000000/)
  real, dimension(50), parameter :: intpol_vgd_exp =  &
            (/9.57142830,7.42857170,5.28571463,3.14285755,1.00000000, &
              1.00000000,1.00000000,1.00000000,1.00000000,1.00000000, &
              1.00000000,1.00000000,1.00000000,1.00000000,1.00000000, &
              1.00000000,1.00000000,1.00000000,1.00000000,1.00000000, &
              0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
              0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
              0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
              0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
              0.00000000,0.00000000,0.00000000,0.00000000,0.00000000, &
              0.00000000,0.00000000,0.00000000,0.00000000,0.00000000/)

! prfile_test variables
  real, dimension(51), parameter :: prfile_uf_exp = (/1.00000000,         &
              0.800000012,0.600000024,0.399999976,0.199999988,0.00000000, &
              0.200000003,0.400000006,0.600000024,0.800000012,1.00000000, &
              0.800000012,0.600000024,0.399999976,0.199999988,0.00000000, &
              0.200000003,0.400000006,0.600000024,0.800000012,1.00000000, &
              1.00000000,1.00000000,1.00000000,1.00000000,1.00000000,     &
              1.00000000,1.00000000,1.00000000,1.00000000,1.00000000,     &
              1.00000000,1.00000000,1.00000000,1.00000000,1.00000000,     &
              1.00000000,1.00000000,1.00000000,1.00000000,1.00000000,     &
              1.00000000,1.00000000,1.00000000,1.00000000,1.00000000,     &
              0.800000012,0.600000024,0.399999976,0.199999988,0.00000000/)
  real, dimension(51), parameter :: prfile_vf_exp = (/1.00000000,         &
              0.800000012,0.600000024,0.399999976,0.199999988,0.00000000, &
              0.200000003,0.400000006,0.600000024,0.800000012,1.00000000, &
              0.800000012,0.600000024,0.399999976,0.199999988,0.00000000, &
              0.200000003,0.400000006,0.600000024,0.800000012,1.00000000, &
              1.00000000,1.00000000,1.00000000,1.00000000,1.00000000,     &
              1.00000000,1.00000000,1.00000000,1.00000000,1.00000000,     &
              1.00000000,1.00000000,1.00000000,1.00000000,1.00000000,     &
              1.00000000,1.00000000,1.00000000,1.00000000,1.00000000,     &
              1.00000000,1.00000000,1.00000000,1.00000000,1.00000000,     &
              0.800000012,0.600000024,0.399999976,0.199999988,0.00000000/)
  real, dimension(51), parameter :: prfile_tf_exp = (/1.00000000,0.800000012, &
        0.600000024,0.399999976,0.199999988,0.00000000,0.00000000,0.00000000, &
        0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,    &
        0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,    &
        0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,    &
        0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,    &
        0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,    &
        0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,    &
        0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,    &
        0.00000000/)
  real, dimension(51), parameter :: prfile_qf_exp = (/1.00000000,0.800000012, &
        0.600000024,0.399999976,0.199999988,0.00000000,0.00000000,0.00000000, &
        0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,    &
        0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,    &
        0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,    &
        0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,    &
        0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,    &
        0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,    &
        0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,    &
        0.00000000/)
  real, parameter :: prfile_awind_exp = 1.41421354

! f_control variables
  integer, parameter :: iyr_expected = 89
  integer, parameter :: imo_expected = 8
  integer, parameter :: iday_expected = 4
  real, parameter :: tz_expected = 6
  real, parameter :: xlat_expected = 39.25
  real, parameter :: xlong_expected = 96.34
  real, parameter :: strtim_expected = 0530
  real, parameter :: timend_expected = 2330
  real, parameter :: outtt_expected = 30
  real, parameter :: slope_expected = 0
  real, parameter :: aspect_expected = 0
  real, parameter :: f_expected = 0.5
  real, parameter :: fsub_expected = 0.75
  real, parameter :: wmax_expected = 0.34
  real, parameter :: btemp_expected = 24.63
  real, parameter :: tp_expected = 13
  character, parameter :: dual_ti_expected = 'Y'
  real, parameter :: ti_a_expected = 12.0
  real, parameter :: ti_b_expected = 12.0
!  character, parameter :: albedo_gflag_expected = 'F'
  real, parameter :: albg_expected = 0.1
  real, parameter :: epsi_expected = 0.96
!  integer, parameter :: index_soils_expected = 5
  real, parameter :: omega_expected = 3.13
  real, parameter :: zo_expected = 0.05
!  real, parameter :: obst_hgt_expected = 1
!  logical, parameter :: cloud_flag_expected = .true.
  real, parameter :: cld_fract_expected = 1.4
  real, parameter :: frveg_expected = 1
  real, parameter :: xlai_expected = 0
  real, parameter :: epsf_expected = 0.96
!  character, parameter :: albedo_fflag_expected = 'N'
  real, parameter :: albf_expected = 0.1
  character, parameter :: stmtype_expected =  'L'
!  integer, parameter :: index_veggies_expected = 4
  real, parameter :: volrel_expected = 10
  real, parameter :: rmin_expected = 150
  real, parameter :: rcut_expected = 1000
  real, parameter :: wilt_expected = 0.08
  real, parameter :: vegheight_expected = 0.5
  real, parameter :: width_expected = 0.12
  character, parameter :: steady_expected = 'Y'
  real, parameter :: ci_expected = 300
  real, parameter :: co_expected = 330
  real, parameter :: coz_sfc_expected = 0
  real, parameter :: coz_air_expected = 0.08
!  real, parameter :: NOBS_pTq_expected = 12
!  real, parameter :: NOBS_wind_expected = 11
!  real, parameter :: station_height_expected = 0.886
!  real, parameter :: ugs_expected = 4.649
!  real, parameter :: vgs_expected = 8.445
!  real, parameter :: ps_expected(12) = (/967,24,5,949,24,5,900,23,3,850,24,8 /)
!  real, parameter :: ts_expected(12) = (/700,11,7,677,9,11,628,4,6,612,2,7/)
!  real, parameter :: dep_expected(12) = (/606,2,15,530,-7,17,460,-10,30,250,-40,30/)
!  real, parameter :: dd0_expected(11) = (/180,7,0,185,10,1,225,35,3,240,25/)
!  real, parameter :: ff0_expected(11) = (/5,225,15,7,215,15,9,230,30,14,240/)
!  real, parameter :: zh_expected(11) = (/25,20,245,44,30,255,43,41,195,14,54/)
  real, parameter :: ozone_flux_plant_expected = 581.332764
  real, parameter :: ozone_fglobal_expected = 10028.6299

! ieee_arithmetic setup -- choices: up, down, nearest, to zero
  if (ieee_support_rounding(IEEE_NEAREST)) then
    call ieee_set_rounding_mode(IEEE_NEAREST)
  end if

! mod_testing variable setup
  n = 1
  ntests = 48
  call initialize_tests(tests,ntests)
! end  mod_testing variable setup

!
! splint_test
!
  if (splint_test) then
    do i = 1,splint_max_array
      splint_arg1(i) = 2*i
      splint_arg2(i) = 2*i
    end do
    splint_arg4 = 12
    if (.not. allocated(splint_arg3)) then
      allocate(splint_arg3(splint_arg4))
      do i = 1,splint_arg4
        splint_arg3(i) = 2*i
      end do
    end if
    splint_arg5 = 3
    splint_output=splint(splint_arg1, splint_arg2, splint_arg3, splint_arg4, splint_arg5)
    tests(n) = assert(eq(splint_output,splint_expected), 'SPLINT')
    n = n + 1
    if (allocated(splint_arg3)) then
      deallocate(splint_arg3)
    end if
  end if

!
! spline_test
!
  if (spline_test) then
    do i = 1,62
      spline_arg1(i) = 2*i
      spline_arg2(i) = 2*i
    end do
    spline_arg3 = 62
    spline_arg4 = 1.0
    spline_arg5 = 2.5
    if (.not. allocated(spline_output)) then
      allocate(spline_output(spline_arg3))
    end if
    tests(n) = assert(eq(                                                   &
       spline(spline_arg1,spline_arg2,spline_arg3,spline_arg4,spline_arg5), &
       spline_expected) .eqv. .true., 'SPLINE')
    n = n + 1
    if (allocated(spline_output)) then
      deallocate(spline_output)
    end if
  end if

!
! transm_ftabs_test
!
! Test the "else" branch in ftabsT()

  if (transm_ftabs_test) then
    call ftabsT_init
    transm_ftabs_arg1 = 2.75401473
    tests(n) = assert(eq(ftabst(transm_ftabs_arg1),transm_ftabs_expected), 'ftabs() else')
    n = n + 1
  end if

!
! transm_ftabs_test2
!
! Test the "if" branch in ftabsT()

  if (transm_ftabs_test2) then
    call ftabsT_init
    transm_ftabs_arg1 = 11.0
    tests(n) = assert(eq(ftabst(transm_ftabs_arg1),transm_ftabs_expected2), 'ftabs() if')
    n = n + 1
  end if

!
! transm_ftscatT_test
!
! Test the "else" branch in ftscatT()

  if (transm_ftscatT_test) then
    call ftscatT_init
    transm_ftscatT_arg1 = 2.75401473
    tests(n) = assert(eq(ftscatT(transm_ftscatT_arg1),transm_ftscatT_expected), 'ftscatT() else')
    n = n + 1
  end if

!
! transm_ftscatT_test2
!
! Test the "if" branch in ftscatT()

  if (transm_ftscatT_test2) then
    call ftscatT_init
    transm_ftscatT_arg1 = 11.0
    tests(n) = assert(eq(ftscatT(transm_ftscatT_arg1),transm_ftscatT_expected2), 'ftscatT() if')
    n = n + 1
  end if

!
! transm_fbscatT_test
!
! Test the "else" branch in fbscatT()

  if (transm_fbscatT_test) then
    call fbscatT_init
    transm_fbscatT_arg1 = 2.75401473
    tests(n) = assert(eq(fbscatT(transm_fbscatT_arg1),transm_fbscatT_expected), 'fbscatT() else')
    n = n + 1
  end if

!
! transm_fbscatT_test2
!
! Test the "if" branch in fbscatT()

  if (transm_fbscatT_test2) then
    call fbscatT_init
    transm_fbscatT_arg1 = 11.0
    tests(n) = assert(eq(fbscatT(transm_fbscatT_arg1),transm_fbscatT_expected2), 'fbscatT() if')
    n = n + 1
  end if

!
! advect_func_test
!
  if (advect_func_test) then
    call advect_init
    tests(n) = assert(eq(advect(),advect_test_expected), 'advect()')
    n = n + 1
  end if

!
! cond_test
!
  if (cond_test) then
    call cond_init
    tests(n) = assert(eq(cond(),cond_test_expected),'cond()')
    n = n + 1
  end if

!
! stomfs_test
!
  if (stomfs_test) then
    call stomfs_init
    tests(n) = assert(eq(stomfs(sc, sol),stomfs_test_expected), 'stomfs()')
    n = n + 1
  end if

!
! stomrs_test_hi_temp
!
  if (stomrs_test_hi_temp) then
    call stomrs_init
    TF = MAXTEMP + 1.0
    tests(n) = assert(eq( &
               stomrs(ft,tf,rmin,mintemp,maxtemp,psisup,psiwc,b1,b2,psie,psice,fs), &
               stomrs_test_hi_temp_expected), 'stomrs_test_hi_temp_expected')
    n = n + 1
  end if

!
! stomrs_test_lo_temp
!
  if (stomrs_test_lo_temp) then
    call stomrs_init
    TF = MINTEMP - 1.0
    tests(n) = assert(eq( &
               stomrs(ft,tf,rmin,mintemp,maxtemp,psisup,psiwc,b1,b2,psie,psice,fs), &
               stomrs_test_lo_temp_expected), 'stomrs_test_lo_temp_expected')
    n = n + 1
  end if

!
! stomrs_test_hi_psi
!
  if (stomrs_test_hi_psi) then
    call stomrs_init
    TF = (MINTEMP + MAXTEMP)/2
    PSISUP = PSIWC + 1.0
    tests(n) = assert(eq( &
               stomrs(ft,tf,rmin,mintemp,maxtemp,psisup,psiwc,b1,b2,psie,psice,fs), &
               stomrs_test_hi_psi_expected), 'stomrs_test_hi_psi')
    n = n + 1
  end if

!
! stomrs_test_lo_psi
!
  if (stomrs_test_lo_psi) then
    call stomrs_init
    TF = (MINTEMP + MAXTEMP)/2
    PSISUP = PSIWC - 1.0
    tests(n) = assert(eq( &
               stomrs(ft,tf,rmin,mintemp,maxtemp,psisup,psiwc,b1,b2,psie,psice,fs), &
               stomrs_test_lo_psi_expected), 'stomrs_test_lo_psi')
    n = n + 1
  end if

!
! stomc_test
!

  if (stomc_test) then
    call stomc_init
    tests(n) = assert(eq(stomc(ft,rmin,b1,psice,fs),stomc_test_expected),'stomc_test')
    n = n + 1
  end if

!
! you_star_test
!

  if (you_star_test) then
    call you_star_init
    tests(n) = assert(eq( &
               you_star(you_star_arg1,you_star_arg2,you_star_arg3,you_star_arg4,you_star_arg5), &
               you_star_test_expected),'you_star_test')
    n = n + 1
  end if

!
! r_ohms_test
!

  if (r_ohms_test) then
    call r_ohms_init
    tests(n) = assert(eq( &
               r_ohms(r_ohms_arg1,r_ohms_arg2,r_ohms_arg3,r_ohms_arg4,r_ohms_arg5), &
               r_ohms_test_expected),'r_ohms_test')
    n = n + 1
  end if

!
! windf_test
!

  if (windf_test) then
    call windf_init
    tests(n) = assert(eq(windf(windf_arg1,windf_arg2,windf_arg3,windf_arg4,windf_arg5),windf_test_expected),'windf_test')
    n = n + 1
  end if

!
! stab_test
!

  if (stab_test) then
    call stab_init
    tests(n) = assert(eq(stab(stab_arg1,stab_arg2),stab_test_expected),'stab_test')
    n = n + 1
  end if

!
! stabh_test
!

  if (stabh_test) then
    ! stabh takes same argument as stab so we'll recycle initialization
    call stab_init
    tests(n) = assert(eq(stabh(stab_arg1,stab_arg2),stabh_test_expected),'stabh_test')
    n = n + 1
  end if

!
! fstabh_test
!

  if (fstabh_test) then
    call fstabh_init
    tests(n) = assert(eq(fstabh(fstabh_arg1,fstabh_arg2),fstabh_test_expected),'fstabh_test')
    n = n + 1
  end if

!
! fstabm_test
!

  if (fstabm_test) then
    ! fstabm takes the same arguments as fstabh so we'll recycle the initialization
    call fstabh_init
    tests(n) = assert(eq(fstabm(fstabh_arg1,fstabh_arg2),fstabm_test_expected,eq_tol),'fstabm_test')
    n = n + 1
  end if

!
! restrn_test
!

  if (restrn_test) then
    call restrn_init
    tests(n) = assert(eq(restrn(restrn_arg1,restrn_arg2,restrn_arg3,restrn_arg4),restrn_test_expected),'restrn()')
    n = n + 1
  end if

!
! psgcal_test
!

  if (psgcal_test) then
    call psgcal_init
    psig = psgcal(thmax,thv,cosbyb,psis)
    tests(n) = assert(eq(PSIG,psgcal_test_expected), 'psgcal_test')
    n = n + 1
  end if

!
! vegflx_test
!

  if (vegflx_test) then
    call vegflx_init
    call vegflx(vegflx_arg1)
    tests(n) = assert(eq(QAF,vegflx_test_expected), 'vegflx_test')
    n = n + 1
  end if

!
! avr_test_init1
!

  if (avr_test_init1) then
    call avr_init
    do i=1,10
    call avr(avr_arg1,avr_arg2,avr_arg3)
    avr_arg1 = avr_arg1+1
    end do
!    write(6,*) 'avr after 10 passes: avr_arg1, avr_arg2, avr_arg3: ',avr_arg1, &
!      avr_arg2, avr_arg3
    tests(n) = assert(eq(avr_arg2,avr_test_init1_expected), 'avr_test_init1')
    n = n + 1
  end if

!
! avr_test_init_not1
!

!  if (avr_test_init_not1) then
!    call avr_init
!    call avr(avr_arg1,avr_arg2,avr_arg3)
!    call avr(avr_arg1,avr_arg2,avr_arg3)
!    call avr(avr_arg1,avr_arg2,avr_arg3)
!    call avr(avr_arg1,avr_arg2,avr_arg3)
!    tests(n) = assert(eq(avr_arg2,avr_test_init_not1_expected), 'avr_test_init_not1')
!    n = n + 1
!  end if

!
! ozone_test
!

  if (ozone_test) then
    call ozone_init
    call ozone
    tests(n) = assert(eq(flux_plant,ozone_flux_plant_expected), 'ozone flux_plant')
    n = n + 1
    tests(n) = assert(eq(fglobal,ozone_fglobal_expected), 'ozone fglobal')
    n = n + 1
  end if

!
! co2flx_test
!

  if (co2flx_test) then
    call co2flx_init
    call co2flx
    tests(n) = assert(eq(ccan,co2flx_ccan_expected), 'co2flx ccan')
    n = n + 1
    tests(n) = assert(eq(fco2,co2flx_fco2_expected), 'co2flx fco2')
    n = n + 1
  end if

!
! veghot_test
!

  if (veghot_test) then
    call veghot_init
    call veghot(veghot_arg1, veghot_arg2)
    tests(n) = assert(eq(veghot_arg2,veghot_heatv_expected), 'veghot heatv')
    n = n + 1
  end if

!
! hot_test
!

  if (hot_test) then
    call hot_init
    ! Case I: FRVEG == 1 AND RNET > 0
    FRVEG = 1.0
    RNET = 1.0
    call hot(hot_arg1,hot_arg2,hot_arg3,hot_arg4)
    tests(n) = assert(eq(heat, hot_caseI_expected), 'hot case I')
    n = n + 1
    ! Case II: 1 > FRVEG > 0 AND RNET > 0
    FRVEG = 0.5
    RNET = 1.0
    call hot(hot_arg1,hot_arg2,hot_arg3,hot_arg4)
    tests(n) = assert(eq(heat, hot_caseII_expected), 'hot case II')
    n = n + 1
    ! Case III: All others
    FRVEG = 1.5
    RNET = 1.0
    call hot(hot_arg1,hot_arg2,hot_arg3,hot_arg4)
    tests(n) = assert(eq(heat, hot_caseIII_expected), 'hot case III')
    n = n + 1
  end if

!
! albedo_test (4 cases)
!

  if (albedo_test) then
    call albedo_init

!  Case I
    ALBG = 0.0
    ALBF = 0.0
    call albedo(albedo_test_arg1)
    tests(n) = assert(eq(ALBDOE,albedo_caseI_expected), 'albedo case I')
    n = n + 1

!  Case II
    ALBG = 20.0
    ALBF = 0.0
    call albedo(albedo_test_arg1)
    tests(n) = assert(eq(ALBDOE,albedo_caseII_expected), 'albedo case II')
    n = n + 1

!  Case III
    ALBG = 0.0
    ALBF = 20.0
    call albedo(albedo_test_arg1)
    tests(n) = assert(eq(ALBDOE,albedo_caseIII_expected), 'albedo case III')
    n = n + 1

!  Case IV
    ALBG = 46.0
    ALBF = 71.0
    call albedo(albedo_test_arg1)
    tests(n) = assert(eq(ALBDOE,albedo_caseIV_expected), 'albedo case IV')
    n = n + 1
  end if

!
! intpol_test (4 array variables)
!
  if (intpol_test) then
    call intpol_init
    call intpol
    tests(n) = assert(eq(ug,intpol_ug_exp) .eqv. .true., 'Intpol UG(50)')
    n = n + 1
    tests(n) = assert(eq(vg,intpol_vg_exp) .eqv. .true., 'Intpol VG(50)')
    n = n + 1
    tests(n) = assert(eq(ugd,intpol_ugd_exp) .eqv. .true., 'Intpol UGD(50)')
    n = n + 1
    tests(n) = assert(eq(vgd,intpol_vgd_exp) .eqv. .true., 'Intpol VGD(50)')
    n = n + 1
  end if

!
! prfile_test (4 array variables and scalar variable)
!

  if (prfile_test) then
    call prfile_init
    call prfile
    tests(n) = assert(eq(u_fine,prfile_uf_exp) .eqv. .true., 'PRFILE u_fine')
    n = n + 1
    tests(n) = assert(eq(v_fine,prfile_vf_exp) .eqv. .true., 'PRFILE v_fine')
    n = n + 1
    tests(n) = assert(eq(t_fine,prfile_tf_exp) .eqv. .true., 'PRFILE t_fine')
    n = n + 1
    tests(n) = assert(eq(q_fine,prfile_qf_exp) .eqv. .true., 'PRFILE q_fine')
    n = n + 1
    tests(n) = assert(eq(awind,prfile_awind_exp), 'PRFILE awind')
    n = n + 1
  end if

!
! f_control_test
!

  if (f_control_test) then
!    call start
    tests(n) = assert(iyr == iyr_expected, 'iyr')
    n = n + 1
    tests(n) = assert(imo == imo_expected, 'imo')
    n = n + 1
    tests(n) = assert(iday == iday_expected, 'iday')
    n = n + 1
    tests(n) = assert(eq(tz,tz_expected), 'tz')
    n = n + 1
    tests(n) = assert(eq(xlat,xlat_expected), 'xlat')
    n = n + 1
    tests(n) = assert(eq(xlong,xlong_expected), 'xlong')
    n = n + 1
    tests(n) = assert(eq(strtim,strtim_expected), 'strtim')
    n = n + 1
    tests(n) = assert(eq(timend,timend_expected), 'timend')
    n = n + 1
    tests(n) = assert(eq(outtt,outtt_expected), 'outtt')
    n = n + 1
    tests(n) = assert(eq(slope,slope_expected), 'slope')
    n = n + 1
    tests(n) = assert(eq(aspect,aspect_expected), 'aspect')
    n = n + 1
    tests(n) = assert(eq(f,f_expected), 'f')
    n = n + 1
    tests(n) = assert(eq(fsub,fsub_expected), 'f_sub')
    n = n + 1
    tests(n) = assert(eq(wmax,wmax_expected), 'wmax')
    n = n + 1
    tests(n) = assert(eq(btemp,btemp_expected), 'btemp')
    n = n + 1
    tests(n) = assert(eq(tp,tp_expected), 'tp')
    n = n + 1
    tests(n) = assert(dual_ti == dual_ti_expected, 'dual_ti')
    n = n + 1
    tests(n) = assert(eq(ti_a,ti_a_expected), 'ti_a')
    n = n + 1
    tests(n) = assert(eq(ti_b,ti_b_expected), 'ti_b')
    n = n + 1
!    tests(n) = assert(albedo_gflag == albedo_gflag_expected, 'albedo_gflag')
!    n = n + 1
    tests(n) = assert(eq(albg,albg_expected), 'albg')
    n = n + 1
    tests(n) = assert(eq(epsi,epsi_expected), 'epsi')
    n = n + 1
!    tests(n) = assert(index_soils == index_soils_expected, 'index_soils')
!    n = n + 1
    tests(n) = assert(eq(omega,omega_expected), 'omega')
    n = n + 1
    tests(n) = assert(eq(zo,zo_expected), 'zo')
    n = n + 1
!    tests(n) = assert(obst_hgt == obst_hgt_expected, 'obst_hgt')
!    n = n + 1
!    tests(n) = assert(cloud_flag == cloud_flag_expected, 'cloud_flag')
!    n = n + 1
    tests(n) = assert(cld_fract == cld_fract_expected, 'cld_fract')
    n = n + 1
    tests(n) = assert(eq(frveg,frveg_expected), 'frveg')
    n = n + 1
    tests(n) = assert(eq(xlai,xlai_expected), 'xlai')
    n = n + 1
    tests(n) = assert(eq(epsf,epsf_expected), 'epsf')
    n = n + 1
!    tests(n) = assert(albedo_fflag == albedo_fflag_expected, 'albedo_flag')
!    n = n + 1
    tests(n) = assert(eq(albf,albf_expected), 'albf')
    n = n + 1
    tests(n) = assert(stmtype == stmtype_expected, 'stmtype')
    n = n + 1
!    tests(n) = assert(index_veggies == index_veggies_expected, 'index_veggies')
!    n = n + 1
    tests(n) = assert(eq(volrel,volrel_expected), 'volrel')
    n = n + 1
    tests(n) = assert(eq(rmin,rmin_expected), 'rmin')
    n = n + 1
    tests(n) = assert(eq(rcut,rcut_expected), 'rcut')
    n = n + 1
    tests(n) = assert(eq(wilt,wilt_expected), 'wilt')
    n = n + 1
    tests(n) = assert(eq(vegheight,vegheight_expected), 'vegheight')
    n = n + 1
    tests(n) = assert(eq(width,width_expected), 'width')
    n = n + 1
    tests(n) = assert(steady == steady_expected, 'steady')
    n = n + 1
    tests(n) = assert(eq(ci,ci_expected), 'ci')
    n = n + 1
    tests(n) = assert(eq(co,co_expected), 'co')
    n = n + 1
    tests(n) = assert(eq(coz_sfc,coz_sfc_expected), 'coz_sfc')
    n = n + 1
    tests(n) = assert(eq(coz_air,coz_air_expected), 'coz_air')
    n = n + 1
!    call snding
!    tests(n) = assert(NOBS_pTq == NOBS_pTq_expected, 'NOBS_pTq')
!    n = n + 1
!    tests(n) = assert(NOBS_wind == NOBS_wind_expected, 'NOBS_wind')
!    n = n + 1
!    tests(n) = assert(station_height == station_height_expected, 'station_height')
!    n = n + 1
!    tests(n) = assert(ugs == ugs_expected, 'ugs')
!    n = n + 1
!    tests(n) = assert(vgs == vgs_expected, 'vgs')
!    n = n + 1
!    tests(n) = assert(ps == ps_expected, 'ps')
!    n = n + 1
!    tests(n) = assert(ts == ts_expected, 'ts')
!    n = n + 1
!    tests(n) = assert(dep == dep_expected, 'dep')
!    n = n + 1
!    tests(n) = assert(dd0 == dd0_expected, 'dd0')
!    n = n + 1
!    tests(n) = assert(ff0 == ff0_expected, 'ff0')
!    n = n + 1
!    tests(n) = assert(zh == zh_expected, 'zh')
!    n = n + 1
  end if

! Report test summary
  test_failed = .false.
  call report_tests(tests,test_failed)
  if (test_failed) stop 1

contains
  subroutine ftabsT_init
    OMEGA = 3.13
    ABSTBL(9) = 0.718309104
    ABSTBL(10) = 0.707092881
    ABSTBL(46) = 0.437107921
    PS1 = 967.0
    return
  end subroutine ftabsT_init

  subroutine ftscatT_init
    OMEGA = 3.13
    SCATBL(9) = 0.548576415
    SCATBL(10) = 0.527300596
    SCATBL(46) = 0.173223004
    PS1 = 967.0
    return
  end subroutine ftscatT_init

  subroutine fbscatT_init
    OMEGA = 3.13
    BSCTBL(9) = 0.307897508
    BSCTBL(10) = 0.307446688
    BSCTBL(46) = 0.304252803
    PS1 = 967.0
    return
  end subroutine fbscatT_init

  subroutine advect_init
    CF = 9.19953891E-05
    OTEMP = 295.149994
    VGD(5) = 8.41951942
    VGD(3) = 8.43165302
    VGD(1) = 8.44378662
    UGD(5) = 13.0658665
    UGD(3) = 9.05783463
    UGD(1) = 5.04980326
    return
  end subroutine advect_init

  subroutine cond_init
    RKS = 6.60699987
    THV = 1.00000000
    THMAX = 0.338999987
    COSBYB = 2.78999996
    return
  end subroutine cond_init

  subroutine stomfs_init
    sc = 1.0
    sol = 1.0
    return
  end subroutine stomfs_init

  subroutine stomrs_init
    RMIN = 1.0
    MINTEMP = 1.0
    MAXTEMP = 100.0
    PSIWC = 1.0
    PSIE = 1.0
    PSICE = 2.0
    B1 = 1.0
    B2 = 2.0
    call stomfs_init
    FS = stomfs(sc,sol)
    return
  end subroutine stomrs_init

  subroutine stomc_init
    B1 = 1.0
    PSICE = 2.0
    RMIN = 1.0
    call stomfs_init
    FS = stomfs(sc,sol)
    return
  end subroutine stomc_init

  subroutine you_star_init
    you_star_arg1 = 0.1
    you_star_arg2 = 1.0
    you_star_arg3 = 2.0
    you_star_arg4 = 3.0
    you_star_arg5 = karman
    return
  end subroutine you_star_init

  subroutine r_ohms_init
    r_ohms_arg1 = 0.1
    r_ohms_arg2 = 1.0
    r_ohms_arg3 = 2.0
    r_ohms_arg4 = 3.0
    r_ohms_arg5 = karman
    return
  end subroutine r_ohms_init

  subroutine windf_init
    windf_arg1 = 1.0
    windf_arg2 = 2.0
    windf_arg3 = 3.0
    windf_arg4 = 4.0
    windf_arg5 = karman
    return
  end subroutine windf_init

  subroutine stab_init
    ! also used to test stabh
    stab_arg1 = 0.1
    stab_arg2 = 2.0
    return
  end subroutine stab_init

  subroutine fstabh_init
    ! also used to test fstabm
    fstabh_arg1 = 1.0
    fstabh_arg2 = 2.0
    return
  end subroutine fstabh_init

  subroutine restrn_init
    restrn_arg1 = 1.0
    restrn_arg2 = 2.0
    restrn_arg3 = 3.0
    restrn_arg4 = karman
    return
  end subroutine

  subroutine psgcal_init
    THMAX = 0.338999987
    THV = 1.0
    COSBYB = 2.78999996
    PSIS = 1.0
    return
  end subroutine psgcal_init

  subroutine vegflx_init
    vegflx_arg1 = 0.3
    TAF = 10.0
    CHF = 10.0
    F = 100.0
    CHG = 0.1
    QAF = 0.5
    HG = 1.0
    CHG = 100.5
    XLEF = 1.2
    CHA = 100.0
    TA = 0.4
    CHF = 10.0
    TF = 1.0
    RST = 1.0
    QA = 100.0
    QSTF = 10.0
    return
  end subroutine vegflx_init

  subroutine avr_init
    avr_arg1 = 293.0
    avr_arg3 = 1
    return
  end subroutine avr_init

  subroutine ozone_init
    USTAR = 1.0
    UTEN = 1.0
    UAF = 0.5
    XLAI = 2.0
    RAF = 1.0
    RCUT = 1.0
    RST = 1.0
    CHG = 2.0
    RZASCR = 1.0
    COZ_AIR = 1.0
    frveg = 0.5
    sumo3 = 0.1
    return
  end subroutine ozone_init

  subroutine co2flx_init
    USTAR = 1.0
    UTEN = 2.0
    UAF = 1.0
    XLAI = 1.0
    RZASCR = 2.0
    RAF = 1.0
    RST = 1.0
    CO = 2.0
    CI = 1.0
    FRVEG = 0.5
    return
  end subroutine co2flx_init

  subroutine veghot_init
    TF = 1.0
    TAF = 1.0
    CHF = 1.0
    CHG = 1.0
    HF = 1.0
    TT(2) = 1.0
    Z(2) = 1.0
    LAMBDA = 1.0
    RNETG = 3.0
    XLEG = 1.0
    veghot_arg1 = 1.0
    veghot_arg2 = 1.0
  end subroutine veghot_init

  subroutine hot_init
    hot_arg2 = 4.0
    hot_arg3 = 2.0
    hot_arg4 = 1.0
    LAMBDA = 1.0
    ATEMP = 2.0
    TT(2) = 1.0
    Z(2) = 1.0
    TF = 1.0
    TAF = 1.0
    CHF = 1.0
    CHG = 1.0
    HF = 1.0
    RNETG = 3.0
    XLEG = 1.0
    GBL_sum = 1.0
    return
  end subroutine hot_init

  subroutine albedo_init
    albedo_test_arg1 = -0.510662317
    XLAI = 1.0
    FRVEG = 1.0
    WGG = 100.0
    WMAX = 1230.0
    return
  end subroutine albedo_init

  subroutine intpol_init
    NTRP = 10
    VGS = 10.0
    UGS = 10.0
    VD = 1.0
    UD = 1.0
    VGD(NTRP) = 1.0
    UGD(NTRP) = 1.0
    return
  end subroutine intpol_init

  subroutine prfile_init
    integer :: i

    call advect_init
    call intpol_init

    ud = 0.0
    vd = 0.0
    ud(1) = 1.0
    vd(1) = 1.0
    ud(3) = 1.0
    vd(3) = 1.0
    do i = 5, NTRP
      ud(i) = 1.0
      vd(i) = 1.0
    end do
    td(1) = 1.0
    qd(1) = 1.0
    return
  end subroutine prfile_init

end program test_simsphere


!  Copyright 2020 University of Warwick

!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at

!      http://www.apache.org/licenses/LICENSE-2.0

!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.

!******************************************************************************
! This module contains the boundary conditions for the entire code
! Any new boundary conditions should be added here
!******************************************************************************

MODULE boundary

  USE shared_data
  USE mpiboundary

  IMPLICIT NONE

CONTAINS

  !****************************************************************************
  ! Set up any necessary variables for the chosen boundary conditions
  !****************************************************************************

  SUBROUTINE set_boundary_conditions

    any_open = .FALSE.
    IF (xbc_min == BC_OPEN .OR. xbc_max == BC_OPEN &
        .OR. ybc_min == BC_OPEN .OR. ybc_max == BC_OPEN &
        .OR. zbc_min == BC_OPEN .OR. zbc_max == BC_OPEN) any_open = .TRUE.

  END SUBROUTINE set_boundary_conditions

  !****************************************************************************
  ! Call all of the boundaries needed by the core Lagrangian solver
  !****************************************************************************

  SUBROUTINE boundary_conditions

    CALL bfield_bcs
    CALL energy_bcs
    CALL density_bcs
    CALL velocity_bcs

  END SUBROUTINE boundary_conditions

  !****************************************************************************
  ! Boundary conditions for magnetic field through plane
  !****************************************************************************

  SUBROUTINE bfield_bcs

    CALL bfield_mpi

    IF (proc_x_min == MPI_PROC_NULL .AND. xbc_min == BC_USER) THEN
      bx(-1,:,:) = bx(1,:,:)
      bx(-2,:,:) = bx(2,:,:)
      by( 0,:,:) = by(1,:,:)
      by(-1,:,:) = by(2,:,:)
      bz( 0,:,:) = bz(1,:,:)
      bz(-1,:,:) = bz(2,:,:)
    END IF

    IF (proc_x_max == MPI_PROC_NULL .AND. xbc_max == BC_USER) THEN
      bx(nx+1,:,:) = bx(nx-1,:,:)
      bx(nx+2,:,:) = bx(nx-2,:,:)
      by(nx+1,:,:) = by(nx  ,:,:)
      by(nx+2,:,:) = by(nx-1,:,:)
      bz(nx+1,:,:) = bz(nx  ,:,:)
      bz(nx+2,:,:) = bz(nx-1,:,:)
    END IF

    IF (proc_y_min == MPI_PROC_NULL .AND. ybc_min == BC_USER) THEN
      bx(:, 0,:) = bx(:,1,:)
      bx(:,-1,:) = bx(:,2,:)
      by(:,-1,:) = by(:,1,:)
      by(:,-2,:) = by(:,2,:)
      bz(:, 0,:) = bz(:,1,:)
      bz(:,-1,:) = bz(:,2,:)
    END IF

    IF (proc_y_max == MPI_PROC_NULL .AND. ybc_max == BC_USER) THEN
      bx(:,ny+1,:) = bx(:,ny  ,:)
      bx(:,ny+2,:) = bx(:,ny-1,:)
      by(:,ny+1,:) = by(:,ny-1,:)
      by(:,ny+2,:) = by(:,ny-2,:)
      bz(:,ny+1,:) = bz(:,ny  ,:)
      bz(:,ny+2,:) = bz(:,ny-1,:)
    END IF

    IF (proc_z_min == MPI_PROC_NULL .AND. zbc_min == BC_USER) THEN
      bx(:,:, 0) = bx(:,:,1)
      bx(:,:,-1) = bx(:,:,2)
      by(:,:, 0) = by(:,:,1)
      by(:,:,-1) = by(:,:,2)
      bz(:,:,-1) = bz(:,:,1)
      bz(:,:,-2) = bz(:,:,2)
    END IF

    IF (proc_z_max == MPI_PROC_NULL .AND. zbc_max == BC_USER) THEN
      bx(:,:,nz+1) = bx(:,:,nz  )
      bx(:,:,nz+2) = bx(:,:,nz-1)
      by(:,:,nz+1) = by(:,:,nz  )
      by(:,:,nz+2) = by(:,:,nz-1)
      bz(:,:,nz+1) = bz(:,:,nz-1)
      bz(:,:,nz+2) = bz(:,:,nz-2)
    END IF

  END SUBROUTINE bfield_bcs

  !****************************************************************************
  ! Boundary conditions for specific internal energy
  !****************************************************************************

  SUBROUTINE energy_bcs

    CALL energy_mpi

    IF (proc_x_min == MPI_PROC_NULL .AND. xbc_min == BC_USER) THEN
      energy( 0,:,:) = energy(1,:,:)
      energy(-1,:,:) = energy(2,:,:)
    END IF

    IF (proc_x_max == MPI_PROC_NULL .AND. xbc_max == BC_USER) THEN
      energy(nx+1,:,:) = energy(nx  ,:,:)
      energy(nx+2,:,:) = energy(nx-1,:,:)
    END IF

    IF (proc_y_min == MPI_PROC_NULL .AND. ybc_min == BC_USER) THEN
      energy(:, 0,:) = energy(:,1,:)
      energy(:,-1,:) = energy(:,2,:)
    END IF

    IF (proc_y_max == MPI_PROC_NULL .AND. ybc_max == BC_USER) THEN
      energy(:,ny+1,:) = energy(:,ny  ,:)
      energy(:,ny+2,:) = energy(:,ny-1,:)
    END IF

    IF (proc_z_min == MPI_PROC_NULL .AND. zbc_min == BC_USER) THEN
      energy(:,:, 0) = energy(:,:,1)
      energy(:,:,-1) = energy(:,:,2)
    END IF

    IF (proc_z_max == MPI_PROC_NULL .AND. zbc_max == BC_USER) THEN
      energy(:,:,nz+1) = energy(:,:,nz  )
      energy(:,:,nz+2) = energy(:,:,nz-1)
    END IF

  END SUBROUTINE energy_bcs

  !****************************************************************************
  ! Boundary conditions for density
  !****************************************************************************

  SUBROUTINE density_bcs

    CALL density_mpi

    IF (proc_x_min == MPI_PROC_NULL .AND. xbc_min == BC_USER) THEN
      rho( 0,:,:) = rho(1,:,:)
      rho(-1,:,:) = rho(2,:,:)
    END IF

    IF (proc_x_max == MPI_PROC_NULL .AND. xbc_max == BC_USER) THEN
      rho(nx+1,:,:) = rho(nx  ,:,:)
      rho(nx+2,:,:) = rho(nx-1,:,:)
    END IF

    IF (proc_y_min == MPI_PROC_NULL .AND. ybc_min == BC_USER) THEN
      rho(:, 0,:) = rho(:,1,:)
      rho(:,-1,:) = rho(:,2,:)
    END IF

    IF (proc_y_max == MPI_PROC_NULL .AND. ybc_max == BC_USER) THEN
      rho(:,ny+1,:) = rho(:,ny  ,:)
      rho(:,ny+2,:) = rho(:,ny-1,:)
    END IF

    IF (proc_z_min == MPI_PROC_NULL .AND. zbc_min == BC_USER) THEN
      rho(:,:, 0) = rho(:,:,1)
      rho(:,:,-1) = rho(:,:,2)
    END IF

    IF (proc_z_max == MPI_PROC_NULL .AND. zbc_max == BC_USER) THEN
      rho(:,:,nz+1) = rho(:,:,nz  )
      rho(:,:,nz+2) = rho(:,:,nz-1)
    END IF

  END SUBROUTINE density_bcs

  !****************************************************************************
  ! Boundary conditions for temperature
  !****************************************************************************

  SUBROUTINE temperature_bcs

    CALL density_mpi

    IF (proc_x_min == MPI_PROC_NULL .AND. xbc_min == BC_USER) THEN
      temperature( 0,:,:) = temperature(1,:,:)
      temperature(-1,:,:) = temperature(2,:,:)
    END IF

    IF (proc_x_max == MPI_PROC_NULL .AND. xbc_max == BC_USER) THEN
      temperature(nx+1,:,:) = temperature(nx  ,:,:)
      temperature(nx+2,:,:) = temperature(nx-1,:,:)
    END IF

    IF (proc_y_min == MPI_PROC_NULL .AND. ybc_min == BC_USER) THEN
      temperature(:, 0,:) = temperature(:,1,:)
      temperature(:,-1,:) = temperature(:,2,:)
    END IF

    IF (proc_y_max == MPI_PROC_NULL .AND. ybc_max == BC_USER) THEN
      temperature(:,ny+1,:) = temperature(:,ny  ,:)
      temperature(:,ny+2,:) = temperature(:,ny-1,:)
    END IF

    IF (proc_z_min == MPI_PROC_NULL .AND. zbc_min == BC_USER) THEN
      temperature(:,:, 0) = temperature(:,:,1)
      temperature(:,:,-1) = temperature(:,:,2)
    END IF

    IF (proc_z_max == MPI_PROC_NULL .AND. zbc_max == BC_USER) THEN
      temperature(:,:,nz+1) = temperature(:,:,nz  )
      temperature(:,:,nz+2) = temperature(:,:,nz-1)
    END IF

  END SUBROUTINE temperature_bcs

  !****************************************************************************
  ! Full timestep velocity boundary conditions
  !****************************************************************************

  SUBROUTINE velocity_bcs

    CALL velocity_mpi

    IF (proc_x_min == MPI_PROC_NULL .AND. xbc_min == BC_USER) THEN
      vx(-2:0,:,:) = 0.0_num
      vy(-2:0,:,:) = 0.0_num
      vz(-2:0,:,:) = 0.0_num
    END IF

    IF (proc_x_max == MPI_PROC_NULL .AND. xbc_max == BC_USER) THEN
      vx(nx:nx+2,:,:) = 0.0_num
      vy(nx:nx+2,:,:) = 0.0_num
      vz(nx:nx+2,:,:) = 0.0_num
    END IF

    IF (proc_y_min == MPI_PROC_NULL .AND. ybc_min == BC_USER) THEN
      vx(:,-2:0,:) = 0.0_num
      vy(:,-2:0,:) = 0.0_num
      vz(:,-2:0,:) = 0.0_num
    END IF

    IF (proc_y_max == MPI_PROC_NULL .AND. ybc_max == BC_USER) THEN
      vx(:,ny:ny+2,:) = 0.0_num
      vy(:,ny:ny+2,:) = 0.0_num
      vz(:,ny:ny+2,:) = 0.0_num
    END IF

    IF (proc_z_min == MPI_PROC_NULL .AND. zbc_min == BC_USER) THEN
      vx(:,:,-2:0) = 0.0_num
      vy(:,:,-2:0) = 0.0_num
      vz(:,:,-2:0) = 0.0_num
    END IF

    IF (proc_z_max == MPI_PROC_NULL .AND. zbc_max == BC_USER) THEN
      vx(:,:,nz:nz+2) = 0.0_num
      vy(:,:,nz:nz+2) = 0.0_num
      vz(:,:,nz:nz+2) = 0.0_num
    END IF

  END SUBROUTINE velocity_bcs

  !****************************************************************************
  ! Half timestep velocity boundary conditions
  !****************************************************************************

  SUBROUTINE remap_v_bcs

    CALL remap_v_mpi

    IF (proc_x_min == MPI_PROC_NULL .AND. xbc_min == BC_USER) THEN
      vx1(-2:0,:,:) = 0.0_num
      vy1(-2:0,:,:) = 0.0_num
      vz1(-2:0,:,:) = 0.0_num
    END IF

    IF (proc_x_max == MPI_PROC_NULL .AND. xbc_max == BC_USER) THEN
      vx1(nx:nx+2,:,:) = 0.0_num
      vy1(nx:nx+2,:,:) = 0.0_num
      vz1(nx:nx+2,:,:) = 0.0_num
    END IF

    IF (proc_y_min == MPI_PROC_NULL .AND. ybc_min == BC_USER) THEN
      vx1(:,-2:0,:) = 0.0_num
      vy1(:,-2:0,:) = 0.0_num
      vz1(:,-2:0,:) = 0.0_num
    END IF

    IF (proc_y_max == MPI_PROC_NULL .AND. ybc_max == BC_USER) THEN
      vx1(:,ny:ny+2,:) = 0.0_num
      vy1(:,ny:ny+2,:) = 0.0_num
      vz1(:,ny:ny+2,:) = 0.0_num
    END IF

    IF (proc_z_min == MPI_PROC_NULL .AND. zbc_min == BC_USER) THEN
      vx1(:,:,-2:0) = 0.0_num
      vy1(:,:,-2:0) = 0.0_num
      vz1(:,:,-2:0) = 0.0_num
    END IF

    IF (proc_z_max == MPI_PROC_NULL .AND. zbc_max == BC_USER) THEN
      vx1(:,:,nz:nz+2) = 0.0_num
      vy1(:,:,nz:nz+2) = 0.0_num
      vz1(:,:,nz:nz+2) = 0.0_num
    END IF

  END SUBROUTINE remap_v_bcs

END MODULE boundary


text:
#include "eiscor.h"
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! z_unifact_2x2diagblock
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! This routine computes a two by two diagonal block of a unitary upper
! hessenberg matrix that is stored as a product of givens rotations
! and a complex diagonal matrix.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! INPUT VARIABLES:
!
!  TOP             LOGICAL
!                    .TRUE.: top block is computed
!                    .FALSE.: bottom block is computed
!
!  Q               REAL(8) array of dimension (6)
!                    array of generators for givens rotations
!                    generators must be orthogonal to working precision
!
!  D               REAL(8) array of dimension (4)
!                    array of generators for complex diagonal matrix
!                    on output contains the eigenvalues
!
! OUTPUT VARIABLES:
!
!  H              COMPLEX(8) array of dimension (2,2)
!                   on exit contains the desired 2x2 block
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine z_unifact_2x2diagblock(TOP,Q,D,H)

  implicit none

  ! input variables
  logical, intent(in) :: TOP
  real(8), intent(in) :: Q(6), D(4)
  complex(8), intent(inout) :: H(2,2)

  ! TOP == .TRUE.
  if (TOP) then

    ! initialize H
    H(1,1) = cmplx(Q(1),Q(2),kind=8)
    H(2,1) = cmplx(Q(3),0d0,kind=8)
    H(1,2) = -H(2,1)
    H(2,2) = conjg(H(1,1))

    ! second rotation
    H(:,2) = H(:,2)*cmplx(Q(4),Q(5),kind=8)

  ! TOP == .FALSE.
  else

    ! initialize H
    H(1,1) = cmplx(Q(4),Q(5),kind=8)
    H(2,1) = cmplx(Q(6),0d0,kind=8)
    H(1,2) = -H(2,1)
    H(2,2) = conjg(H(1,1))

    ! second rotation
    H(1,:) = H(1,:)*cmplx(Q(1),-Q(2),kind=8)

  end if

  ! diagaonal
  H(:,1) = H(:,1)*cmplx(D(1),D(2),kind=8)
  H(:,2) = H(:,2)*cmplx(D(3),D(4),kind=8)

end subroutine z_unifact_2x2diagblock

max_stars_repo_name: trcameron/FPML

max_issues_repo_name: trcameron/FPML

max_forks_repo_name: trcameron/FPML


text:
MODULE environ

	use constantes
	implicit none

	contains

subroutine magneticf(V,B)
	!----------------------------------------------------------
	! Magnetic field computation.
	! CALLING SEQUENCE: call magneticf(V,B)
	! INPUTS:  V: position vector in cartesian coordinates
	! OUTPUTS: B: magnetic field vector in cartesian coordinate (Tesla)
	!----------------------------------------------------------
	real(kind=8), dimension(:,:),intent(in)  :: V
	real(kind=8), dimension(:,:),intent(out) :: B

	B(1,:)=0.
	B(2,:)=0.
	B(3,:)=B0*(V(3,:)/z0B)**(-3)

end subroutine magneticf

subroutine read_environ()
        !----------------------------------------------------
        ! Lit le fichier init_environ.txt et range les donnees
        !  dans un fichier
        !----------------------------------------------------
        use constantes

        open(40,file='init_environ.txt')
        read(40,*)n0
        read(40,*)z0
        read(40,*)B0
        read(40,*)z0B
        close(40)
        write (*,*)"n0 ",n0," z0 ",z0," B0 ",B0," z0B ",z0B

end subroutine read_environ

subroutine density(V,Ne)
	!-----------------------------------------------------
	! Electron density computation
	! CALLING SEQUENCE: call density(V,Ne)
	! INPUTS:   V: position vector in cartesian coordinates
	! OUTPUTS: Ne: electron density (cm-3)
	!-----------------------------------------------------
        use constantes
        real(kind=8),dimension(:,:),intent(in)    :: V
	real(kind=8), dimension(:), intent(out)   :: Ne

!	real(kind=8)                              :: n0,z0

!	open(40, file='init_environ.txt')
!	read(40,*) n0
!	read(40,*) z0

	Ne(:)=n0*(V(3,:)/z0)**(-2)

!	close(40)

end subroutine density

END MODULE environ

max_stars_repo_name: maserlib/ARTEMIS-P

max_issues_repo_name: maserlib/ARTEMIS-P

max_forks_repo_name: maserlib/ARTEMIS-P


text:
      SUBROUTINE STEPDS(F,NEQN,Y,X,H,EPS,WT,START,HOLD,K,KOLD,
     1   CRASH,PHI,P,YP,ALPHA,W,G,KSTEPS,XOLD,IVC,IV,KGI,GI,
     2   FPWA1,FPWA2,FPWA3,IFPC1,IFPWA1,FPWA4,FPWA5,IFPC2,IFPC3,
     3   PAR,IPAR)
C
C
C   WRITTEN BY L. F. SHAMPINE AND M. K. GORDON
C
C   ABSTRACT
C
C   SUBROUTINE  STEPS  IS NORMALLY USED INDIRECTLY THROUGH SUBROUTINE
C   DEABM .  BECAUSE  DEABM  SUFFICES FOR MOST PROBLEMS AND IS MUCH
C   EASIER TO USE, USING IT SHOULD BE CONSIDERED BEFORE USING  STEPS
C   ALONE.
C
C   SUBROUTINE STEPS INTEGRATES A SYSTEM OF  NEQN  FIRST ORDER ORDINARY
C   DIFFERENTIAL EQUATIONS ONE STEP, NORMALLY FROM X TO X+H, USING A
C   MODIFIED DIVIDED DIFFERENCE FORM OF THE ADAMS PECE FORMULAS.  LOCAL
C   EXTRAPOLATION IS USED TO IMPROVE ABSOLUTE STABILITY AND ACCURACY.
C   THE CODE ADJUSTS ITS ORDER AND STEP SIZE TO CONTROL THE LOCAL ERROR
C   PER UNIT STEP IN A GENERALIZED SENSE.  SPECIAL DEVICES ARE INCLUDED
C   TO CONTROL ROUNDOFF ERROR AND TO DETECT WHEN THE USER IS REQUESTING
C   TOO MUCH ACCURACY.
C
C   THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT,
C   COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS, THE INITIAL
C   VALUE PROBLEM  BY L. F. SHAMPINE AND M. K. GORDON.
C   FURTHER DETAILS ON USE OF THIS CODE ARE AVAILABLE IN *SOLVING
C   ORDINARY DIFFERENTIAL EQUATIONS WITH ODE, STEP, AND INTRP*,
C   BY L. F. SHAMPINE AND M. K. GORDON, SLA-73-1060.
C
C
C   THE PARAMETERS REPRESENT --
C      F -- SUBROUTINE TO EVALUATE DERIVATIVES
C      NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED
C      Y(*) -- SOLUTION VECTOR AT X
C      X -- INDEPENDENT VARIABLE
C      H -- APPROPRIATE STEP SIZE FOR NEXT STEP.  NORMALLY DETERMINED BY
C           CODE
C      EPS -- LOCAL ERROR TOLERANCE
C      WT(*) -- VECTOR OF WEIGHTS FOR ERROR CRITERION
C      START -- LOGICAL VARIABLE SET .TRUE. FOR FIRST STEP,  .FALSE.
C           OTHERWISE
C      HOLD -- STEP SIZE USED FOR LAST SUCCESSFUL STEP
C      K -- APPROPRIATE ORDER FOR NEXT STEP (DETERMINED BY CODE)
C      KOLD -- ORDER USED FOR LAST SUCCESSFUL STEP
C      CRASH -- LOGICAL VARIABLE SET .TRUE. WHEN NO STEP CAN BE TAKEN,
C           .FALSE. OTHERWISE.
C      YP(*) -- DERIVATIVE OF SOLUTION VECTOR AT  X  AFTER SUCCESSFUL
C           STEP
C      KSTEPS -- COUNTER ON ATTEMPTED STEPS
C
C   THE VARIABLES X,XOLD,KOLD,KGI AND IVC AND THE ARRAYS Y,PHI,ALPHA,G,
C   W,P,IV AND GI ARE REQUIRED FOR THE INTERPOLATION SUBROUTINE SINTRP.
C   THE ARRAYS FPWA* AND IFPWA1 AND INTEGER CONSTANTS IFPC* ARE
C   WORKING STORAGE PASSED DIRECTLY THROUGH TO  FODEDS.  THE ARRAYS
C   PAR AND IPAR ARE USER PARAMETERS PASSED THROUGH TO RHOA AND RHOJS.
C
C   INPUT TO STEPS
C
C      FIRST CALL --
C
C   THE USER MUST PROVIDE STORAGE IN HIS CALLING PROGRAM FOR ALL ARRAYS
C   IN THE CALL LIST, NAMELY
C
C     DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),
C    1  ALPHA(12),W(12),G(13),GI(11),IV(10),   FPWA1(NEQN),
C    2  FPWA2(NEQN-1),FPWA3(NEQN-1,NEQN),FPWA4(NEQN-1),
C    3  FPWA5(NEQN),IFPWA1(NEQN)
C                              --                --    **NOTE**
C
C   THE USER MUST ALSO DECLARE  START  AND  CRASH
C   LOGICAL VARIABLES AND  F  AN EXTERNAL SUBROUTINE, SUPPLY THE
C   SUBROUTINE  F(X,Y,YP,FPWA1,FPWA2,FPWA3,IFPC1,IFPWA1,FPWA4,FPWA5,
C                 IFPC2,NEQN-1,IFPC3,PAR,IPAR) TO EVALUATE
C      DY(I)/DX = YP(I) = F(X,Y(1),Y(2),...,Y(NEQN))
C   AND INITIALIZE ONLY THE FOLLOWING PARAMETERS.
C      NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED
C      Y(*) -- VECTOR OF INITIAL VALUES OF DEPENDENT VARIABLES
C      X -- INITIAL VALUE OF THE INDEPENDENT VARIABLE
C      H -- NOMINAL STEP SIZE INDICATING DIRECTION OF INTEGRATION
C           AND MAXIMUM SIZE OF STEP.  MUST BE VARIABLE
C      EPS -- LOCAL ERROR TOLERANCE PER STEP.  MUST BE VARIABLE
C      WT(*) -- VECTOR OF NON-ZERO WEIGHTS FOR ERROR CRITERION
C      START -- .TRUE.
C      KSTEPS -- SET KSTEPS TO ZERO
C   DEFINE U TO BE THE MACHINE UNIT ROUNDOFF QUANTITY BY CALLING
C   THE FUNCTION ROUTINE  D1MACH,  U = D1MACH(3), OR BY
C   COMPUTING U SO THAT U IS THE SMALLEST POSITIVE NUMBER SUCH
C   THAT 1.0+U .GT. 1.0.
C
C   STEPS  REQUIRES THAT THE L2 NORM OF THE VECTOR WITH COMPONENTS
C   LOCAL ERROR(L)/WT(L)  BE LESS THAN  EPS  FOR A SUCCESSFUL STEP.  THE
C   ARRAY  WT  ALLOWS THE USER TO SPECIFY AN ERROR TEST APPROPRIATE
C   FOR HIS PROBLEM.  FOR EXAMPLE,
C      WT(L) = 1.0  SPECIFIES ABSOLUTE ERROR,
C            = ABS(Y(L))  ERROR RELATIVE TO THE MOST RECENT VALUE OF THE
C                 L-TH COMPONENT OF THE SOLUTION,
C            = ABS(YP(L))  ERROR RELATIVE TO THE MOST RECENT VALUE OF
C                 THE L-TH COMPONENT OF THE DERIVATIVE,
C            = MAX(WT(L),ABS(Y(L)))  ERROR RELATIVE TO THE LARGEST
C                 MAGNITUDE OF L-TH COMPONENT OBTAINED SO FAR,
C            = ABS(Y(L))*RELERR/EPS + ABSERR/EPS  SPECIFIES A MIXED
C                 RELATIVE-ABSOLUTE TEST WHERE  RELERR  IS RELATIVE
C                 ERROR,  ABSERR  IS ABSOLUTE ERROR AND  EPS =
C                 MAX(RELERR,ABSERR) .
C
C      SUBSEQUENT CALLS --
C
C   SUBROUTINE  STEPS  IS DESIGNED SO THAT ALL INFORMATION NEEDED TO
C   CONTINUE THE INTEGRATION, INCLUDING THE STEP SIZE  H  AND THE ORDER
C   K , IS RETURNED WITH EACH STEP.  WITH THE EXCEPTION OF THE STEP
C   SIZE, THE ERROR TOLERANCE, AND THE WEIGHTS, NONE OF THE PARAMETERS
C   SHOULD BE ALTERED.  THE ARRAY  WT  MUST BE UPDATED AFTER EACH STEP
C   TO MAINTAIN RELATIVE ERROR TESTS LIKE THOSE ABOVE.  NORMALLY THE
C   INTEGRATION IS CONTINUED JUST BEYOND THE DESIRED ENDPOINT AND THE
C   SOLUTION INTERPOLATED THERE WITH SUBROUTINE  SINTRP .  IF IT IS
C   IMPOSSIBLE TO INTEGRATE BEYOND THE ENDPOINT, THE STEP SIZE MAY BE
C   REDUCED TO HIT THE ENDPOINT SINCE THE CODE WILL NOT TAKE A STEP
C   LARGER THAN THE  H  INPUT.  CHANGING THE DIRECTION OF INTEGRATION,
C   I.E., THE SIGN OF  H , REQUIRES THE USER SET  START = .TRUE. BEFORE
C   CALLING  STEPS  AGAIN.  THIS IS THE ONLY SITUATION IN WHICH  START
C   SHOULD BE ALTERED.
C
C   OUTPUT FROM STEPS
C
C      SUCCESSFUL STEP --
C
C   THE SUBROUTINE RETURNS AFTER EACH SUCCESSFUL STEP WITH  START  AND
C   CRASH  SET .FALSE. .  X  REPRESENTS THE INDEPENDENT VARIABLE
C   ADVANCED ONE STEP OF LENGTH  HOLD  FROM ITS VALUE ON INPUT AND  Y
C   THE SOLUTION VECTOR AT THE NEW VALUE OF  X .  ALL OTHER PARAMETERS
C   REPRESENT INFORMATION CORRESPONDING TO THE NEW  X  NEEDED TO
C   CONTINUE THE INTEGRATION.
C
C      UNSUCCESSFUL STEP --
C
C   WHEN THE ERROR TOLERANCE IS TOO SMALL FOR THE MACHINE PRECISION,
C   THE SUBROUTINE RETURNS WITHOUT TAKING A STEP AND  CRASH = .TRUE. .
C   AN APPROPRIATE STEP SIZE AND ERROR TOLERANCE FOR CONTINUING ARE
C   ESTIMATED AND ALL OTHER INFORMATION IS RESTORED AS UPON INPUT
C   BEFORE RETURNING.  TO CONTINUE WITH THE LARGER TOLERANCE, THE USER
C   JUST CALLS THE CODE AGAIN.  A RESTART IS NEITHER REQUIRED NOR
C   DESIRABLE.
C***REFERENCES  SHAMPINE L.F., GORDON M.K., *SOLVING ORDINARY
C                 DIFFERENTIAL EQUATIONS WITH ODE, STEP, AND INTRP*,
C                 SLA-73-1060, SANDIA LABORATORIES, 1973.
C
      DOUBLE PRECISION ABSH,ALPHA,BETA,D1MACH,EPS,ERK,ERKM1,ERKM2,
     1  ERKP1,ERR,FOURU,FPWA1,FPWA2,FPWA3,FPWA4,FPWA5,G,GI,GSTR,H,
     2  HNEW,HOLD,P,PAR,P5EPS,PHI,PSI,R,REALI,REALNS,RHO,ROUND,SIG,
     3  SUM,TAU,TEMP1,TEMP2,TEMP3,TEMP4,TEMP5,TEMP6,TWO,TWOU,V,
     4  W,WT,X,XOLD,Y,YP
      INTEGER I,IFAIL,IFPC1,IFPC2,IFPC3,IFPWA1,IM1,IPAR,IP1,IQ,IV,
     1  IVC,J,JV,K,KGI,KM1,KM2,KNEW,KOLD,KP1,KP2,KPREV,KSTEPS,
     2  L,LIMIT1,LIMIT2,NEQN,NS,NSM2,NSP1,NSP2
      LOGICAL START,CRASH,PHASE1,NORND
C
      DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12),
     1  ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),GI(11),IV(10),
     2  FPWA1(NEQN),FPWA2(NEQN-1),FPWA3(IFPC1),FPWA4(NEQN-1),
     3  FPWA5(6*NEQN+IFPC1),IFPWA1(NEQN+1),PAR(1),IPAR(1)
      DIMENSION TWO(13),GSTR(13)
C
C   ALL LOCAL VARIABLES ARE SAVED, RATHER THAN PASSED, IN THIS
C   SPECIALIZED VERSION OF STEPS.
C
      SAVE
C
      EXTERNAL F
C
      DATA TWO/2d0,4d0,8d0,16d0,32d0,64d0,128d0,256d0,512d0,1024d0,
     1  2048d0,4096d0,8192d0/
      DATA GSTR/0.500d0,0.0833d0,0.0417d0,0.0264d0,0.0188d0,0.0143d0,
     1  0.0114d0,0.00936d0, 0.00789d0,0.00679d0,0.00592d0,0.00524d0,
     2  0.00468d0/
C
C
C       ***     BEGIN BLOCK 0     ***
C   CHECK IF STEP SIZE OR ERROR TOLERANCE IS TOO SMALL FOR MACHINE
C   PRECISION.  IF FIRST STEP, INITIALIZE PHI ARRAY AND ESTIMATE A
C   STARTING STEP SIZE.
C                   ***
C
C   IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE
C
C***FIRST EXECUTABLE STATEMENT
      TWOU = 2.0 * D1MACH(4)
      FOURU = TWOU + TWOU
      CRASH = .TRUE.
      IF(ABS(H) .GE. FOURU*ABS(X)) GO TO 5
      H = SIGN(FOURU*ABS(X),H)
      RETURN
 5    P5EPS = 0.5*EPS
C
C   IF ERROR TOLERANCE IS TOO SMALL, INCREASE IT TO AN ACCEPTABLE VALUE
C
      ROUND = 0.0
      DO 10 L = 1,NEQN
 10     ROUND = ROUND + (Y(L)/WT(L))**2
      ROUND = TWOU*SQRT(ROUND)
      IF(P5EPS .GE. ROUND) GO TO 15
      EPS = 2.0*ROUND*(1.0 + FOURU)
      RETURN
 15   CRASH = .FALSE.
      G(1) = 1.0
      G(2) = 0.5
      SIG(1) = 1.0
      IF(.NOT.START) GO TO 99
C
C   INITIALIZE.  COMPUTE APPROPRIATE STEP SIZE FOR FIRST STEP
C
      CALL F(X,Y,YP,FPWA1,FPWA2,FPWA3,IFPC1,IFPWA1,FPWA4,FPWA5,
     $       IFPC2,NEQN-1,IFPC3,PAR,IPAR)
      IF (IFPC3 .GT. 0) RETURN
      SUM = 0.0
      DO 20 L = 1,NEQN
        PHI(L,1) = YP(L)
        PHI(L,2) = 0.0
 20     SUM = SUM + (YP(L)/WT(L))**2
      SUM = SQRT(SUM)
      ABSH = ABS(H)
      IF(EPS .LT. 16.0*SUM*H*H) ABSH = 0.25*SQRT(EPS/SUM)
      H = SIGN(MAX(ABSH,FOURU*ABS(X)),H)
C
C*      U = D1MACH(3)
C*      BIG = SQRT(D1MACH(2))
C*      CALL HSTART (F,NEQN,X,X+H,Y,YP,WT,1,U,BIG,
C*     1             PHI(1,3),PHI(1,4),PHI(1,5),PHI(1,6),RPAR,IPAR,H)
C
      HOLD = 0.0
      K = 1
      KOLD = 0
      KPREV = 0
      START = .FALSE.
      PHASE1 = .TRUE.
      NORND = .TRUE.
      IF(P5EPS .GT. 100.0*ROUND) GO TO 99
      NORND = .FALSE.
      DO 25 L = 1,NEQN
 25     PHI(L,15) = 0.0
 99   IFAIL = 0
C       ***     END BLOCK 0     ***
C
C       ***     BEGIN BLOCK 1     ***
C   COMPUTE COEFFICIENTS OF FORMULAS FOR THIS STEP.  AVOID COMPUTING
C   THOSE QUANTITIES NOT CHANGED WHEN STEP SIZE IS NOT CHANGED.
C                   ***
C
 100  KP1 = K+1
      KP2 = K+2
      KM1 = K-1
      KM2 = K-2
C
C   NS IS THE NUMBER OF STEPS TAKEN WITH SIZE H, INCLUDING THE CURRENT
C   ONE.  WHEN K.LT.NS, NO COEFFICIENTS CHANGE
C
      IF(H .NE. HOLD) NS = 0
      IF (NS.LE.KOLD) NS = NS+1
      NSP1 = NS+1
      IF (K .LT. NS) GO TO 199
C
C   COMPUTE THOSE COMPONENTS OF ALPHA(*),BETA(*),PSI(*),SIG(*) WHICH
C   ARE CHANGED
C
      BETA(NS) = 1.0
      REALNS = NS
      ALPHA(NS) = 1.0/REALNS
      TEMP1 = H*REALNS
      SIG(NSP1) = 1.0
      IF(K .LT. NSP1) GO TO 110
      DO 105 I = NSP1,K
        IM1 = I-1
        TEMP2 = PSI(IM1)
        PSI(IM1) = TEMP1
        BETA(I) = BETA(IM1)*PSI(IM1)/TEMP2
        TEMP1 = TEMP2 + H
        ALPHA(I) = H/TEMP1
        REALI = I
 105    SIG(I+1) = REALI*ALPHA(I)*SIG(I)
 110  PSI(K) = TEMP1
C
C   COMPUTE COEFFICIENTS G(*)
C
C   INITIALIZE V(*) AND SET W(*).
C
      IF(NS .GT. 1) GO TO 120
      DO 115 IQ = 1,K
        TEMP3 = IQ*(IQ+1)
        V(IQ) = 1.0/TEMP3
 115    W(IQ) = V(IQ)
      IVC = 0
      KGI = 0
      IF (K .EQ. 1) GO TO 140
      KGI = 1
      GI(1) = W(2)
      GO TO 140
C
C   IF ORDER WAS RAISED, UPDATE DIAGONAL PART OF V(*)
C
 120  IF(K .LE. KPREV) GO TO 130
      IF (IVC .EQ. 0) GO TO 122
      JV = KP1 - IV(IVC)
      IVC = IVC - 1
      GO TO 123
 122  JV = 1
      TEMP4 = K*KP1
      V(K) = 1.0/TEMP4
      W(K) = V(K)
      IF (K .NE. 2) GO TO 123
      KGI = 1
      GI(1) = W(2)
 123  NSM2 = NS-2
      IF(NSM2 .LT. JV) GO TO 130
      DO 125 J = JV,NSM2
        I = K-J
        V(I) = V(I) - ALPHA(J+1)*V(I+1)
 125    W(I) = V(I)
      IF (I .NE. 2) GO TO 130
      KGI = NS - 1
      GI(KGI) = W(2)
C
C   UPDATE V(*) AND SET W(*)
C
 130  LIMIT1 = KP1 - NS
      TEMP5 = ALPHA(NS)
      DO 135 IQ = 1,LIMIT1
        V(IQ) = V(IQ) - TEMP5*V(IQ+1)
 135    W(IQ) = V(IQ)
      G(NSP1) = W(1)
      IF (LIMIT1 .EQ. 1) GO TO 137
      KGI = NS
      GI(KGI) = W(2)
 137  W(LIMIT1+1) = V(LIMIT1+1)
      IF (K .GE. KOLD) GO TO 140
      IVC = IVC + 1
      IV(IVC) = LIMIT1 + 2
C
C   COMPUTE THE G(*) IN THE WORK VECTOR W(*)
C
 140  NSP2 = NS + 2
      KPREV = K
      IF(KP1 .LT. NSP2) GO TO 199
      DO 150 I = NSP2,KP1
        LIMIT2 = KP2 - I
        TEMP6 = ALPHA(I-1)
        DO 145 IQ = 1,LIMIT2
 145      W(IQ) = W(IQ) - TEMP6*W(IQ+1)
 150    G(I) = W(1)
 199    CONTINUE
C       ***     END BLOCK 1     ***
C
C       ***     BEGIN BLOCK 2     ***
C   PREDICT A SOLUTION P(*), EVALUATE DERIVATIVES USING PREDICTED
C   SOLUTION, ESTIMATE LOCAL ERROR AT ORDER K AND ERRORS AT ORDERS K,
C   K-1, K-2 AS IF CONSTANT STEP SIZE WERE USED.
C                   ***
C
C   INCREMENT COUNTER ON ATTEMPTED STEPS
C
      KSTEPS = KSTEPS + 1
C
C   CHANGE PHI TO PHI STAR
C
      IF(K .LT. NSP1) GO TO 215
      DO 210 I = NSP1,K
        TEMP1 = BETA(I)
        DO 205 L = 1,NEQN
 205      PHI(L,I) = TEMP1*PHI(L,I)
 210    CONTINUE
C
C   PREDICT SOLUTION AND DIFFERENCES
C
 215  DO 220 L = 1,NEQN
        PHI(L,KP2) = PHI(L,KP1)
        PHI(L,KP1) = 0.0
 220    P(L) = 0.0
      DO 230 J = 1,K
        I = KP1 - J
        IP1 = I+1
        TEMP2 = G(I)
        DO 225 L = 1,NEQN
          P(L) = P(L) + TEMP2*PHI(L,I)
 225      PHI(L,I) = PHI(L,I) + PHI(L,IP1)
 230    CONTINUE
      IF(NORND) GO TO 240
      DO 235 L = 1,NEQN
        TAU = H*P(L) - PHI(L,15)
        P(L) = Y(L) + TAU
 235    PHI(L,16) = (P(L) - Y(L)) - TAU
      GO TO 250
 240  DO 245 L = 1,NEQN
 245    P(L) = Y(L) + H*P(L)
 250  XOLD = X
      X = X + H
      ABSH = ABS(H)
      CALL F(X,P,YP,FPWA1,FPWA2,FPWA3,IFPC1,IFPWA1,FPWA4,FPWA5,
     $       IFPC2,NEQN-1,IFPC3,PAR,IPAR)
      IF (IFPC3 .GT. 0) RETURN
C
C   ESTIMATE ERRORS AT ORDERS K,K-1,K-2
C
      ERKM2 = 0.0
      ERKM1 = 0.0
      ERK = 0.0
      DO 265 L = 1,NEQN
        TEMP3 = 1.0/WT(L)
        TEMP4 = YP(L) - PHI(L,1)
        IF(KM2)265,260,255
 255    ERKM2 = ERKM2 + ((PHI(L,KM1)+TEMP4)*TEMP3)**2
 260    ERKM1 = ERKM1 + ((PHI(L,K)+TEMP4)*TEMP3)**2
 265    ERK = ERK + (TEMP4*TEMP3)**2
      IF(KM2)280,275,270
 270  ERKM2 = ABSH*SIG(KM1)*GSTR(KM2)*SQRT(ERKM2)
 275  ERKM1 = ABSH*SIG(K)*GSTR(KM1)*SQRT(ERKM1)
 280  TEMP5 = ABSH*SQRT(ERK)
      ERR = TEMP5*(G(K)-G(KP1))
      ERK = TEMP5*SIG(KP1)*GSTR(K)
      KNEW = K
C
C   TEST IF ORDER SHOULD BE LOWERED
C
      IF(KM2)299,290,285
 285  IF(MAX(ERKM1,ERKM2) .LE. ERK) KNEW = KM1
      GO TO 299
 290  IF(ERKM1 .LE. 0.5*ERK) KNEW = KM1
C
C   TEST IF STEP SUCCESSFUL
C
 299  IF(ERR .LE. EPS) GO TO 400
C       ***     END BLOCK 2     ***
C
C       ***     BEGIN BLOCK 3     ***
C   THE STEP IS UNSUCCESSFUL.  RESTORE  X, PHI(*,*), PSI(*) .
C   IF THIRD CONSECUTIVE FAILURE, SET ORDER TO ONE.  IF STEP FAILS MORE
C   THAN THREE TIMES, CONSIDER AN OPTIMAL STEP SIZE.  DOUBLE ERROR
C   TOLERANCE AND RETURN IF ESTIMATED STEP SIZE IS TOO SMALL FOR MACHINE
C   PRECISION.
C                   ***
C
C   RESTORE X, PHI(*,*) AND PSI(*)
C
      PHASE1 = .FALSE.
      X = XOLD
      DO 310 I = 1,K
        TEMP1 = 1.0/BETA(I)
        IP1 = I+1
        DO 305 L = 1,NEQN
 305      PHI(L,I) = TEMP1*(PHI(L,I) - PHI(L,IP1))
 310    CONTINUE
      IF(K .LT. 2) GO TO 320
      DO 315 I = 2,K
 315    PSI(I-1) = PSI(I) - H
C
C   ON THIRD FAILURE, SET ORDER TO ONE.  THEREAFTER, USE OPTIMAL STEP
C   SIZE
C
 320  IFAIL = IFAIL + 1
      TEMP2 = 0.5
      IF(IFAIL - 3) 335,330,325
 325  IF(P5EPS .LT. 0.25*ERK) TEMP2 = SQRT(P5EPS/ERK)
 330  KNEW = 1
 335  H = TEMP2*H
      K = KNEW
      NS = 0
      IF(ABS(H) .GE. FOURU*ABS(X)) GO TO 340
      CRASH = .TRUE.
      H = SIGN(FOURU*ABS(X),H)
      EPS = EPS + EPS
      RETURN
 340  GO TO 100
C       ***     END BLOCK 3     ***
C
C       ***     BEGIN BLOCK 4     ***
C   THE STEP IS SUCCESSFUL.  CORRECT THE PREDICTED SOLUTION, EVALUATE
C   THE DERIVATIVES USING THE CORRECTED SOLUTION AND UPDATE THE
C   DIFFERENCES.  DETERMINE BEST ORDER AND STEP SIZE FOR NEXT STEP.
C                   ***
 400  KOLD = K
      HOLD = H
C
C   CORRECT AND EVALUATE
C
      TEMP1 = H*G(KP1)
      IF(NORND) GO TO 410
      DO 405 L = 1,NEQN
        TEMP3 = Y(L)
        RHO = TEMP1*(YP(L) - PHI(L,1)) - PHI(L,16)
        Y(L) = P(L) + RHO
        PHI(L,15) = (Y(L) - P(L)) - RHO
 405    P(L) = TEMP3
      GO TO 420
 410  DO 415 L = 1,NEQN
        TEMP3 = Y(L)
        Y(L) = P(L) + TEMP1*(YP(L) - PHI(L,1))
 415    P(L) = TEMP3
 420  CALL F(X,Y,YP,FPWA1,FPWA2,FPWA3,IFPC1,IFPWA1,FPWA4,FPWA5,
     $       IFPC2,NEQN-1,IFPC3,PAR,IPAR)
      IF (IFPC3 .GT. 0) RETURN
C
C   UPDATE DIFFERENCES FOR NEXT STEP
C
      DO 425 L = 1,NEQN
        PHI(L,KP1) = YP(L) - PHI(L,1)
 425    PHI(L,KP2) = PHI(L,KP1) - PHI(L,KP2)
      DO 435 I = 1,K
        DO 430 L = 1,NEQN
 430      PHI(L,I) = PHI(L,I) + PHI(L,KP1)
 435    CONTINUE
C
C   ESTIMATE ERROR AT ORDER K+1 UNLESS:
C     IN FIRST PHASE WHEN ALWAYS RAISE ORDER,
C     ALREADY DECIDED TO LOWER ORDER,
C     STEP SIZE NOT CONSTANT SO ESTIMATE UNRELIABLE
C
      ERKP1 = 0.0
      IF(KNEW .EQ. KM1  .OR.  K .EQ. 12) PHASE1 = .FALSE.
      IF(PHASE1) GO TO 450
      IF(KNEW .EQ. KM1) GO TO 455
      IF(KP1 .GT. NS) GO TO 460
      DO 440 L = 1,NEQN
 440    ERKP1 = ERKP1 + (PHI(L,KP2)/WT(L))**2
      ERKP1 = ABSH*GSTR(KP1)*SQRT(ERKP1)
C
C   USING ESTIMATED ERROR AT ORDER K+1, DETERMINE APPROPRIATE ORDER
C   FOR NEXT STEP
C
      IF(K .GT. 1) GO TO 445
      IF(ERKP1 .GE. 0.5*ERK) GO TO 460
      GO TO 450
 445  IF(ERKM1 .LE. MIN(ERK,ERKP1)) GO TO 455
      IF(ERKP1 .GE. ERK  .OR.  K .EQ. 12) GO TO 460
C
C   HERE ERKP1 .LT. ERK .LT. MAX(ERKM1,ERKM2) ELSE ORDER WOULD HAVE
C   BEEN LOWERED IN BLOCK 2.  THUS ORDER IS TO BE RAISED
C
C   RAISE ORDER
C
 450  K = KP1
      ERK = ERKP1
      GO TO 460
C
C   LOWER ORDER
C
 455  K = KM1
      ERK = ERKM1
C
C   WITH NEW ORDER DETERMINE APPROPRIATE STEP SIZE FOR NEXT STEP
C
 460  HNEW = H + H
      IF(PHASE1) GO TO 465
      IF(P5EPS .GE. ERK*TWO(K+1)) GO TO 465
      HNEW = H
      IF(P5EPS .GE. ERK) GO TO 465
      TEMP2 = K+1
      R = (P5EPS/ERK)**(1.0/TEMP2)
      HNEW = ABSH*MAX(0.5D0,MIN(0.9D0,R))
      HNEW = SIGN(MAX(HNEW,FOURU*ABS(X)),H)
 465  H = HNEW
      RETURN
C       ***     END BLOCK 4     ***
      END

ext: f

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg