text:
module test_atmos_dyn

  !-----------------------------------------------------------------------------
  use scale_precision
  use scale_atmos_grid_cartesC_index
  use scale_index
  use scale_tracer
  use scale_atmos_dyn, only: &
     ATMOS_DYN
  use dc_test, only: &
     AssertEqual, &
     AssertGreaterThan, &
     AssertLessThan
  use scale_io, only: &
     H_SHORT
  use scale_comm_cartesC, only: &
     COMM_vars8, &
     COMM_wait
  use scale_atmos_grid_cartesC, only: &
     DOMAIN_CENTER_Y => ATMOS_GRID_CARTESC_DOMAIN_CENTER_Y, &
     CY              => ATMOS_GRID_CARTESC_CY,              &
     CZ              => ATMOS_GRID_CARTESC_CZ,              &
     FZ              => ATMOS_GRID_CARTESC_FZ,              &
     CDZ             => ATMOS_GRID_CARTESC_CDZ,             &
     CDX             => ATMOS_GRID_CARTESC_CDX,             &
     CDY             => ATMOS_GRID_CARTESC_CDY,             &
     FDZ             => ATMOS_GRID_CARTESC_FDZ,             &
     FDX             => ATMOS_GRID_CARTESC_FDX,             &
     FDY             => ATMOS_GRID_CARTESC_FDY,             &
     RCDZ            => ATMOS_GRID_CARTESC_RCDZ,            &
     RCDX            => ATMOS_GRID_CARTESC_RCDX,            &
     RCDY            => ATMOS_GRID_CARTESC_RCDY,            &
     RFDZ            => ATMOS_GRID_CARTESC_RFDZ,            &
     RFDX            => ATMOS_GRID_CARTESC_RFDX,            &
     RFDY            => ATMOS_GRID_CARTESC_RFDY
  use scale_atmos_hydrometeor, only: &
     I_QV
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  public :: test_atmos_dyn_run

  !-----------------------------------------------------------------------------
  real(RP), allocatable :: DENS(:,:,:)
  real(RP), allocatable :: MOMZ(:,:,:)
  real(RP), allocatable :: MOMX(:,:,:)
  real(RP), allocatable :: MOMY(:,:,:)
  real(RP), allocatable :: RHOT(:,:,:)
  real(RP), allocatable :: QTRC(:,:,:,:)
  real(RP), allocatable :: PROG(:,:,:,:)

  real(RP), allocatable :: DENS_av(:,:,:)
  real(RP), allocatable :: MOMZ_av(:,:,:)
  real(RP), allocatable :: MOMX_av(:,:,:)
  real(RP), allocatable :: MOMY_av(:,:,:)
  real(RP), allocatable :: RHOT_av(:,:,:)
  real(RP), allocatable :: QTRC_av(:,:,:,:)

  real(RP), allocatable :: DENS_tp(:,:,:)
  real(RP), allocatable :: MOMZ_tp(:,:,:)
  real(RP), allocatable :: MOMX_tp(:,:,:)
  real(RP), allocatable :: MOMY_tp(:,:,:)
  real(RP), allocatable :: RHOT_tp(:,:,:)
  real(RP), allocatable :: QTRC_tp(:,:,:,:)

  real(RP), allocatable :: DENS_o(:,:,:)
  real(RP), allocatable :: MOMZ_o(:,:,:)
  real(RP), allocatable :: MOMX_o(:,:,:)
  real(RP), allocatable :: MOMY_o(:,:,:)
  real(RP), allocatable :: RHOT_o(:,:,:)
  real(RP), allocatable :: QTRC_o(:,:,:,:)

  real(RP), allocatable :: CORIOLIS(:,:)

  real(RP), allocatable :: REF_dens(:,:,:)
  real(RP), allocatable :: REF_pott(:,:,:)
  real(RP), allocatable :: REF_qv  (:,:,:)
  real(RP), allocatable :: REF_pres(:,:,:)

  real(RP), allocatable :: DAMP_var(:,:,:,:)
  real(RP), allocatable :: DAMP_alpha(:,:,:,:)
  real(RP), allocatable :: MFLUX_OFFSET_X(:,:,:)
  real(RP), allocatable :: MFLUX_OFFSET_Y(:,:,:)

  real(RP), allocatable :: PHI(:,:,:)
  real(RP), allocatable :: GSQRT(:,:,:,:)
  real(RP), allocatable :: J13G(:,:,:,:)
  real(RP), allocatable :: J23G(:,:,:,:)
  real(RP) :: J33G
  real(RP), allocatable :: MAPF(:,:,:,:)

  real(RP), allocatable :: AQ_R(:)
  real(RP), allocatable :: AQ_CV(:)
  real(RP), allocatable :: AQ_CP(:)
  real(RP), allocatable :: AQ_MASS(:)

  integer  :: BND_QA
  integer, allocatable  :: BND_IQ(:)
  real(RP) :: BND_SMOOTHER_FACT

  integer  :: nd_order
  real(RP) :: nd_coef
  real(RP) :: nd_sfc_fact
  logical  :: nd_use_rs

  real(RP) :: divdmp_coef

  logical  :: flag_tracer_split_tend = .false.

  logical  :: flag_fct_momentum = .true.
  logical  :: flag_fct_t        = .true.
  logical  :: flag_fct_tracer   = .true.
  logical  :: flag_fct_along_stream = .true.

  real(RP), allocatable :: ZERO(:,:,:)

  integer :: KME ! end of main region

  character(len=H_SHORT) :: DYN_TYPE

  character(len=H_SHORT) :: DYN_Tinteg_Short_TYPE
  character(len=H_SHORT) :: DYN_Tinteg_Tracer_TYPE
  character(len=H_SHORT) :: DYN_Tinteg_Large_TYPE
  character(len=H_SHORT) :: DYN_Tstep_Tracer_TYPE
  character(len=H_SHORT) :: DYN_Tstep_Large_TYPE
  character(len=H_SHORT) :: DYN_FVM_FLUX_TYPE
  character(len=H_SHORT) :: DYN_FVM_FLUX_TYPE_TRACER
  real(RP)               :: wdamp_tau
  real(RP)               :: wdamp_height

  integer :: k, i, j, iq
  character(len=11) :: message
  !-----------------------------------------------------------------------------
contains

  subroutine test_atmos_dyn_run
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_io
  use scale_atmos_dyn, only: &
     ATMOS_DYN_setup
  use scale_atmos_dyn_Tstep_short, only: &
     ATMOS_DYN_Tstep_short_regist
  use scale_atmos_grid_cartesC, only: &
     CBFZ => ATMOS_GRID_CARTESC_CBFZ
  use scale_const, only: &
     GRAV => CONST_GRAV
  !-----------------------------------------------------------------------------
  implicit none
  !-----------------------------------------------------------------------------
  !
  !++ parameters & variables
  !
  !-----------------------------------------------------------------------------
  real(RP) :: lat(IA,JA)
  character(len=H_SHORT) :: CSDUMMY(1)
  character(len=H_MID)   :: CMDUMMY(1)
  integer :: j
  !=============================================================================

  !########## Initial setup ##########

  ! allocate
  allocate( DENS(KA,IA,JA) )
  allocate( MOMZ(KA,IA,JA) )
  allocate( MOMX(KA,IA,JA) )
  allocate( MOMY(KA,IA,JA) )
  allocate( RHOT(KA,IA,JA) )
  allocate( QTRC(KA,IA,JA,QA) )

  allocate( DENS_av(KA,IA,JA) )
  allocate( MOMZ_av(KA,IA,JA) )
  allocate( MOMX_av(KA,IA,JA) )
  allocate( MOMY_av(KA,IA,JA) )
  allocate( RHOT_av(KA,IA,JA) )
  allocate( QTRC_av(KA,IA,JA,QA) )

  allocate( DENS_tp(KA,IA,JA) )
  allocate( MOMZ_tp(KA,IA,JA) )
  allocate( MOMX_tp(KA,IA,JA) )
  allocate( MOMY_tp(KA,IA,JA) )
  allocate( RHOT_tp(KA,IA,JA) )
  allocate( QTRC_tp(KA,IA,JA,QA) )

  allocate( DENS_o(KA,IA,JA) )
  allocate( MOMZ_o(KA,IA,JA) )
  allocate( MOMX_o(KA,IA,JA) )
  allocate( MOMY_o(KA,IA,JA) )
  allocate( RHOT_o(KA,IA,JA) )
  allocate( QTRC_o(KA,IA,JA,QA) )

  allocate( CORIOLIS(IA,JA) )

  allocate( REF_dens(KA,IA,JA) )
  allocate( REF_pott(KA,IA,JA) )
  allocate( REF_qv  (KA,IA,JA) )
  allocate( REF_pres(KA,IA,JA) )

  allocate( DAMP_var(KA,IA,JA,5+QA) )
  allocate( DAMP_alpha(KA,IA,JA,5+QA) )
  allocate( MFLUX_OFFSET_X(KA,JA,2) )
  allocate( MFLUX_OFFSET_Y(KA,IA,2) )

  allocate( PHI(KA,IA,JA) )
  allocate( GSQRT(KA,IA,JA,7) )
  allocate( J13G(KA,IA,JA,7) )
  allocate( J23G(KA,IA,JA,7) )

  allocate( MAPF(IA,JA,2,4) )

  allocate( AQ_R(QA) )
  allocate( AQ_CV(QA) )
  allocate( AQ_CP(QA) )
  allocate( AQ_MASS(QA) )

  allocate( ZERO(KA,IA,JA) )

  allocate( PROG(KA,IA,JA,1) )

  allocate( BND_IQ(QA) )
  BND_QA            = 0
  BND_IQ(:)         = -1 ! dummy
  BND_SMOOTHER_FACT = 0.2_RP

  ZERO(:,:,:) = 0.0_RP

  CORIOLIS(:,:) = 0.0_RP

  nd_order = 2
  nd_coef = 0.01_RP
  nd_sfc_fact = 1.0_RP
  nd_use_rs = .true.
  do j = 1, JA
     lat(:,j) = real(j, RP)
  end do

  DYN_TYPE = "FVM-HEVE"
  call ATMOS_DYN_Tstep_short_regist( DYN_TYPE, & !(in)
                                     VA, CSDUMMY, CMDUMMY, CSDUMMY ) ! (out)

  DYN_Tinteg_Short_TYPE = "RK4"
  DYN_Tinteg_Tracer_TYPE = "RK3WS2002"
  DYN_Tinteg_Large_TYPE = "EULER"
  DYN_Tstep_Tracer_TYPE = "FVM-HEVE"
  DYN_Tstep_Large_TYPE = "FVM-HEVE"
  DYN_FVM_FLUX_TYPE = "CD4"
  DYN_FVM_FLUX_TYPE_TRACER = "UD3KOREN1993"
  wdamp_tau = 10.0_RP
  wdamp_height = 0.0_RP

  call ATMOS_DYN_setup( &
       DYN_Tinteg_Short_TYPE,              & ! (in)
       DYN_Tinteg_Tracer_TYPE,             & ! (in)
       DYN_Tinteg_Large_TYPE,              & ! (in)
       DYN_Tstep_Tracer_TYPE,              & ! (in)
       DYN_Tstep_Large_TYPE,               & ! (in)
       DYN_TYPE,                           & ! (in)
       DYN_FVM_FLUX_TYPE,                  & ! (in)
       DYN_FVM_FLUX_TYPE_TRACER,           & ! (in)
       DENS, MOMZ, MOMX, MOMY, RHOT, QTRC, & ! (in)
       PROG,                               & ! (in)
       CDZ, CDX, CDY, FDZ, FDX, FDY,       & ! (in)
       wdamp_tau, wdamp_height, FZ         ) ! (in)

  do k = KS+1, KE
     if ( CBFZ(k) > 0.0_RP ) then
        KME = k - 1
        exit
     end if
  end do

  MOMZ(KE,:,:) = 0.0_RP

  AQ_R(:) = 1.0_RP
  AQ_CV(:) = 1.0_RP
  AQ_CP(:) = 1.0_RP
  AQ_MASS(:) = 1.0_RP

  do j = 1, JA
  do i = 1, IA
     PHI(:,i,j) = CZ(:) * GRAV
  end do
  end do
  GSQRT(:,:,:,:) = 1.0_RP
  J13G(:,:,:,:) = 0.0_RP
  J23G(:,:,:,:) = 0.0_RP
  J33G          = 1.0_RP
  MAPF(:,:,:,:) = 1.0_RP

  divdmp_coef = 0.0_RP

  DENS_tp(:,:,:) = 0.0_RP
  MOMZ_tp(:,:,:) = 0.0_RP
  MOMX_tp(:,:,:) = 0.0_RP
  MOMY_tp(:,:,:) = 0.0_RP
  RHOT_tp(:,:,:) = 0.0_RP
  QTRC_tp(:,:,:,:) = 0.0_RP

  !########## test ##########

  call test_undef

  call test_const

  call test_conserve

  call test_cwc

  call test_fctminmax

end subroutine test_atmos_dyn_run
!=============================================================================

subroutine test_undef

  real(RP) :: BIG(KA,IA,JA)

  BIG(:,:,:) = 9.99E9_RP

  write(*,*) "Test undef"

  do j = 1, JA
  do i = 1, IA
  do k = 1, KA
     MOMZ(k,i,j) = k + i + j
     MOMX(k,i,j) = k + i - j
     MOMY(k,i,j) = k - i + j
     RHOT(k,i,j) = k - i - j + KA + IA + JA
     DENS(k,i,j) = k + i + j + 1.0_RP
     do iq = 1, QA
        QTRC(k,i,j,iq) = 1.0_RP*(k + i + j + iq) / (KA + IA + JA + QA) / QA
     end do
  end do
  end do
  end do

  do j = 1, JA
  do i = 1, IA
  do k = 1, KA
     REF_dens(k,i,j) = KA - k + 1
     REF_pott(k,i,j) = k
     REF_qv  (k,i,j) = KA - k + 1
     REF_pres(k,i,j) = KA - k + 1
  end do
  end do
  end do

  DAMP_var  (:,:,:,:) = -9.999E30_RP
  DAMP_alpha(:,:,:,:) = 0.0_RP

  MFLUX_OFFSET_X(:,:,:) = 0.0_RP
  MFLUX_OFFSET_Y(:,:,:) = 0.0_RP

  do i = 1, 2
     call ATMOS_DYN( &
          DENS, MOMZ, MOMX, MOMY, RHOT, QTRC,          & ! (inout)
          PROG,                                        & ! (inout)
          DENS_av, MOMZ_av, MOMX_av, MOMY_av, RHOT_av, QTRC_av, & ! (out)
          DENS_tp, MOMZ_tp, MOMX_tp, MOMY_tp, RHOT_tp, QTRC_tp, & ! (in)
          CORIOLIS,                                    & ! (in)
          CDZ, CDX, CDY, FDZ, FDX, FDY,                & ! (in)
          RCDZ, RCDX, RCDY, RFDZ, RFDX, RFDY,          & ! (in)
          PHI, GSQRT, J13G, J23G, J33G, MAPF,          & ! (in)
          AQ_R, AQ_CV, AQ_CP, AQ_MASS,                 & ! (in)
          REF_dens, REF_pott, REF_qv, REF_pres,        & ! (in)
          nd_coef, nd_coef, nd_order, nd_sfc_fact, nd_use_rs, & ! (in)
          BND_QA, BND_IQ, BND_SMOOTHER_FACT,                  & ! (in)
          DAMP_var(:,:,:,1), DAMP_var(:,:,:,2), DAMP_var(:,:,:,3), DAMP_var(:,:,:,4), DAMP_var(:,:,:,5), DAMP_var(:,:,:,6:6+QA-1), & ! (in)
          DAMP_alpha(:,:,:,1), DAMP_alpha(:,:,:,2), DAMP_alpha(:,:,:,3), DAMP_alpha(:,:,:,4), DAMP_alpha(:,:,:,5), & ! (in)
          DAMP_alpha(:,:,:,6:6+QA-1),                  & ! (in)
          MFLUX_OFFSET_X(:,:,:), MFLUX_OFFSET_Y(:,:,:), & ! (in)
          divdmp_coef,                                 & ! (in)
          flag_tracer_split_tend,                      & ! (in)
          flag_fct_momentum, flag_fct_t, flag_fct_tracer, & ! (in)
          flag_fct_along_stream,                       & ! (in)
          .false.,                                     & ! (in)
          I_QV,                                        & ! (in)
          1.0_DP, 1.0_DP                               ) ! (in)

  end do

  call AssertLessThan("MOMZ", BIG(KS:KE,IS:IE,JS:JE), MOMZ(KS:KE,IS:IE,JS:JE))
  call AssertLessThan("MOMX", BIG(KS:KE,IS:IE,JS:JE), MOMX(KS:KE,IS:IE,JS:JE))
  call AssertLessThan("MOMY", BIG(KS:KE,IS:IE,JS:JE), MOMY(KS:KE,IS:IE,JS:JE))
  call AssertLessThan("RHOT", BIG(KS:KE,IS:IE,JS:JE), RHOT(KS:KE,IS:IE,JS:JE))
  message = "iq = ??"
  do iq = 1, QA
     write(message(6:7), "(i2)") iq
     call AssertLessThan(message, BIG(KS:KE,IS:IE,JS:JE), QTRC(KS:KE,IS:IE,JS:JE,iq))
  end do

end subroutine test_undef

subroutine test_const
  real(RP) :: answer(KA,IA,JA)

  write(*,*) "Test constant"

  DENS(:,:,:) = 1.0_RP
  MOMZ(:,:,:) = 2.0_RP
  MOMX(:,:,:) = 3.0_RP
  MOMY(:,:,:) = 4.0_RP
  RHOT(:,:,:) = 300.0_RP
  QTRC(:,:,:,:) = 0.1_RP

  REF_dens(:,:,:) = 1.0_RP
  REF_pott(:,:,:) = 300.0_RP
  REF_qv(:,:,:)   = 0.001_RP
  REF_pres(:,:,:) = 1000._RP

  DAMP_var  (:,:,:,:) = -9.999E30_RP
  DAMP_alpha(:,:,:,:) = 0.0_RP

  MFLUX_OFFSET_X(:,:,:) = 0.0_RP
  MFLUX_OFFSET_Y(:,:,:) = 0.0_RP

  call ATMOS_DYN( &
       DENS, MOMZ, MOMX, MOMY, RHOT, QTRC,          & ! (inout)
       PROG,                                        & ! (inout)
       DENS_av, MOMZ_av, MOMX_av, MOMY_av, RHOT_av, QTRC_av, & ! (out)
       DENS_tp, MOMZ_tp, MOMX_tp, MOMY_tp, RHOT_tp, QTRC_tp, & ! (in)
       CORIOLIS,                                    & ! (in)
       CDZ, CDX, CDY, FDZ, FDX, FDY,                & ! (in)
       RCDZ, RCDX, RCDY, RFDZ, RFDX, RFDY,          & ! (in)
       PHI, GSQRT, J13G, J23G, J33G, MAPF,          & ! (in)
       AQ_R, AQ_CV, AQ_CP, AQ_MASS,                 & ! (in)
       REF_dens, REF_pott, REF_qv, REF_pres,        & ! (in)
       nd_coef, nd_coef, nd_order, nd_sfc_fact, nd_use_rs, & ! (in)
       BND_QA, BND_IQ, BND_SMOOTHER_FACT,                  & ! (in)
       DAMP_var(:,:,:,1), DAMP_var(:,:,:,2), DAMP_var(:,:,:,3), DAMP_var(:,:,:,4), DAMP_var(:,:,:,5), DAMP_var(:,:,:,6:6+QA-1), & ! (in)
       DAMP_alpha(:,:,:,1), DAMP_alpha(:,:,:,2), DAMP_alpha(:,:,:,3), DAMP_alpha(:,:,:,4), DAMP_alpha(:,:,:,5), & ! (in)
       DAMP_alpha(:,:,:,6:6+QA-1),                  & ! (in)
       MFLUX_OFFSET_X(:,:,:), MFLUX_OFFSET_Y(:,:,:), & ! (in)
       divdmp_coef,                                 & ! (in)
       flag_tracer_split_tend,                      & ! (in)
       flag_fct_momentum, flag_fct_t, flag_fct_tracer, & ! (in)
       flag_fct_along_stream,                       & ! (in)
       .false.,                                     & ! (in)
       I_QV,                                        & ! (in)
       1.0_DP, 1.0_DP                               ) ! (in)

  do k = KS, KE
     answer(k,:,:) = MOMZ(k,IS,JS)
  end do
  call AssertEqual("MOMZ", answer(KS:KE,IS:IE,JS:JE), MOMZ(KS:KE,IS:IE,JS:JE))
  do k = KS, KE
     answer(k,:,:) = MOMX(k,IS,JS)
  end do
  call AssertEqual("MOMX", answer(KS:KE,IS:IE,JS:JE), MOMX(KS:KE,IS:IE,JS:JE))
  do k = KS, KE
     answer(k,:,:) = MOMY(k,IS,JS)
  end do
  call AssertEqual("MOMY", answer(KS:KE,IS:IE,JS:JE), MOMY(KS:KE,IS:IE,JS:JE))
  do k = KS, KE
     answer(k,:,:) = DENS(k,IS,JS)
  end do
  call AssertEqual("DENS", answer(KS:KE,IS:IE,JS:JE), DENS(KS:KE,IS:IE,JS:JE))
  do k = KS, KE
     answer(k,:,:) = RHOT(k,IS,JS)
  end do
  call AssertEqual("RHOT", answer(KS:KE,IS:IE,JS:JE), RHOT(KS:KE,IS:IE,JS:JE))
  message = "iq = ??"
  do iq = 1, QA
     do k = KS, KE
        answer(k,:,:) = QTRC(k,IS,JS,iq)
     end do
     write(message(6:7), "(i2)") iq
     call AssertEqual(message, answer(KS:KE,IS:IE,JS:JE), QTRC(KS:KE,IS:IE,JS:JE,iq))
  end do

end subroutine test_const

subroutine test_conserve

  real(RP) :: total_o, total

  write(*,*) "Test conserve"

  do j = 1, JA
  do i = 1, IA
  do k = KS, KE
     MOMZ(k,i,j) = (k + i + j) * 1.0_RP
     MOMX(k,i,j) = (k*2 + i + j) * 1.0_RP
     MOMY(k,i,j) = (k + i*2 + j) * 1.0_RP
     RHOT(k,i,j) = (k + i + j*2) * 1.0_RP
     DENS(k,i,j) = (k*2 + i*2 + j) * 1.0_RP + 1.0_RP
     do iq = 1, QA
        QTRC(k,i,j,iq) = (k + i + j + iq) * 0.1_RP / (KE + IE + JE + QA)
     end do
  end do
  end do
  end do

  REF_dens(:,:,:) = 1.0_RP
  REF_pott(:,:,:) = 1.0_RP
  REF_qv(:,:,:)   = 1.0_RP
  REF_pres(:,:,:) = 1.0_RP

  DAMP_var  (:,:,:,:) = -9.999E30_RP
  DAMP_alpha(:,:,:,:) = 0.0_RP

  MFLUX_OFFSET_X(:,:,:) = 0.0_RP
  MFLUX_OFFSET_Y(:,:,:) = 0.0_RP

  call COMM_vars8( DENS(:,:,:), 1 )
  call COMM_vars8( MOMZ(:,:,:), 2 )
  call COMM_vars8( MOMX(:,:,:), 3 )
  call COMM_vars8( MOMY(:,:,:), 4 )
  call COMM_vars8( RHOT(:,:,:), 5 )
  call COMM_wait ( DENS(:,:,:), 1 )
  call COMM_wait ( MOMZ(:,:,:), 2 )
  call COMM_wait ( MOMX(:,:,:), 3 )
  call COMM_wait ( MOMY(:,:,:), 4 )
  call COMM_wait ( RHOT(:,:,:), 5 )

  do iq = 1, QA
     call COMM_vars8( QTRC(:,:,:,iq), iq )
  enddo
  do iq = 1, QA
     call COMM_wait ( QTRC(:,:,:,iq), iq )
  enddo

  call copy

  call ATMOS_DYN( &
         DENS, MOMZ, MOMX, MOMY, RHOT, QTRC,          & ! (out)
         PROG,                                        & ! (inout)
         DENS_av, MOMZ_av, MOMX_av, MOMY_av, RHOT_av, QTRC_av, & ! (inout)
         DENS_tp, MOMZ_tp, MOMX_tp, MOMY_tp, RHOT_tp, QTRC_tp, & ! (in)
         CORIOLIS,                                    & ! (in)
         CDZ, CDX, CDY, FDZ, FDX, FDY,                & ! (in)
         RCDZ, RCDX, RCDY, RFDZ, RFDX, RFDY,          & ! (in)
         PHI, GSQRT, J13G, J23G, J33G, MAPF,          & ! (in)
         AQ_R, AQ_CV, AQ_CP, AQ_MASS,                 & ! (in)
         REF_dens, REF_pott, REF_qv, REF_pres,        & ! (in)
         nd_coef, nd_coef, nd_order, nd_sfc_fact, nd_use_rs, & ! (in)
         BND_QA, BND_IQ, BND_SMOOTHER_FACT,                  & ! (in)
         DAMP_var(:,:,:,1), DAMP_var(:,:,:,2), DAMP_var(:,:,:,3), DAMP_var(:,:,:,4), DAMP_var(:,:,:,5), DAMP_var(:,:,:,6:6+QA-1), & ! (in)
         DAMP_alpha(:,:,:,1), DAMP_alpha(:,:,:,2), DAMP_alpha(:,:,:,3), DAMP_alpha(:,:,:,4), DAMP_alpha(:,:,:,5), & ! (in)
         DAMP_alpha(:,:,:,6:6+QA-1),                  & ! (in)
         MFLUX_OFFSET_X(:,:,:), MFLUX_OFFSET_Y(:,:,:), & ! (in)
         divdmp_coef,                                 & ! (in)
         flag_tracer_split_tend,                      & ! (in)
         flag_fct_momentum, flag_fct_t, flag_fct_tracer, & ! (in)
         flag_fct_along_stream,                       & ! (in)
         .true.,                                      & ! (in)
         I_QV,                                        & ! (in)
         1.0_DP, 1.0_DP                               ) ! (in)

  total_o = 0.0_RP
  total = 0.0_RP
  do j = JS, JE
  do i = IS, IE
  do k = KS, KE
     total_o = total_o + DENS_o(k,i,j) * CDZ(k) * CDX(i) * CDY(j)
     total   = total   + DENS  (k,i,j) * CDZ(k) * CDX(i) * CDY(j)
  end do
  end do
  end do
  call AssertEqual("DENS", total_o, total, RP*2-2, -10)

  total_o = 0.0_RP
  total = 0.0_RP
  do j = JS, JE
  do i = IS, IE
  do k = KS, KE
     total_o = total_o + MOMX_o(k,i,j) * CDZ(k) * CDX(i) * CDY(j)
     total   = total   + MOMX  (k,i,j) * CDZ(k) * CDX(i) * CDY(j)
  end do
  end do
  end do
  call AssertEqual("MOMX", total_o, total, RP*2-2, -10)

  total_o = 0.0_RP
  total = 0.0_RP
  do j = JS, JE
  do i = IS, IE
  do k = KS, KE
     total_o = total_o + MOMY_o(k,i,j) * CDZ(k) * CDX(i) * CDY(j)
     total   = total   + MOMY  (k,i,j) * CDZ(k) * CDX(i) * CDY(j)
  end do
  end do
  end do
  call AssertEqual("MOMY", total_o, total, RP*2-2, -10)

  total_o = 0.0_RP
  total = 0.0_RP
  do j = JS, JE
  do i = IS, IE
  do k = KS, KE
     total_o = total_o + RHOT_o(k,i,j) * CDZ(k) * CDX(i) * CDY(j)
     total   = total   + RHOT  (k,i,j) * CDZ(k) * CDX(i) * CDY(j)
  end do
  end do
  end do
  call AssertEqual("RHOT", total_o, total, RP*2-2, -10)

  message = "iq = ??"
  do iq = 1, QA
     write(message(6:7), "(i2)") iq
     total_o = 0.0_RP
     total = 0.0_RP
     do j = JS, JE
     do i = IS, IE
     do k = KS, KE
        total_o = total_o + RHOT_o(k,i,j) * CDZ(k) * CDX(i) * CDY(j)
        total   = total   + RHOT  (k,i,j) * CDZ(k) * CDX(i) * CDY(j)
     end do
     end do
     end do
     call AssertEqual(message, total_o, total, RP*2-2, -10)
  end do

end subroutine test_conserve

subroutine test_cwc
  real(RP), parameter :: POTT = 300._RP
  real(RP), parameter :: Q = 0.01_RP
  real(RP) :: answer(KA,IA,JA)

  write(*,*) "Test CWC"

  call random_number(DENS)
  DENS(:,:,:) = DENS(:,:,:)*0.1 + 1.0_RP
  call random_number(MOMZ)
  call random_number(MOMX)
  call random_number(MOMY)

  RHOT(:,:,:) = POTT * DENS(:,:,:)
  QTRC(:,:,:,:) = Q

  REF_dens(:,:,:) = 1.0_RP
  REF_pott(:,:,:) = POTT
  REF_qv(:,:,:)   = Q
  REF_pres(:,:,:) = 1000._RP

  DAMP_var  (:,:,:,:) = -9.999E30_RP
  DAMP_alpha(:,:,:,:) = 0.0_RP

  MFLUX_OFFSET_X(:,:,:) = 0.0_RP
  MFLUX_OFFSET_Y(:,:,:) = 0.0_RP

  call COMM_vars8( DENS(:,:,:), 1 )
  call COMM_vars8( MOMZ(:,:,:), 2 )
  call COMM_vars8( MOMX(:,:,:), 3 )
  call COMM_vars8( MOMY(:,:,:), 4 )
  call COMM_vars8( RHOT(:,:,:), 5 )
  call COMM_wait ( DENS(:,:,:), 1 )
  call COMM_wait ( MOMZ(:,:,:), 2 )
  call COMM_wait ( MOMX(:,:,:), 3 )
  call COMM_wait ( MOMY(:,:,:), 4 )
  call COMM_wait ( RHOT(:,:,:), 5 )

  do iq = 1, QA
     call COMM_vars8( QTRC(:,:,:,iq), iq )
  enddo
  do iq = 1, QA
     call COMM_wait ( QTRC(:,:,:,iq), iq )
  enddo

  call ATMOS_DYN( &
       DENS, MOMZ, MOMX, MOMY, RHOT, QTRC,          & ! (inout)
       PROG,                                        & ! (inout)
       DENS_av, MOMZ_av, MOMX_av, MOMY_av, RHOT_av, QTRC_av, & ! (out)
       DENS_tp, MOMZ_tp, MOMX_tp, MOMY_tp, RHOT_tp, QTRC_tp, & ! (in)
       CORIOLIS,                                    & ! (in)
       CDZ, CDX, CDY, FDZ, FDX, FDY,                & ! (in)
       RCDZ, RCDX, RCDY, RFDZ, RFDX, RFDY,          & ! (in)
       PHI, GSQRT, J13G, J23G, J33G, MAPF,          & ! (in)
       AQ_R, AQ_CV, AQ_CP, AQ_MASS,                 & ! (in)
       REF_dens, REF_pott, REF_qv, REF_pres,        & ! (in)
       nd_coef, nd_coef, nd_order, nd_sfc_fact, nd_use_rs, & ! (in)
       BND_QA, BND_IQ, BND_SMOOTHER_FACT,                  & ! (in)
       DAMP_var(:,:,:,1), DAMP_var(:,:,:,2), DAMP_var(:,:,:,3), DAMP_var(:,:,:,4), DAMP_var(:,:,:,5), DAMP_var(:,:,:,6:6+QA-1), & ! (in)
       DAMP_alpha(:,:,:,1), DAMP_alpha(:,:,:,2), DAMP_alpha(:,:,:,3), DAMP_alpha(:,:,:,4), DAMP_alpha(:,:,:,5), & ! (in)
       DAMP_alpha(:,:,:,6:6+QA-1),                  & ! (in)
       MFLUX_OFFSET_X(:,:,:), MFLUX_OFFSET_Y(:,:,:), & ! (in)
       divdmp_coef,                                 & ! (in)
       flag_tracer_split_tend,                      & ! (in)
       flag_fct_momentum, flag_fct_t, flag_fct_tracer, & ! (in)
       flag_fct_along_stream,                       & ! (in)
       .false.,                                     & ! (in)
       I_QV,                                        & ! (in)
       1.0_DP, 1.0_DP                               ) ! (in)

  answer(:,:,:) = POTT
  call AssertEqual("POTT", answer(KS:KE,IS:IE,JS:JE), RHOT(KS:KE,IS:IE,JS:JE)/DENS(KS:KE,IS:IE,JS:JE), RP*2-2, -10)

  answer(:,:,:) = Q
  message = "iq = ??"
  do iq = 1, QA
     if ( .not. TRACER_ADVC(iq) ) cycle
     write(message(6:7), "(i2)") iq
     call AssertEqual(message, answer(KS:KE,IS:IE,JS:JE), QTRC(KS:KE,IS:IE,JS:JE,iq), RP*2-2, -10)
  end do

end subroutine test_cwc

subroutine test_fctminmax
  real(RP), parameter :: POTT_MAX = 350_RP
  real(RP), parameter :: POTT_MIN = 300_RP
  real(RP), parameter :: Q_MAX = 0.1_RP
  real(RP), parameter :: Q_MIN = 0_RP
  real(RP) :: MINMAX(KA,IA,JA)
  real(RP) :: epsilon
  integer :: k, i, j

  epsilon = 0.1_RP**(RP*2-2)

  write(*,*) "Test FCT"

  call random_number(DENS)
  DENS(:,:,:) = DENS(:,:,:)*0.1 + 1.0_RP
  call random_number(MOMZ)
  call random_number(MOMX)
  call random_number(MOMY)

  RHOT(:,:,:) = POTT_MIN * DENS(:,:,:)
  QTRC(:,:,:,:) = Q_MIN
  do j = JS, JS+JMAX/2
  do i = IS, IS+IMAX/2
  do k = KS, KS+KMAX/2
     RHOT(k,i,j) = POTT_MAX * DENS(k,i,j)
     QTRC(k,i,j,:) = Q_MAX
  end do
  end do
  end do

  REF_dens(:,:,:) = 1.0_RP
  REF_pott(:,:,:) = POTT_MIN
  REF_qv(:,:,:)   = Q_MIN
  REF_pres(:,:,:) = 1000._RP

  DAMP_var  (:,:,:,:) = -9.999E30_RP
  DAMP_alpha(:,:,:,:) = 0.0_RP

  MFLUX_OFFSET_X(:,:,:) = 0.0_RP
  MFLUX_OFFSET_Y(:,:,:) = 0.0_RP

  call COMM_vars8( DENS(:,:,:), 1 )
  call COMM_vars8( MOMZ(:,:,:), 2 )
  call COMM_vars8( MOMX(:,:,:), 3 )
  call COMM_vars8( MOMY(:,:,:), 4 )
  call COMM_vars8( RHOT(:,:,:), 5 )
  call COMM_wait ( DENS(:,:,:), 1 )
  call COMM_wait ( MOMZ(:,:,:), 2 )
  call COMM_wait ( MOMX(:,:,:), 3 )
  call COMM_wait ( MOMY(:,:,:), 4 )
  call COMM_wait ( RHOT(:,:,:), 5 )

  do iq = 1, QA
     call COMM_vars8( QTRC(:,:,:,iq), iq )
  enddo
  do iq = 1, QA
     call COMM_wait ( QTRC(:,:,:,iq), iq )
  enddo

  call ATMOS_DYN( &
       DENS, MOMZ, MOMX, MOMY, RHOT, QTRC,          & ! (inout)
       PROG,                                        & ! (inout)
       DENS_av, MOMZ_av, MOMX_av, MOMY_av, RHOT_av, QTRC_av, & ! (out)
       DENS_tp, MOMZ_tp, MOMX_tp, MOMY_tp, RHOT_tp, QTRC_tp, & ! (in)
       CORIOLIS,                                    & ! (in)
       CDZ, CDX, CDY, FDZ, FDX, FDY,                & ! (in)
       RCDZ, RCDX, RCDY, RFDZ, RFDX, RFDY,          & ! (in)
       PHI, GSQRT, J13G, J23G, J33G, MAPF,          & ! (in)
       AQ_R, AQ_CV, AQ_CP, AQ_MASS,                 & ! (in)
       REF_dens, REF_pott, REF_qv, REF_pres,        & ! (in)
       0.0_RP, 0.0_RP, nd_order, nd_sfc_fact, nd_use_rs, & ! (in)
       BND_QA, BND_IQ, BND_SMOOTHER_FACT,                & ! (in)
       DAMP_var(:,:,:,1), DAMP_var(:,:,:,2), DAMP_var(:,:,:,3), DAMP_var(:,:,:,4), DAMP_var(:,:,:,5), DAMP_var(:,:,:,6:6+QA-1), & ! (in)
       DAMP_alpha(:,:,:,1), DAMP_alpha(:,:,:,2), DAMP_alpha(:,:,:,3), DAMP_alpha(:,:,:,4), DAMP_alpha(:,:,:,5), & ! (in)
       DAMP_alpha(:,:,:,6:6+QA-1),                  & ! (in)
       MFLUX_OFFSET_X(:,:,:), MFLUX_OFFSET_Y(:,:,:), & ! (in)
       divdmp_coef,                                 & ! (in)
       flag_tracer_split_tend,                      & ! (in)
       flag_fct_momentum, flag_fct_t, flag_fct_tracer, & ! (in)
       flag_fct_along_stream,                       & ! (in)
       .false.,                                     & ! (in)
       I_QV,                                        & ! (in)
       1.0_DP, 1.0_DP                               ) ! (in)

  message = "iq = ??"
  do iq = 1, QA
     if ( .not. TRACER_ADVC(iq) ) cycle
     write(message(6:7), "(i2)") iq
     write(message(9:11),"(a3)") "MAX"
     MINMAX(:,:,:) = Q_MAX * ( 1_RP + epsilon )
     call AssertLessThan(message, MINMAX(KS:KE,IS:IE,JS:JE), QTRC(KS:KE,IS:IE,JS:JE,iq))
     write(message(9:11),"(a3)") "MIN"
     MINMAX(:,:,:) = Q_MIN - epsilon
     call AssertGreaterThan(message, MINMAX(KS:KE,IS:IE,JS:JE), QTRC(KS:KE,IS:IE,JS:JE,iq))
  end do

end subroutine test_fctminmax

! private
subroutine copy
  DENS_o(:,:,:) = DENS(:,:,:)
  MOMZ_o(:,:,:) = MOMZ(:,:,:)
  MOMX_o(:,:,:) = MOMX(:,:,:)
  MOMY_o(:,:,:) = MOMY(:,:,:)
  RHOT_o(:,:,:) = RHOT(:,:,:)
  QTRC_o(:,:,:,:) = QTRC(:,:,:,:)
end subroutine copy

end module test_atmos_dyn

max_stars_repo_name: slayoo/scale

max_issues_repo_name: slayoo/scale

max_forks_repo_name: slayoo/scale


text:
great_rsd(p1,ee1).
great_rsd(g1,dd1).
great_rsd(e1,x1).
great_rsd(e1,f1).
great_rsd(t1,m1).
great_rsd(u1,t1).

ext: f

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
      function ilog2(l)
      ilog2 = dlog(real(l, 8))/dlog(2.0d0)
      end

      subroutine rotate(theta, x1, x2, x1p, x2p)
      double precision theta, x1, x2, x1p, x2p, c, s
      c = cos(theta)
      s = sin(theta)
      x1p = c*x1-s*x2
      x2p = s*x1+c*x2
      end

      subroutine rot1(theta, n, x)
      double precision theta, x, x1p, x2p
      dimension x(n)
      do k = 1, n-1, 2
         call rotate(theta, x(k), x(k+1), x1p, x2p)
         x(k) = x1p
         x(k+1) = x2p
      end do
      end

      subroutine rot2(theta, n, x)
      double precision theta, x, x1p, x2p
      dimension x(n)
      do k = 2, n-2, 2
         call rotate(theta, x(k), x(k+1), x1p, x2p)
         x(k) = x1p
         x(k+1) = x2p
      end do
      end

      subroutine magsqr(n, x, y)
      double precision x, y
      dimension x(n)
      y = 0.0d0
      do k = 1, n
         y = y+x(k)*x(k)
      end do
      end

      subroutine f(n, x, l, phi, y)
      double precision phi, x, y, x1
      dimension x(n), x1(1000)
      do k = 1, n
         x1(k) = x(k)
      end do
      do i = 1, l
         m = 2**(ilog2(l)-
     +           ilog2(1+int(mod(1013.0d0*3.0d0**phi*real(i, 8),
     +                           real(l, 8)))))
         do j = 1, m
            call magsqr(n, x1, y)
            y = sqrt(y)
            call rot1(1.2d0*y, n, x1)
            call rot2(1.4d0*y, n, x1)
         end do
      end do
      call magsqr(n, x1, y)
      y = y/2.0d0
      end

      program main
      double precision phi, x, xb, y, yb
      dimension x(1000), xb(1000)
      read *, n
      read *, l
      read *, phi
      do k = 1, n
         x(k) = n+1-k
         xb(k) = 0.0d0
      end do
      yb = 1.0d0
      call f(n, x, l, phi, y)
      print *, y
      call f_b(n, x, xb, l, phi, y, yb)
      do k = 1, n
         print *, xb(k)
      end do
      end

max_stars_repo_name: qobi/checkpoint-VLAD

max_issues_repo_name: qobi/checkpoint-VLAD

max_forks_repo_name: qobi/checkpoint-VLAD


text:
!Crown Copyright 2012 AWE.
!
! This file is part of CloverLeaf.
!
! CloverLeaf is free software: you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the
! Free Software Foundation, either version 3 of the License, or (at your option)
! any later version.
!
! CloverLeaf is distributed in the hope that it will be useful, but
! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
! FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! CloverLeaf. If not, see http://www.gnu.org/licenses/.

!>  @brief Communication Utilities
!>  @author Wayne Gaudin
!>  @details Contains all utilities required to run CloverLeaf in a distributed
!>  environment, including initialisation, mesh decompostion, reductions and
!>  halo exchange using explicit buffers.
!>
!>  Note the halo exchange is currently coded as simply as possible and no
!>  optimisations have been implemented, such as post receives before sends or packing
!>  buffers with multiple data fields. This is intentional so the effect of these
!>  optimisations can be measured on large systems, as and when they are added.
!>
!>  Even without these modifications CloverLeaf weak scales well on moderately sized
!>  systems of the order of 10K cores.

MODULE clover_module

  USE data_module
  USE definitions_module

  IMPLICIT NONE

CONTAINS

SUBROUTINE clover_barrier

  INTEGER :: err

END SUBROUTINE clover_barrier

SUBROUTINE clover_abort

  INTEGER :: ierr,err

END SUBROUTINE clover_abort

SUBROUTINE clover_finalize

  INTEGER :: err

  CLOSE(g_out)
  CALL FLUSH(0)
  CALL FLUSH(6)
  CALL FLUSH(g_out)

END SUBROUTINE clover_finalize

SUBROUTINE clover_init_comms

  IMPLICIT NONE

  INTEGER :: err,rank,size

  rank=0
  size=1

  parallel%parallel=.TRUE.
  parallel%task=rank

  IF(rank.EQ.0) THEN
    parallel%boss=.TRUE.
  ENDIF

  parallel%boss_task=0
  parallel%max_task=size

END SUBROUTINE clover_init_comms

SUBROUTINE clover_get_num_chunks(count)

  IMPLICIT NONE

  INTEGER :: count

! Should be changed so there can be more than one chunk per mpi task

  count=parallel%max_task

END SUBROUTINE clover_get_num_chunks

SUBROUTINE clover_decompose(x_cells,y_cells,left,right,bottom,top)

  ! This decomposes the mesh into a number of chunks.
  ! The number of chunks may be a multiple of the number of mpi tasks
  ! Doesn't always return the best split if there are few factors
  ! All factors need to be stored and the best picked. But its ok for now

  IMPLICIT NONE

  INTEGER :: x_cells,y_cells,left(:),right(:),top(:),bottom(:)
  INTEGER :: c,delta_x,delta_y

  REAL(KIND=8) :: mesh_ratio,factor_x,factor_y
  INTEGER  :: chunk_x,chunk_y,mod_x,mod_y,split_found

  INTEGER  :: cx,cy,chunk,add_x,add_y,add_x_prev,add_y_prev

  ! 2D Decomposition of the mesh

  mesh_ratio=real(x_cells)/real(y_cells)

  chunk_x=number_of_chunks
  chunk_y=1

  split_found=0 ! Used to detect 1D decomposition
  DO c=1,number_of_chunks
    IF (MOD(number_of_chunks,c).EQ.0) THEN
      factor_x=number_of_chunks/real(c)
      factor_y=c
      !Compare the factor ratio with the mesh ratio
      IF(factor_x/factor_y.LE.mesh_ratio) THEN
        chunk_y=c
        chunk_x=number_of_chunks/c
        split_found=1
        EXIT
      ENDIF
    ENDIF
  ENDDO

  IF(split_found.EQ.0.OR.chunk_y.EQ.number_of_chunks) THEN ! Prime number or 1D decomp detected
    IF(mesh_ratio.GE.1.0) THEN
      chunk_x=number_of_chunks
      chunk_y=1
    ELSE
      chunk_x=1
      chunk_y=number_of_chunks
    ENDIF
  ENDIF

  delta_x=x_cells/chunk_x
  delta_y=y_cells/chunk_y
  mod_x=MOD(x_cells,chunk_x)
  mod_y=MOD(y_cells,chunk_y)

  ! Set up chunk mesh ranges and chunk connectivity

  add_x_prev=0
  add_y_prev=0
  chunk=1
  DO cy=1,chunk_y
    DO cx=1,chunk_x
      add_x=0
      add_y=0
      IF(cx.LE.mod_x)add_x=1
      IF(cy.LE.mod_y)add_y=1
      left(chunk)=(cx-1)*delta_x+1+add_x_prev
      right(chunk)=left(chunk)+delta_x-1+add_x
      bottom(chunk)=(cy-1)*delta_y+1+add_y_prev
      top(chunk)=bottom(chunk)+delta_y-1+add_y
      chunks(chunk)%chunk_neighbours(chunk_left)=chunk_x*(cy-1)+cx-1
      chunks(chunk)%chunk_neighbours(chunk_right)=chunk_x*(cy-1)+cx+1
      chunks(chunk)%chunk_neighbours(chunk_bottom)=chunk_x*(cy-2)+cx
      chunks(chunk)%chunk_neighbours(chunk_top)=chunk_x*(cy)+cx
      IF(cx.EQ.1)chunks(chunk)%chunk_neighbours(chunk_left)=external_face
      IF(cx.EQ.chunk_x)chunks(chunk)%chunk_neighbours(chunk_right)=external_face
      IF(cy.EQ.1)chunks(chunk)%chunk_neighbours(chunk_bottom)=external_face
      IF(cy.EQ.chunk_y)chunks(chunk)%chunk_neighbours(chunk_top)=external_face
      IF(cx.LE.mod_x)add_x_prev=add_x_prev+1
      chunk=chunk+1
    ENDDO
    add_x_prev=0
    IF(cy.LE.mod_y)add_y_prev=add_y_prev+1
  ENDDO

  IF(parallel%boss)THEN
    WRITE(g_out,*)
    WRITE(g_out,*)"Mesh ratio of ",mesh_ratio
    WRITE(g_out,*)"Decomposing the mesh into ",chunk_x," by ",chunk_y," chunks"
    WRITE(g_out,*)
  ENDIF

END SUBROUTINE clover_decompose

SUBROUTINE clover_allocate_buffers(chunk)

  IMPLICIT NONE

  INTEGER      :: chunk

  ! Unallocated buffers for external boundaries caused issues on some systems so they are now
  !  all allocated
  IF(parallel%task.EQ.chunks(chunk)%task)THEN
    !IF(chunks(chunk)%chunk_neighbours(chunk_left).NE.external_face) THEN
      ALLOCATE(chunks(chunk)%left_snd_buffer(2*(chunks(chunk)%field%y_max+5)))
      ALLOCATE(chunks(chunk)%left_rcv_buffer(2*(chunks(chunk)%field%y_max+5)))
    !ENDIF
    !IF(chunks(chunk)%chunk_neighbours(chunk_right).NE.external_face) THEN
      ALLOCATE(chunks(chunk)%right_snd_buffer(2*(chunks(chunk)%field%y_max+5)))
      ALLOCATE(chunks(chunk)%right_rcv_buffer(2*(chunks(chunk)%field%y_max+5)))
    !ENDIF
    !IF(chunks(chunk)%chunk_neighbours(chunk_bottom).NE.external_face) THEN
      ALLOCATE(chunks(chunk)%bottom_snd_buffer(2*(chunks(chunk)%field%x_max+5)))
      ALLOCATE(chunks(chunk)%bottom_rcv_buffer(2*(chunks(chunk)%field%x_max+5)))
    !ENDIF
    !IF(chunks(chunk)%chunk_neighbours(chunk_top).NE.external_face) THEN
      ALLOCATE(chunks(chunk)%top_snd_buffer(2*(chunks(chunk)%field%x_max+5)))
      ALLOCATE(chunks(chunk)%top_rcv_buffer(2*(chunks(chunk)%field%x_max+5)))
    !ENDIF
  ENDIF

END SUBROUTINE clover_allocate_buffers

SUBROUTINE clover_exchange(fields,depth)

  IMPLICIT NONE

  INTEGER      :: fields(:),depth

  ! Assuming 1 patch per task, this will be changed
  ! Also, not packing all fields for each communication, doing one at a time

  IF(fields(FIELD_DENSITY0).EQ.1) THEN
    CALL clover_exchange_message(parallel%task+1,chunks(parallel%task+1)%field%density0,      &
                                 chunks(parallel%task+1)%left_snd_buffer,                      &
                                 chunks(parallel%task+1)%left_rcv_buffer,                      &
                                 chunks(parallel%task+1)%right_snd_buffer,                     &
                                 chunks(parallel%task+1)%right_rcv_buffer,                     &
                                 chunks(parallel%task+1)%bottom_snd_buffer,                    &
                                 chunks(parallel%task+1)%bottom_rcv_buffer,                    &
                                 chunks(parallel%task+1)%top_snd_buffer,                       &
                                 chunks(parallel%task+1)%top_rcv_buffer,                       &
                                 depth,CELL_DATA)
  ENDIF

  IF(fields(FIELD_DENSITY1).EQ.1) THEN
    CALL clover_exchange_message(parallel%task+1,chunks(parallel%task+1)%field%density1,      &
                                 chunks(parallel%task+1)%left_snd_buffer,                      &
                                 chunks(parallel%task+1)%left_rcv_buffer,                      &
                                 chunks(parallel%task+1)%right_snd_buffer,                     &
                                 chunks(parallel%task+1)%right_rcv_buffer,                     &
                                 chunks(parallel%task+1)%bottom_snd_buffer,                    &
                                 chunks(parallel%task+1)%bottom_rcv_buffer,                    &
                                 chunks(parallel%task+1)%top_snd_buffer,                       &
                                 chunks(parallel%task+1)%top_rcv_buffer,                       &
                                 depth,CELL_DATA)
  ENDIF

  IF(fields(FIELD_ENERGY0).EQ.1) THEN
    CALL clover_exchange_message(parallel%task+1,chunks(parallel%task+1)%field%energy0,       &
                                 chunks(parallel%task+1)%left_snd_buffer,                      &
                                 chunks(parallel%task+1)%left_rcv_buffer,                      &
                                 chunks(parallel%task+1)%right_snd_buffer,                     &
                                 chunks(parallel%task+1)%right_rcv_buffer,                     &
                                 chunks(parallel%task+1)%bottom_snd_buffer,                    &
                                 chunks(parallel%task+1)%bottom_rcv_buffer,                    &
                                 chunks(parallel%task+1)%top_snd_buffer,                       &
                                 chunks(parallel%task+1)%top_rcv_buffer,                       &
                                 depth,CELL_DATA)
  ENDIF

  IF(fields(FIELD_ENERGY1).EQ.1) THEN
    CALL clover_exchange_message(parallel%task+1,chunks(parallel%task+1)%field%energy1,       &
                                 chunks(parallel%task+1)%left_snd_buffer,                      &
                                 chunks(parallel%task+1)%left_rcv_buffer,                      &
                                 chunks(parallel%task+1)%right_snd_buffer,                     &
                                 chunks(parallel%task+1)%right_rcv_buffer,                     &
                                 chunks(parallel%task+1)%bottom_snd_buffer,                    &
                                 chunks(parallel%task+1)%bottom_rcv_buffer,                    &
                                 chunks(parallel%task+1)%top_snd_buffer,                       &
                                 chunks(parallel%task+1)%top_rcv_buffer,                       &
                                 depth,CELL_DATA)
  ENDIF

  IF(fields(FIELD_PRESSURE).EQ.1) THEN
    CALL clover_exchange_message(parallel%task+1,chunks(parallel%task+1)%field%pressure,      &
                                 chunks(parallel%task+1)%left_snd_buffer,                      &
                                 chunks(parallel%task+1)%left_rcv_buffer,                      &
                                 chunks(parallel%task+1)%right_snd_buffer,                     &
                                 chunks(parallel%task+1)%right_rcv_buffer,                     &
                                 chunks(parallel%task+1)%bottom_snd_buffer,                    &
                                 chunks(parallel%task+1)%bottom_rcv_buffer,                    &
                                 chunks(parallel%task+1)%top_snd_buffer,                       &
                                 chunks(parallel%task+1)%top_rcv_buffer,                       &
                                 depth,CELL_DATA)
  ENDIF

  IF(fields(FIELD_VISCOSITY).EQ.1) THEN
    CALL clover_exchange_message(parallel%task+1,chunks(parallel%task+1)%field%viscosity,     &
                                 chunks(parallel%task+1)%left_snd_buffer,                      &
                                 chunks(parallel%task+1)%left_rcv_buffer,                      &
                                 chunks(parallel%task+1)%right_snd_buffer,                     &
                                 chunks(parallel%task+1)%right_rcv_buffer,                     &
                                 chunks(parallel%task+1)%bottom_snd_buffer,                    &
                                 chunks(parallel%task+1)%bottom_rcv_buffer,                    &
                                 chunks(parallel%task+1)%top_snd_buffer,                       &
                                 chunks(parallel%task+1)%top_rcv_buffer,                       &
                                 depth,CELL_DATA)
  ENDIF

  IF(fields(FIELD_SOUNDSPEED).EQ.1) THEN
    CALL clover_exchange_message(parallel%task+1,chunks(parallel%task+1)%field%soundspeed,    &
                                 chunks(parallel%task+1)%left_snd_buffer,                      &
                                 chunks(parallel%task+1)%left_rcv_buffer,                      &
                                 chunks(parallel%task+1)%right_snd_buffer,                     &
                                 chunks(parallel%task+1)%right_rcv_buffer,                     &
                                 chunks(parallel%task+1)%bottom_snd_buffer,                    &
                                 chunks(parallel%task+1)%bottom_rcv_buffer,                    &
                                 chunks(parallel%task+1)%top_snd_buffer,                       &
                                 chunks(parallel%task+1)%top_rcv_buffer,                       &
                                 depth,CELL_DATA)
  ENDIF

  IF(fields(FIELD_XVEL0).EQ.1) THEN
    CALL clover_exchange_message(parallel%task+1,chunks(parallel%task+1)%field%xvel0,         &
                                 chunks(parallel%task+1)%left_snd_buffer,                      &
                                 chunks(parallel%task+1)%left_rcv_buffer,                      &
                                 chunks(parallel%task+1)%right_snd_buffer,                     &
                                 chunks(parallel%task+1)%right_rcv_buffer,                     &
                                 chunks(parallel%task+1)%bottom_snd_buffer,                    &
                                 chunks(parallel%task+1)%bottom_rcv_buffer,                    &
                                 chunks(parallel%task+1)%top_snd_buffer,                       &
                                 chunks(parallel%task+1)%top_rcv_buffer,                       &
                                 depth,VERTEX_DATA)
  ENDIF

  IF(fields(FIELD_XVEL1).EQ.1) THEN
    CALL clover_exchange_message(parallel%task+1,chunks(parallel%task+1)%field%xvel1,         &
                                 chunks(parallel%task+1)%left_snd_buffer,                      &
                                 chunks(parallel%task+1)%left_rcv_buffer,                      &
                                 chunks(parallel%task+1)%right_snd_buffer,                     &
                                 chunks(parallel%task+1)%right_rcv_buffer,                     &
                                 chunks(parallel%task+1)%bottom_snd_buffer,                    &
                                 chunks(parallel%task+1)%bottom_rcv_buffer,                    &
                                 chunks(parallel%task+1)%top_snd_buffer,                       &
                                 chunks(parallel%task+1)%top_rcv_buffer,                       &
                                 depth,VERTEX_DATA)
  ENDIF

  IF(fields(FIELD_YVEL0).EQ.1) THEN
    CALL clover_exchange_message(parallel%task+1,chunks(parallel%task+1)%field%yvel0,         &
                                 chunks(parallel%task+1)%left_snd_buffer,                      &
                                 chunks(parallel%task+1)%left_rcv_buffer,                      &
                                 chunks(parallel%task+1)%right_snd_buffer,                     &
                                 chunks(parallel%task+1)%right_rcv_buffer,                     &
                                 chunks(parallel%task+1)%bottom_snd_buffer,                    &
                                 chunks(parallel%task+1)%bottom_rcv_buffer,                    &
                                 chunks(parallel%task+1)%top_snd_buffer,                       &
                                 chunks(parallel%task+1)%top_rcv_buffer,                       &
                                 depth,VERTEX_DATA)
  ENDIF

  IF(fields(FIELD_YVEL1).EQ.1) THEN
    CALL clover_exchange_message(parallel%task+1,chunks(parallel%task+1)%field%yvel1,         &
                                 chunks(parallel%task+1)%left_snd_buffer,                      &
                                 chunks(parallel%task+1)%left_rcv_buffer,                      &
                                 chunks(parallel%task+1)%right_snd_buffer,                     &
                                 chunks(parallel%task+1)%right_rcv_buffer,                     &
                                 chunks(parallel%task+1)%bottom_snd_buffer,                    &
                                 chunks(parallel%task+1)%bottom_rcv_buffer,                    &
                                 chunks(parallel%task+1)%top_snd_buffer,                       &
                                 chunks(parallel%task+1)%top_rcv_buffer,                       &
                                 depth,VERTEX_DATA)
  ENDIF

  IF(fields(FIELD_VOL_FLUX_X).EQ.1) THEN
    CALL clover_exchange_message(parallel%task+1,chunks(parallel%task+1)%field%vol_flux_x,    &
                                 chunks(parallel%task+1)%left_snd_buffer,                      &
                                 chunks(parallel%task+1)%left_rcv_buffer,                      &
                                 chunks(parallel%task+1)%right_snd_buffer,                     &
                                 chunks(parallel%task+1)%right_rcv_buffer,                     &
                                 chunks(parallel%task+1)%bottom_snd_buffer,                    &
                                 chunks(parallel%task+1)%bottom_rcv_buffer,                    &
                                 chunks(parallel%task+1)%top_snd_buffer,                       &
                                 chunks(parallel%task+1)%top_rcv_buffer,                       &
                                 depth,X_FACE_DATA)
  ENDIF

  IF(fields(FIELD_VOL_FLUX_Y).EQ.1) THEN
    CALL clover_exchange_message(parallel%task+1,chunks(parallel%task+1)%field%vol_flux_y,    &
                                 chunks(parallel%task+1)%left_snd_buffer,                      &
                                 chunks(parallel%task+1)%left_rcv_buffer,                      &
                                 chunks(parallel%task+1)%right_snd_buffer,                     &
                                 chunks(parallel%task+1)%right_rcv_buffer,                     &
                                 chunks(parallel%task+1)%bottom_snd_buffer,                    &
                                 chunks(parallel%task+1)%bottom_rcv_buffer,                    &
                                 chunks(parallel%task+1)%top_snd_buffer,                       &
                                 chunks(parallel%task+1)%top_rcv_buffer,                       &
                                 depth,Y_FACE_DATA)
  ENDIF

  IF(fields(FIELD_MASS_FLUX_X).EQ.1) THEN
    CALL clover_exchange_message(parallel%task+1,chunks(parallel%task+1)%field%mass_flux_x,   &
                                 chunks(parallel%task+1)%left_snd_buffer,                      &
                                 chunks(parallel%task+1)%left_rcv_buffer,                      &
                                 chunks(parallel%task+1)%right_snd_buffer,                     &
                                 chunks(parallel%task+1)%right_rcv_buffer,                     &
                                 chunks(parallel%task+1)%bottom_snd_buffer,                    &
                                 chunks(parallel%task+1)%bottom_rcv_buffer,                    &
                                 chunks(parallel%task+1)%top_snd_buffer,                       &
                                 chunks(parallel%task+1)%top_rcv_buffer,                       &
                                 depth,X_FACE_DATA)
  ENDIF

  IF(fields(FIELD_MASS_FLUX_Y).EQ.1) THEN
    CALL clover_exchange_message(parallel%task+1,chunks(parallel%task+1)%field%mass_flux_y,   &
                                 chunks(parallel%task+1)%left_snd_buffer,                      &
                                 chunks(parallel%task+1)%left_rcv_buffer,                      &
                                 chunks(parallel%task+1)%right_snd_buffer,                     &
                                 chunks(parallel%task+1)%right_rcv_buffer,                     &
                                 chunks(parallel%task+1)%bottom_snd_buffer,                    &
                                 chunks(parallel%task+1)%bottom_rcv_buffer,                    &
                                 chunks(parallel%task+1)%top_snd_buffer,                       &
                                 chunks(parallel%task+1)%top_rcv_buffer,                       &
                                 depth,Y_FACE_DATA)
  ENDIF

END SUBROUTINE clover_exchange

SUBROUTINE clover_exchange_message(chunk,field,                            &
                                   left_snd_buffer,                        &
                                   left_rcv_buffer,                        &
                                   right_snd_buffer,                       &
                                   right_rcv_buffer,                       &
                                   bottom_snd_buffer,                      &
                                   bottom_rcv_buffer,                      &
                                   top_snd_buffer,                         &
                                   top_rcv_buffer,                         &
                                   depth,field_type)

  USE pack_kernel_module

  IMPLICIT NONE

  REAL(KIND=8) :: field(-1:,-1:) ! This seems to work for any type of mesh data
  REAL(KIND=8) :: left_snd_buffer(:),left_rcv_buffer(:),right_snd_buffer(:),right_rcv_buffer(:)
  REAL(KIND=8) :: bottom_snd_buffer(:),bottom_rcv_buffer(:),top_snd_buffer(:),top_rcv_buffer(:)

  INTEGER      :: chunk,depth,field_type

  INTEGER      :: size,err,request(8),tag,message_count,x_inc,y_inc
  INTEGER      :: receiver,sender

  ! Field type will either be cell, vertex, x_face or y_face to get the message limits correct

  !  but this will do for now

  ! I am also sending buffers to chunks with the same task id for now.
  ! This can be improved in the future but at the moment there is just 1 chunk per task anyway

  ! The tag will be a function of the sending chunk and the face it is coming from
  !  like chunk 6 sending the left face

  ! No open mp in here either. May be beneficial will packing and unpacking in the future, though I am not sure.

  ! Change this so it will allow more than 1 chunk per task

  request=0
  message_count=0

  ! Pack and send

  ! These array modifications still need to be added on, plus the donor data location changes as in update_halo
  IF(field_type.EQ.CELL_DATA) THEN
    x_inc=0
    y_inc=0
  ENDIF
  IF(field_type.EQ.VERTEX_DATA) THEN
    x_inc=1
    y_inc=1
  ENDIF
  IF(field_type.EQ.X_FACE_DATA) THEN
    x_inc=1
    y_inc=0
  ENDIF
  IF(field_type.EQ.Y_FACE_DATA) THEN
    x_inc=0
    y_inc=1
  ENDIF

  ! Pack real data into buffers
  IF(parallel%task.EQ.chunks(chunk)%task) THEN
    size=(1+(chunks(chunk)%field%y_max+y_inc+depth)-(chunks(chunk)%field%y_min-depth))*depth
    IF(use_fortran_kernels) THEN
      CALL pack_left_right_buffers(chunks(chunk)%field%x_min,chunks(chunk)%field%x_max, &
                                   chunks(chunk)%field%y_min,chunks(chunk)%field%y_max, &
                                   chunks(chunk)%chunk_neighbours(chunk_left),          &
                                   chunks(chunk)%chunk_neighbours(chunk_right),         &
                                   external_face,                                       &
                                   x_inc,y_inc,depth,size,                              &
                                   field,left_snd_buffer,right_snd_buffer)
    ELSEIF(use_C_kernels)THEN
      CALL pack_left_right_buffers_c(chunks(chunk)%field%x_min,chunks(chunk)%field%x_max, &
                                     chunks(chunk)%field%y_min,chunks(chunk)%field%y_max, &
                                     chunks(chunk)%chunk_neighbours(chunk_left),          &
                                     chunks(chunk)%chunk_neighbours(chunk_right),         &
                                     external_face,                                       &
                                     x_inc,y_inc,depth,size,                              &
                                     field,left_snd_buffer,right_snd_buffer)
    ENDIF

    ! Send/receive the data
    IF(chunks(chunk)%chunk_neighbours(chunk_left).NE.external_face) THEN
      tag=4*(chunk)+1 ! 4 because we have 4 faces, 1 because it is leaving the left face
      receiver=chunks(chunks(chunk)%chunk_neighbours(chunk_left))%task
      tag=4*(chunks(chunk)%chunk_neighbours(chunk_left))+2 ! 4 because we have 4 faces, 1 because it is coming from the right face of the left neighbour
      sender=chunks(chunks(chunk)%chunk_neighbours(chunk_left))%task
      message_count=message_count+2
    ENDIF

    IF(chunks(chunk)%chunk_neighbours(chunk_right).NE.external_face) THEN
      tag=4*chunk+2 ! 4 because we have 4 faces, 2 because it is leaving the right face
      receiver=chunks(chunks(chunk)%chunk_neighbours(chunk_right))%task
      tag=4*(chunks(chunk)%chunk_neighbours(chunk_right))+1 ! 4 because we have 4 faces, 1 because it is coming from the left face of the right neighbour
      sender=chunks(chunks(chunk)%chunk_neighbours(chunk_right))%task
      message_count=message_count+2
    ENDIF
  ENDIF

  ! Wait for the messages

  ! Unpack buffers in halo cells
  IF(parallel%task.EQ.chunks(chunk)%task) THEN
    IF(use_fortran_kernels) THEN
      CALL unpack_left_right_buffers(chunks(chunk)%field%x_min,chunks(chunk)%field%x_max, &
                                     chunks(chunk)%field%y_min,chunks(chunk)%field%y_max, &
                                     chunks(chunk)%chunk_neighbours(chunk_left),          &
                                     chunks(chunk)%chunk_neighbours(chunk_right),         &
                                     external_face,                                       &
                                     x_inc,y_inc,depth,size,                              &
                                     field,left_rcv_buffer,right_rcv_buffer)
    ELSEIF(use_C_kernels)THEN
      CALL unpack_left_right_buffers_c(chunks(chunk)%field%x_min,chunks(chunk)%field%x_max, &
                                       chunks(chunk)%field%y_min,chunks(chunk)%field%y_max, &
                                       chunks(chunk)%chunk_neighbours(chunk_left),          &
                                       chunks(chunk)%chunk_neighbours(chunk_right),         &
                                       external_face,                                       &
                                       x_inc,y_inc,depth,size,                              &
                                       field,left_rcv_buffer,right_rcv_buffer)
    ENDIF
  ENDIF

  request=0
  message_count=0

  ! Pack real data into buffers
  IF(parallel%task.EQ.chunks(chunk)%task) THEN
    size=(1+(chunks(chunk)%field%x_max+x_inc+depth)-(chunks(chunk)%field%x_min-depth))*depth
    IF(use_fortran_kernels) THEN
      CALL pack_top_bottom_buffers(chunks(chunk)%field%x_min,chunks(chunk)%field%x_max, &
                                   chunks(chunk)%field%y_min,chunks(chunk)%field%y_max, &
                                   chunks(chunk)%chunk_neighbours(chunk_bottom),        &
                                   chunks(chunk)%chunk_neighbours(chunk_top),           &
                                   external_face,                                       &
                                   x_inc,y_inc,depth,size,                              &
                                   field,bottom_snd_buffer,top_snd_buffer)
    ELSEIF(use_C_kernels)THEN
      CALL pack_top_bottom_buffers_c(chunks(chunk)%field%x_min,chunks(chunk)%field%x_max, &
                                     chunks(chunk)%field%y_min,chunks(chunk)%field%y_max, &
                                     chunks(chunk)%chunk_neighbours(chunk_bottom),        &
                                     chunks(chunk)%chunk_neighbours(chunk_top),           &
                                     external_face,                                       &
                                     x_inc,y_inc,depth,size,                              &
                                     field,bottom_snd_buffer,top_snd_buffer)
    ENDIF

    ! Send/receive the data
    IF(chunks(chunk)%chunk_neighbours(chunk_bottom).NE.external_face) THEN
      tag=4*(chunk)+3 ! 4 because we have 4 faces, 3 because it is leaving the bottom face
      receiver=chunks(chunks(chunk)%chunk_neighbours(chunk_bottom))%task
      tag=4*(chunks(chunk)%chunk_neighbours(chunk_bottom))+4 ! 4 because we have 4 faces, 1 because it is coming from the top face of the bottom neighbour
      sender=chunks(chunks(chunk)%chunk_neighbours(chunk_bottom))%task
      message_count=message_count+2
    ENDIF

    IF(chunks(chunk)%chunk_neighbours(chunk_top).NE.external_face) THEN
      tag=4*(chunk)+4 ! 4 because we have 4 faces, 4 because it is leaving the top face
      receiver=chunks(chunks(chunk)%chunk_neighbours(chunk_top))%task
      tag=4*(chunks(chunk)%chunk_neighbours(chunk_top))+3 ! 4 because we have 4 faces, 4 because it is coming from the left face of the top neighbour
      sender=chunks(chunks(chunk)%chunk_neighbours(chunk_top))%task
      message_count=message_count+2
    ENDIF

  ENDIF

  ! Wait for the messages

  ! Unpack buffers in halo cells
  IF(parallel%task.EQ.chunks(chunk)%task) THEN
    IF(use_fortran_kernels) THEN
      CALL unpack_top_bottom_buffers(chunks(chunk)%field%x_min,chunks(chunk)%field%x_max, &
                                     chunks(chunk)%field%y_min,chunks(chunk)%field%y_max, &
                                     chunks(chunk)%chunk_neighbours(chunk_bottom),        &
                                     chunks(chunk)%chunk_neighbours(chunk_top),           &
                                     external_face,                                       &
                                     x_inc,y_inc,depth,size,                              &
                                     field,bottom_rcv_buffer,top_rcv_buffer)
    ELSEIF(use_C_kernels)THEN
      CALL unpack_top_bottom_buffers_c(chunks(chunk)%field%x_min,chunks(chunk)%field%x_max, &
                                       chunks(chunk)%field%y_min,chunks(chunk)%field%y_max, &
                                       chunks(chunk)%chunk_neighbours(chunk_bottom),        &
                                       chunks(chunk)%chunk_neighbours(chunk_top),           &
                                       external_face,                                       &
                                       x_inc,y_inc,depth,size,                              &
                                       field,bottom_rcv_buffer,top_rcv_buffer)
    ENDIF
  ENDIF

END SUBROUTINE clover_exchange_message

SUBROUTINE clover_sum(value)

  ! Only sums to the master

  IMPLICIT NONE

  REAL(KIND=8) :: value

  REAL(KIND=8) :: total

  INTEGER :: err

  total=value

  value=total

END SUBROUTINE clover_sum

SUBROUTINE clover_min(value)

  IMPLICIT NONE

  REAL(KIND=8) :: value

  REAL(KIND=8) :: minimum

  INTEGER :: err

  minimum=value

  value=minimum

END SUBROUTINE clover_min

SUBROUTINE clover_max(value)

  IMPLICIT NONE

  REAL(KIND=8) :: value

  REAL(KIND=8) :: maximum

  INTEGER :: err

  maximum=value

  value=maximum

END SUBROUTINE clover_max

SUBROUTINE clover_allgather(value,values)

  IMPLICIT NONE

  REAL(KIND=8) :: value

  REAL(KIND=8) :: values(parallel%max_task)

  INTEGER :: err

  values(1)=value ! Just to ensure it will work in serial

END SUBROUTINE clover_allgather

SUBROUTINE clover_check_error(error)

  IMPLICIT NONE

  INTEGER :: error

  INTEGER :: maximum

  INTEGER :: err

  maximum=error

  error=maximum

END SUBROUTINE clover_check_error

END MODULE clover_module

max_stars_repo_name: uwplse/stng

max_issues_repo_name: uwplse/stng

max_forks_repo_name: uwplse/stng


text:
program test

    use omp_lib

    use quadrature2, only: simpson, error_table
    use functions, only: f, fevals, k

    implicit none
    real(kind=8) :: a,b,int_true
    integer :: nvals(12)
    integer :: i, nthreads

    real(kind=8) :: t1, t2, elapsed_time
    integer(kind=8) :: tclock1, tclock2, clock_rate

    nthreads = 1      ! for serial mode
    !$ nthreads = 4   ! for openmp
    !$ call omp_set_num_threads(nthreads)
    print 100, nthreads
100 format("Using ",i2," threads")

    fevals = 0

    k = 1.d3   ! functions module variable for function f2
    a = 0.d0
    b = 2.d0
    int_true = (b-a) + (b**4 - a**4) / 4.d0 - (1.d0/k) * (cos(k*b) - cos(k*a))

    print 10, int_true
 10 format("true integral: ", es22.14)
    print *, " "  ! blank line

    ! values of n to test:   (larger values than before)
    do i=1,12
        nvals(i) = 50 * 2**(i-1)
        enddo

    ! time the call to error_table:
    call system_clock(tclock1)
    call cpu_time(t1)
    call error_table(f, a, b, nvals, int_true, simpson)
    call cpu_time(t2)
    call system_clock(tclock2, clock_rate)

    elapsed_time = float(tclock2 - tclock1) / float(clock_rate)
    print *, " "
    print 11, elapsed_time
 11 format("Elapsed time = ",f12.8, " seconds")

    print 12, t2-t1
 12 format("CPU time = ",f12.8, " seconds")

    ! print the number of function evaluations by each thread:
    do i=0,nthreads-1
        print 101,  i, fevals(i)
101     format("fevals by thread ",i2,": ",i13)
        enddo

    print 102, sum(fevals)
102 format("Total number of fevals: ",i10)

end program test

max_stars_repo_name: philipwangdk/HPC

max_issues_repo_name: philipwangdk/HPC

max_forks_repo_name: philipwangdk/HPC


text:
program complex_05
implicit none

complex, parameter :: x = (1.0,-3.0)
real, parameter :: a = 3.0, b = 4.0
complex, parameter :: y = (a, b)

complex, parameter :: i_ = (0, 1)
complex, parameter :: z = a + i_*b
complex, parameter :: w = a+b + i_*(a-b)

print *, x, y, z, w

end program

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


! test_mod_ls65_OMP.f90   2018/10/27  2017/9/19
!
! ifort -qopenmp -mkl -static -O -o test_mod_ls65_OMPstatic.ex mod_mt19937ar.f90 mod_least65_OMP.f90 test_mod_ls65_OMP.f90
! export OMP_NUM_THREADS=6
! time ./test_mod_ls65_OMP.ex < tsls65_OMP.dat > tsls65_OMP.out6
!
! ifort -qopenmp -mkl -O -o test_mod_ls65_OMP.ex mod_mt19937ar.f90 mod_least65_OMP.f90 test_mod_ls65_OMP.f90
! gfortran -fopenmp -Wall -O -o test_mod_ls65_OMP.ex mod_mt19937ar.f90 mod_least65_OMP.f90 test_mod_ls65_OMP.f90 -llapack -lblas
!
!-----
MODULE least_var
!
  IMPLICIT NONE
  integer, parameter :: DP = kind(1.0D0)
!
  integer(4), parameter :: set_NSMAX=2048
  real(DP), save :: Q(set_NSMAX)
!
END MODULE least_var
!
!-----
program main
!
!$ use omp_lib
  use mod_least, only : NPR,NS,NP,NPFIT,IPFIT,PA,IPRSW,IDRSW,OBS,WSQRT,LEAST,least_alloc,least_dealloc &
  & ,copy_to_init, copy_from_init, copy_to_init_array, copy_from_init_array
  use least_var, only : Q,set_NSMAX
!
  IMPLICIT NONE
  integer, parameter :: DP = kind(1.0D0)
!
  integer(4) :: MAXX,ICON
  real(DP) :: EPSR,SUMS,SIGST
  integer(4) :: i,k,kk,k_begin,k_end
  integer(4) :: imax_init
  integer(4) :: NPR_init,NS_init,NP_init,NPFIT_init,IPRSW_init,IDRSW_init,MAXX_init
  integer(4), ALLOCATABLE :: IPFIT_init(:)
  real(DP) :: EPSR_init
  real(DP), ALLOCATABLE :: PA_init(:,:)
  real(DP), ALLOCATABLE :: OBS_init(:), WSQRT_init(:)
  integer(4), ALLOCATABLE :: ICON_fin(:)
  real(DP), ALLOCATABLE :: PA_fin(:,:),SUMS_fin(:)
  integer(4), parameter :: parallel_do_unit=60, chunk=1   ! parallel_do_unit=48
!  parallel_do_unit=processing unit of  !!$OMP parallel do  e.g. 10 * OMP_NUM_THREADS
!
!$ call omp_set_dynamic(.false.)
!$ write(*,'(A)')'  omp_set_dynamic(.false.)'
!$ write(*,'(A,I10)')'  omp_get_max_threads() = ',omp_get_max_threads() ! export OMP_NUM_THREADS=6
!$ write(*,'(A,I10)')'  omp_get_num_procs() = ',omp_get_num_procs()
!$ write(*,'(A,L10)')'  omp_get_nested() = ',omp_get_nested()           ! omp_get_nested() = F
!$ write(*,'(A,L10)')'  omp_get_dynamic() = ',omp_get_dynamic()         ! omp_get_dynamic() = F
!
  call init_rand_MT    !  Initialize Mersenne-Twister pseudo-random number generator
  READ(*,*) imax_init  !  number of initial parameter set
  WRITE(*,'(A,I10)')'  imax_init = ',imax_init
  if(imax_init < 1) then
    write(*,*)' invalid imax_init'  ;  stop
  endif
!
  NPR = 6    ! write(NPR,'(A)')'  write output'
  IDRSW = 1  ! NUMERICAL DERIVATIVE
!
  READ(*,*) NS, NP, NPFIT
  WRITE(*,'(a,3I5)')'  NS, NP, NPFIT=', NS, NP, NPFIT
  if(NS > set_NSMAX) then
    write(*,*)' invalid input'  ;  stop
  endif
!
  call least_alloc  ! allocate IPFIT(1:NPFIT),OBS(1:NS),WSQRT(1:NS) etc
!
  READ(*,*) IPFIT(1:NPFIT)
  WRITE(*,'(a)')'  IPFIT'
  WRITE(*,'(5I3,2X,5I3)') IPFIT(1:NPFIT)
  READ(*,*) EPSR, MAXX, IPRSW
  WRITE(*,'(a,g15.7,I10)')'  EPSR,MAXX=',EPSR,MAXX
  WRITE(*,'(a,I10)')'  IPRSW=',IPRSW
  WRITE(*,'(a)')'  IDRSW=1   NUMERICAL DERIVATIVE'
  do i=1,NS
    READ(*,*)Q(i),OBS(i),WSQRT(i)
  end do
!
!      copy NPR     ,IDRSW     ,NS     ,NP     ,NPFIT     ,EPSR     ,MAXX     ,IPRSW
!        to NPR_init,IDRSW_init,NS_init,NP_init,NPFIT_init,EPSR_init,MAXX_init,IPRSW_init
  call copy_to_init( &
   &  NPR     ,IDRSW     ,NS     ,NP     ,NPFIT     ,EPSR     ,MAXX     ,IPRSW      &
   & ,NPR_init,IDRSW_init,NS_init,NP_init,NPFIT_init,EPSR_init,MAXX_init,IPRSW_init )
!
  ALLOCATE( IPFIT_init(NPFIT) )
  ALLOCATE(   OBS_init(NS)    )
  ALLOCATE( WSQRT_init(NS)    )
  ALLOCATE( PA_init(NP, parallel_do_unit) )
  ALLOCATE(  PA_fin(NP, parallel_do_unit) )
  ALLOCATE( ICON_fin(parallel_do_unit)    )
  ALLOCATE( SUMS_fin(parallel_do_unit)    )
!
!      copy IPFIT     ,OBS     ,WSQRT
!        to IPFIT_init,OBS_init,WSQRT_init
  call copy_to_init_array( NPFIT, NS       &
       & ,IPFIT     ,OBS     ,WSQRT        &
       & ,IPFIT_init,OBS_init,WSQRT_init   )
!
  call least_dealloc  ! deallocate IPFIT(1:NPFIT),OBS(1:NS),WSQRT(1:NS) etc
!
  write(*,'(A)')'  k  ICON  SUMS  PA(1:NP)'
  do k_begin = 1, imax_init, parallel_do_unit
    k_end = MIN(k_begin + parallel_do_unit -1, imax_init)
    kk = 1 + k_end - k_begin
           ! set initial PA(1:NP) and save in PA_init(1:NP, 1:kk )
    call set_PA_init(NP,PA_init,kk)
!
!$OMP parallel do default(none) &
!$OMP private(k,kk, EPSR,MAXX,SUMS,SIGST,ICON) &
!$OMP shared(k_begin,k_end, NPR_init,IDRSW_init,NS_init,NP_init,NPFIT_init,IPFIT_init) &
!$OMP shared(EPSR_init,MAXX_init,IPRSW_init,OBS_init,WSQRT_init,PA_init) &
!$OMP shared(ICON_fin,SUMS_fin,PA_fin) &
!$OMP schedule(dynamic,chunk)
    do k = k_begin, k_end
      kk= 1 + k - k_begin
!          copy back NPR     ,IDRSW     ,NS     ,NP     ,NPFIT     ,EPSR     ,MAXX     ,IPRSW
!               from NPR_init,IDRSW_init,NS_init,NP_init,NPFIT_init,EPSR_init,MAXX_init,IPRSW_init
      call copy_from_init( &
   &  NPR     ,IDRSW     ,NS     ,NP     ,NPFIT     ,EPSR     ,MAXX     ,IPRSW      &
   & ,NPR_init,IDRSW_init,NS_init,NP_init,NPFIT_init,EPSR_init,MAXX_init,IPRSW_init )
!
      call least_alloc    ! allocate threadprivate IPFIT(1:NPFIT),OBS(1:NS),WSQRT(1:NS) etc
!
!          copy back IPFIT     ,OBS     ,WSQRT
!               from IPFIT_init,OBS_init,WSQRT_init
      call copy_from_init_array( NPFIT, NS       &
             & ,IPFIT     ,OBS     ,WSQRT        &
             & ,IPFIT_init,OBS_init,WSQRT_init   )
!          copy back PA from PA_init
      PA(1:NP) = PA_init(1:NP, kk)
!
      CALL LEAST(EPSR,MAXX,SUMS,SIGST,ICON)
      ICON_fin(kk) = ICON         ! keep ICON in ICON_fin
      SUMS_fin(kk) = SUMS         ! keep SUMS in SUMS_fin
      PA_fin(1:NP,kk) = PA(1:NP)  ! keep PA   in PA_fin
!
      call least_dealloc    ! deallocate threadprivate IPFIT(1:NPFIT),OBS(1:NS),WSQRT(1:NS) etc
    end do
!$OMP end parallel do
!
    do k = k_begin, k_end      ! print results ICON_fin, SUMS_fin, PA_fin
      kk= 1 + k - k_begin
      write(*,'(2I5,100G15.7)') k, ICON_fin(kk), SUMS_fin(kk), PA_fin(1:NP_init,kk)
    end do
!
  end do
!
  DEALLOCATE( PA_init    )
  DEALLOCATE( ICON_fin   )
  DEALLOCATE( SUMS_fin   )
  DEALLOCATE( PA_fin     )
  DEALLOCATE( IPFIT_init )
  DEALLOCATE( OBS_init   )
  DEALLOCATE( WSQRT_init )
!
end program main
!-----
subroutine set_PA_init(NP,PA_init,kk)
!
  use mod_mt19937ar , only : genrand_res53
  IMPLICIT NONE
  integer, parameter :: DP = kind(1.0D0)
!
  integer(4),intent(in ) :: NP,kk
  real(DP),intent(out) :: PA_init(NP,kk)
  integer(4) :: i,j
!
  do i = 1, kk
    PA_init(1, i)= 1.0d0 * (1.0d0 + ( -1.0d0 + 2.0d0 * genrand_res53() ) * 0.1d0 )
    PA_init(2, i)= 1.0d0 * (1.0d0 + ( -1.0d0 + 2.0d0 * genrand_res53() ) * 0.1d0 )
    do j=3,10000,6
      if( j > NP ) exit
      PA_init(j  , i)=  0.25d0 * (1.0d0 + ( -1.0d0 + 2.0d0 * genrand_res53() ) * 0.3d0 )
      PA_init(j+1, i)= -0.90d0 + ((j-3)/6)*0.2d0 + ( -1.0d0 + 2.0d0 * genrand_res53() ) * 0.03d0
      PA_init(j+2, i)=  0.02d0 * (1.0d0 + ( -1.0d0 + 2.0d0 * genrand_res53() ) * 0.3d0 )
      PA_init(j+3, i)= -0.25d0 * (1.0d0 + ( -1.0d0 + 2.0d0 * genrand_res53() ) * 0.3d0 )
      PA_init(j+4, i)= -0.80d0 + ((j-3)/6)*0.2d0 + ( -1.0d0 + 2.0d0 * genrand_res53() ) * 0.03d0
      PA_init(j+5, i)=  0.02d0 * (1.0d0 + ( -1.0d0 + 2.0d0 * genrand_res53() ) * 0.3d0 )
    end do
  end do
!
end subroutine set_PA_init
!
!-----
subroutine FMODL
!
  use mod_least, only : PA,CALC,NS
  use least_var, only : Q
!
  IMPLICIT NONE
  integer, parameter :: DP = kind(1.0D0)
  real(DP) :: xx
  integer(4) :: I,j
!
      DO I=1,NS
        xx=PA(1)+PA(2)*Q(I)
        do j=0,18
         xx=xx+PA(3+j*3)/( ( ( Q(I)-PA(4+j*3) )/PA(5+j*3) )**2 + 1.0d0 )
        end do
        CALC(I)=xx
      end do
!
end subroutine FMODL
!
!-----
subroutine DFMODL
  IMPLICIT NONE
!
  continue
  return
!
end subroutine DFMODL
!
!-----
!  Initialize mod_mt19937ar
!             Mersenne Twister pseudo-random number generator
!    see http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/mt19937ar.html
!        https://en.wikipedia.org/wiki/Mersenne_Twister
subroutine init_rand_MT
!
  use mod_mt19937ar , only : init_by_array, genrand_res53
!
  implicit none
  integer, parameter :: DP = kind(1.0D0)
!
  real(DP) :: xx
  integer(4) :: j
  integer :: init(4)=(/ 291,564,837,1110 /)    ! default of initial state of mod_mt19937ar
  integer :: length=4
!
  read(*,*) init(1:4)                          ! initial state of mod_mt19937ar
  write(*,'(/,A)')'  call init_by_array(init,length)'
  write(*,'(A,4I15)')'     init(1:4) = ',init(1:4)
  write(*,'(A,I10)') '     length = ',length
  call init_by_array(init,length)              ! set initial state of mod_mt19937ar
  write(*,'(/,A)')'  four  genrand_res53() = (0,1) with 53-bit resolution'
  do j=1,4
    xx=genrand_res53()
    write(*,'(10X,G15.7)')xx                   ! first 4 random numbers of genrand_res53()
  end do
!
end subroutine init_rand_MT
!


text:
C Output from Public domain Ratfor, version 1.05
      subroutine rqfnb(n,p,a,y,rhs,d,u,beta,eps,wn,wp,nit,info)
      integer n,p,info,nit(3)
      double precision a(p,n),y(n),rhs(p),d(n),u(n),wn(n,9),wp(p,p+3)
      double precision one,beta,eps
      parameter( one = 1.0d0)
      call lpfnb(n,p,a,y,rhs,d,u,beta,eps,wn(1,1),wn(1,2), wp(1,1),wn(1,
     *3),wn(1,4),wn(1,5),wn(1,6), wp(1,2),wn(1,7),wn(1,8),wn(1,9),wp(1,3
     *),wp(1,4),nit,info)
      return
      end
      subroutine lpfnb(n,p,a,c,b,d,u,beta,eps,x,s,y,z,w, dx,ds,dy,dz,dw,
     *dr,rhs,ada,nit,info)
      integer n,p,pp,i,info,nit(3),maxit
      double precision a(p,n),c(n),b(p)
      double precision zero,one,mone,big,ddot,dmax1,dmin1,dxdz,dsdw
      double precision deltap,deltad,beta,eps,mu,gap,g
      double precision x(n),u(n),s(n),y(p),z(n),w(n),d(n),rhs(p),ada(p,p
     *)
      double precision dx(n),ds(n),dy(p),dz(n),dw(n),dr(n)
      parameter( zero = 0.0d0)
      parameter( one = 1.0d0)
      parameter( mone = -1.0d0)
      parameter( big = 1.0d+20)
      parameter( maxit = 500)
      nit(1)=0
      nit(2)=0
      nit(3)=n
      pp=p*p
      call dgemv('N',p,n,one,a,p,c,1,zero,y,1)
      do23000 i=1,n
      d(i)=one
23000 continue
23001 continue
      call stepy(n,p,a,d,y,ada,info)
      if(info .ne. 0)then
      return
      endif
      call dcopy(n,c,1,s,1)
      call dgemv('T',p,n,mone,a,p,y,1,one,s,1)
      do23004 i=1,n
      if(dabs(s(i)).lt.eps)then
      z(i)=dmax1(s(i), zero) + eps
      w(i)=dmax1(-s(i),zero) + eps
      else
      z(i)=dmax1(s(i), zero)
      w(i)=dmax1(-s(i),zero)
      endif
      s(i)=u(i)-x(i)
23004 continue
23005 continue
      gap = ddot(n,z,1,x,1)+ddot(n,w,1,s,1)
23008 if(gap .gt. eps .and. nit(1).lt.maxit)then
      nit(1)=nit(1)+1
      do23010 i = 1,n
      d(i) = one/(z(i)/x(i) + w(i)/s(i))
      ds(i)=z(i)-w(i)
      dz(i)=d(i)*ds(i)
23010 continue
23011 continue
      call dcopy(p,b,1,dy,1)
      call dgemv('N',p,n,mone,a,p,x,1,one,dy,1)
      call dgemv('N',p,n,one,a,p,dz,1,one,dy,1)
      call dcopy(p,dy,1,rhs,1)
      call stepy(n,p,a,d,dy,ada,info)
      if(info .ne. 0)then
      return
      endif
      call dgemv('T',p,n,one,a,p,dy,1,mone,ds,1)
      deltap=big
      deltad=big
      do23014 i=1,n
      dx(i)=d(i)*ds(i)
      ds(i)=-dx(i)
      dz(i)=-z(i)*(dx(i)/x(i) + one)
      dw(i)=-w(i)*(ds(i)/s(i) + one)
      if(dx(i).lt.0)then
      deltap=dmin1(deltap,-x(i)/dx(i))
      endif
      if(ds(i).lt.0)then
      deltap=dmin1(deltap,-s(i)/ds(i))
      endif
      if(dz(i).lt.0)then
      deltad=dmin1(deltad,-z(i)/dz(i))
      endif
      if(dw(i).lt.0)then
      deltad=dmin1(deltad,-w(i)/dw(i))
      endif
23014 continue
23015 continue
      deltap=dmin1(beta*deltap,one)
      deltad=dmin1(beta*deltad,one)
      if(min(deltap,deltad) .lt. one)then
      nit(2)=nit(2)+1
      mu = ddot(n,x,1,z,1)+ddot(n,s,1,w,1)
      g = mu + deltap*ddot(n,dx,1,z,1)+ deltad*ddot(n,dz,1,x,1) + deltap
     **deltad*ddot(n,dz,1,dx,1)+ deltap*ddot(n,ds,1,w,1)+ deltad*ddot(n,
     *dw,1,s,1) + deltap*deltad*ddot(n,ds,1,dw,1)
      mu = mu * ((g/mu)**3) /dfloat(2*n)
      do23026 i=1,n
      dr(i)=d(i)*(mu*(1/s(i)-1/x(i))+ dx(i)*dz(i)/x(i)-ds(i)*dw(i)/s(i))
23026 continue
23027 continue
      call dswap(p,rhs,1,dy,1)
      call dgemv('N',p,n,one,a,p,dr,1,one,dy,1)
      call dpotrs('U',p,1,ada,p,dy,p,info)
      call dgemv('T',p,n,one,a,p,dy,1,zero,u,1)
      deltap=big
      deltad=big
      do23028 i=1,n
      dxdz = dx(i)*dz(i)
      dsdw = ds(i)*dw(i)
      dx(i)= d(i)*(u(i)-z(i)+w(i))-dr(i)
      ds(i)= -dx(i)
      dz(i)= -z(i)+(mu - z(i)*dx(i) - dxdz)/x(i)
      dw(i)= -w(i)+(mu - w(i)*ds(i) - dsdw)/s(i)
      if(dx(i).lt.0)then
      deltap=dmin1(deltap,-x(i)/dx(i))
      endif
      if(ds(i).lt.0)then
      deltap=dmin1(deltap,-s(i)/ds(i))
      endif
      if(dz(i).lt.0)then
      deltad=dmin1(deltad,-z(i)/dz(i))
      endif
      if(dw(i).lt.0)then
      deltad=dmin1(deltad,-w(i)/dw(i))
      endif
23028 continue
23029 continue
      deltap=dmin1(beta*deltap,one)
      deltad=dmin1(beta*deltad,one)
      endif
      call daxpy(n,deltap,dx,1,x,1)
      call daxpy(n,deltap,ds,1,s,1)
      call daxpy(p,deltad,dy,1,y,1)
      call daxpy(n,deltad,dz,1,z,1)
      call daxpy(n,deltad,dw,1,w,1)
      gap = ddot(n,z,1,x,1)+ddot(n,w,1,s,1)
      goto 23008
      endif
23009 continue
      call daxpy(n,mone,w,1,z,1)
      call dswap(n,z,1,x,1)
      return
      end
      subroutine stepy(n,p,a,d,b,ada,info)
      integer n,p,pp,i,info
      double precision a(p,n),b(p),d(n),ada(p,p),zero
      parameter( zero = 0.0d0)
      pp=p*p
      do23038 j=1,p
      do23040 k=1,p
      ada(j,k)=zero
23040 continue
23041 continue
23038 continue
23039 continue
      do23042 i=1,n
      call dsyr('U',p,d(i),a(1,i),1,ada,p)
23042 continue
23043 continue
      call dposv('U',p,1,ada,p,b,p,info)
      return
      end


text:
*
*     -------------------------------------------------------------
*        B I S E C T
*     -------------------------------------------------------------
*
*
      SUBROUTINE BISECT(N,EPS1,D,E,E2,LB,UB,MM,M,W,IND,IERR,RV4,RV5)
*
      INTEGER I,J,K,L,M,N,P,Q,R,S,II,MM,M1,M2,TAG,IERR,ISTURM
      DOUBLE PRECISION D(N),E(N),E2(N),W(MM),RV4(N),RV5(N)
      DOUBLE PRECISION U,V,LB,T1,T2,UB,XU,X0,X1,EPS1,MACHEP
      DOUBLE PRECISION DABS,DMAX1,DMIN1,DFLOAT
      INTEGER IND(MM)
*
*     THIS SUBROUTINE IS A TRANSLATION OF THE BISECTION TECHNIQUE
*     IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON.
*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
*
*     THIS SUBROUTINE FINDS THOSE EIGENVALUES OF A TRIDIAGONAL
*     SYMMETRIC MATRIX WHICH LIE IN A SPECIFIED INTERVAL,
*     USING BISECTION.
*
*     ON INPUT:
*
*        N IS THE ORDER OF THE MATRIX;
*
*        EPS1 IS AN ABSOLUTE ERROR TOLERANCE FOR THE COMPUTED
*          EIGENVALUES.  IF THE INPUT EPS1 IS NON-POSITIVE,
*          IT IS RESET FOR EACH SUBMATRIX TO A DEFAULT VALUE,
*          NAMELY, MINUS THE PRODUCT OF THE RELATIVE MACHINE
*          PRECISION AND THE 1-NORM OF THE SUBMATRIX;
*
*        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX;
*
*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
*          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY;
*
*        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
*          E2(1) IS ARBITRARY;
*
*        LB AND UB DEFINE THE INTERVAL TO BE SEARCHED FOR EIGENVALUES.
*          IF LB IS NOT LESS THAN UB, NO EIGENVALUES WILL BE FOUND;
*
*        MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF
*          EIGENVALUES IN THE INTERVAL.  WARNING: IF MORE THAN
*          MM EIGENVALUES ARE DETERMINED TO LIE IN THE INTERVAL,
*          AN ERROR RETURN IS MADE WITH NO EIGENVALUES FOUND.
*
*     ON OUTPUT:
*
*        EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS
*          (LAST) DEFAULT VALUE;
*
*        D AND E ARE UNALTERED;
*
*        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED
*          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE
*          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES.
*          E2(1) IS ALSO SET TO ZERO;
*
*        M IS THE NUMBER OF EIGENVALUES DETERMINED TO LIE IN (LB,UB);
*
*        W CONTAINS THE M EIGENVALUES IN ASCENDING ORDER;
*
*        IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES
*          ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
*          1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM
*          THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.;
*
*        IERR IS SET TO
*          ZERO       FOR NORMAL RETURN,
*          3*N+1      IF M EXCEEDS MM;
*
*        RV4 AND RV5 ARE TEMPORARY STORAGE ARRAYS.
*
*     THE ALGOL PROCEDURE STURMCNT CONTAINED IN TRISTURM
*     APPEARS IN BISECT IN-LINE.
*
*     NOTE THAT SUBROUTINE TQL1 OR IMTQL1 IS GENERALLY FASTER THAN
*     BISECT, IF MORE THAN N/4 EIGENVALUES ARE TO BE FOUND.
*
*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
*     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
*
*     ------------------------------------------------------------------
*
*     :::::::::: MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
*                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
*                MACHEP = 16.0D0**(-13) FOR LONG FORM ARITHMETIC
*                ON S360 ::::::::::
      DATA MACHEP/1.D-12/
*
      IERR = 0
      TAG = 0
      T1 = LB
      T2 = UB
*     :::::::::: LOOK FOR SMALL SUB-DIAGONAL ENTRIES ::::::::::
      DO 40 I = 1, N
         IF (I .EQ. 1) GO TO 20
         IF (DABS(E(I)) .GT. MACHEP * (DABS(D(I)) + DABS(D(I-1))))
     :      GO TO 40
   20    E2(I) = 0.0D0
   40 CONTINUE
*     :::::::::: DETERMINE THE NUMBER OF EIGENVALUES
*                IN THE INTERVAL ::::::::::
      P = 1
      Q = N
      X1 = UB
      ISTURM = 1
      GO TO 320
   60 M = S
      X1 = LB
      ISTURM = 2
      GO TO 320
   80 M = M - S
      IF (M .GT. MM) GO TO 980
      Q = 0
      R = 0
*     :::::::::: ESTABLISH AND PROCESS NEXT SUBMATRIX, REFINING
*                INTERVAL BY THE GERSCHGORIN BOUNDS ::::::::::
  100 IF (R .EQ. M) GO TO 1001
      TAG = TAG + 1
      P = Q + 1
      XU = D(P)
      X0 = D(P)
      U = 0.0D0
*
      DO 120 Q = P, N
         X1 = U
         U = 0.0D0
         V = 0.0D0
         IF (Q .EQ. N) GO TO 110
         U = DABS(E(Q+1))
         V = E2(Q+1)
  110    XU = DMIN1(D(Q)-(X1+U),XU)
         X0 = DMAX1(D(Q)+(X1+U),X0)
         IF (V .EQ. 0.0D0) GO TO 140
  120 CONTINUE
*
  140 X1 = DMAX1(DABS(XU),DABS(X0)) * MACHEP
      IF (EPS1 .LE. 0.0D0) EPS1 = -X1
      IF (P .NE. Q) GO TO 180
*     :::::::::: CHECK FOR ISOLATED ROOT WITHIN INTERVAL ::::::::::
      IF (T1 .GT. D(P) .OR. D(P) .GE. T2) GO TO 940
      M1 = P
      M2 = P
      RV5(P) = D(P)
      GO TO 900
  180 X1 = X1 * DFLOAT(Q-P+1)
      LB = DMAX1(T1,XU-X1)
      UB = DMIN1(T2,X0+X1)
      X1 = LB
      ISTURM = 3
      GO TO 320
  200 M1 = S + 1
      X1 = UB
      ISTURM = 4
      GO TO 320
  220 M2 = S
      IF (M1 .GT. M2) GO TO 940
*     :::::::::: FIND ROOTS BY BISECTION ::::::::::
      X0 = UB
      ISTURM = 5
*
      DO 240 I = M1, M2
         RV5(I) = UB
         RV4(I) = LB
  240 CONTINUE
*     :::::::::: LOOP FOR K-TH EIGENVALUE
*                FOR K=M2 STEP -1 UNTIL M1 DO --
*                (-DO- NOT USED TO LEGALIZE COMPUTED-GO-TO) ::::::::::
      K = M2
  250    XU = LB
*     :::::::::: FOR I=K STEP -1 UNTIL M1 DO -- ::::::::::
         DO 260 II = M1, K
            I = M1 + K - II
            IF (XU .GE. RV4(I)) GO TO 260
            XU = RV4(I)
            GO TO 280
  260    CONTINUE
*
  280    IF (X0 .GT. RV5(K)) X0 = RV5(K)
*     :::::::::: NEXT BISECTION STEP ::::::::::
  300    X1 = (XU + X0) * 0.5D0
         IF ((X0 - XU) .LE. (2.0D0 * MACHEP *
     :      (DABS(XU) + DABS(X0)) + DABS(EPS1))) GO TO 420
*     :::::::::: IN-LINE PROCEDURE FOR STURM SEQUENCE ::::::::::
  320    S = P - 1
         U = 1.0D0
*
         DO 340 I = P, Q
            IF (U .NE. 0.0D0) GO TO 325
            V = DABS(E(I)) / MACHEP
            GO TO 330
  325       V = E2(I) / U
  330       U = D(I) - X1 - V
            IF (U .LT. 0.0D0) S = S + 1
  340    CONTINUE
*
         GO TO (60,80,200,220,360), ISTURM
*     :::::::::: REFINE INTERVALS ::::::::::
  360    IF (S .GE. K) GO TO 400
         XU = X1
         IF (S .GE. M1) GO TO 380
         RV4(M1) = X1
         GO TO 300
  380    RV4(S+1) = X1
         IF (RV5(S) .GT. X1) RV5(S) = X1
         GO TO 300
  400    X0 = X1
         GO TO 300
*     :::::::::: K-TH EIGENVALUE FOUND ::::::::::
  420    RV5(K) = X1
      K = K - 1
      IF (K .GE. M1) GO TO 250
*     :::::::::: ORDER EIGENVALUES TAGGED WITH THEIR
*                SUBMATRIX ASSOCIATIONS ::::::::::
  900 S = R
      R = R + M2 - M1 + 1
      J = 1
      K = M1
*
      DO 920 L = 1, R
         IF (J .GT. S) GO TO 910
         IF (K .GT. M2) GO TO 940
         IF (RV5(K) .GE. W(L)) GO TO 915
*
         DO 905 II = J, S
            I = L + S - II
            W(I+1) = W(I)
            IND(I+1) = IND(I)
  905    CONTINUE
*
  910    W(L) = RV5(K)
         IND(L) = TAG
         K = K + 1
         GO TO 920
  915    J = J + 1
  920 CONTINUE
*
  940 IF (Q .LT. N) GO TO 100
      GO TO 1001
*     :::::::::: SET ERROR -- UNDERESTIMATE OF NUMBER OF
*                EIGENVALUES IN INTERVAL ::::::::::
  980 IERR = 3 * N + 1
 1001 LB = T1
      UB = T2
      RETURN
*     :::::::::: LAST CARD OF BISECT ::::::::::
      END


text:
ccalculate loan payment
c     another way of comment
*  yet another way of comment
      PROGRAM LOAN
      WRITE(*,*) 'amount, percentage rate, and years of loan'
      READ(*,*) AMOUNT, PCRATE, NYEARS
      RATE = PCRATE /100.0
      REPAY = RATE * AMOUNT /(1.0 - (1.0+RATE)**(-NYEARS))
      MTPAY = REPAY/12.0
      WRITE(*,*) 'Annual repayments are ', REPAY
      WRITE(*,*) 'Monthly payment is ', MTPAY
      WRITE(*,*) 'End of year  Balance'
      DO 15, IYEAR = 1,NYEARS
        AMOUNT = AMOUNT +(AMOUNT * RATE) - REPAY
        WRITE(*,'(1X,I9,F11.2)') IYEAR, AMOUNT
c        WRITE(*,*) IYEAR, AMOUNT
15    CONTINUE
      END
c test comments after END

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine extrapolate_gradtel(nface,ielfa,xrlfa,gradtel,gradtfa)
!
!     inter/extrapolation of p element values to the faces
!
      implicit none
!
      integer nface,ielfa(4,*),i,iel1,iel2,l
!
      real*8 xrlfa(3,*),gradtel(3,*),gradtfa(3,*),xl1
!
      do i=1,nface
         iel1=ielfa(1,i)
         xl1=xrlfa(1,i)
         iel2=ielfa(2,i)
         if(iel2.gt.0) then
!
!           face in between two elements
!
            do l=1,3
               gradtfa(l,i)=xl1*gradtel(l,iel1)+
     &              xrlfa(2,i)*gradtel(l,iel2)
            enddo
         elseif(ielfa(3,i).gt.0) then
!
!     boundary face; linear extrapolation
!
            do l=1,3
               gradtfa(l,i)=xl1*gradtel(l,iel1)+
     &              xrlfa(3,i)*gradtel(l,ielfa(3,i))
            enddo
         else
!
!     boundary face; constant extrapolation (one element layer)
!
            do l=1,3
               gradtfa(l,i)=gradtel(l,iel1)
            enddo
         endif
c         write(*,*) 'extrapolate_gradtel ',(gradtfa(l,i),l=1,3)
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
subroutine init_genes(jean)
    use numz
    use ran_mod
    use control
    implicit none
    integer,dimension(:,:)::jean
    real(b8) y
    integer i,j,k(2)
    k=ubound(jean) ! return a 2d array with both upper bounds
    do i=1,k(1)
      do j=1,k(2)
        jean(i,j)=int(spread(0.0_b8,real(ncolor,b8))) ! gives us range of 0-ncolor-1
      enddo
    enddo
end subroutine

! our function that returns a vector
! take an integer input vector which represents
! an integer in some base and adds 1
function add1(vector,max) result (rtn)
  integer, dimension(:),intent(in) ::  vector
  integer,dimension(size(vector)) :: rtn
  integer max
  integer len
  logical carry
  len=size(vector)
  rtn=vector
  i=0
  carry=.true.
  do while(carry)
      i=i+1
      rtn(i)=rtn(i)+1
      if(rtn(i) .gt. max)then
          if(i == len)then
              rtn=0
          else
              rtn(i)=0
          endif
      else
          carry=.false.
      endif
  enddo
end function

max_stars_repo_name: timkphd/examples

max_issues_repo_name: timkphd/examples

max_forks_repo_name: timkphd/examples


text:
use s2d, render_board

game_loop(board):
        render_board( board, 150, 100, 50, 10 ).

s2d_render():
    board <- [1 | 0 <= c < 2000 ]
    game_loop(board).

main():
    s2d_window('board test',  800, 600 ).

max_stars_repo_name: diegovalverde/funk

max_issues_repo_name: diegovalverde/funk

max_forks_repo_name: diegovalverde/funk