!****h* root/fortran/test/tH5D.f90
!
! NAME
!  tH5D.f90
!
! FUNCTION
!  Basic testing of Fortran H5D APIs.
!
! COPYRIGHT
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   Copyright by The HDF Group.                                               *
!   Copyright by the Board of Trustees of the University of Illinois.         *
!   All rights reserved.                                                      *
!                                                                             *
!   This file is part of HDF5.  The full HDF5 copyright notice, including     *
!   terms governing use, modification, and redistribution, is contained in    *
!   the COPYING file, which can be found at the root of the source code       *
!   distribution tree, or in https://support.hdfgroup.org/ftp/HDF5/releases.  *
!   If you do not have access to either file, you may request a copy from     *
!   help@hdfgroup.org.                                                        *
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
! NOTES
!  Tests the H5D APIs functionalities of:
!   h5dcreate_f, h5dopen_f, h5dclose_f, h5dget_space_f, h5dget_type_f,
!   h5dread_f, and h5dwrite_f, h5dget_space_status_f
!
!
! CONTAINS SUBROUTINES
!  datasettest, extenddsettest
!
!*****

!
MODULE TH5D

  USE HDF5 ! This module contains all necessary modules
  USE TH5_MISC
  USE TH5_MISC_GEN

CONTAINS
  SUBROUTINE datasettest(cleanup, total_error)

    IMPLICIT NONE
    LOGICAL, INTENT(IN) :: cleanup
    INTEGER, INTENT(OUT) :: total_error

    CHARACTER(LEN=5), PARAMETER :: filename = "dsetf" ! File name
    CHARACTER(LEN=80) :: fix_filename
    CHARACTER(LEN=4), PARAMETER :: dsetname = "dset"     ! Dataset name
    CHARACTER(LEN=9), PARAMETER :: null_dsetname = "null_dset"     ! Dataset name

    INTEGER(HID_T) :: file_id       ! File identifier
    INTEGER(HID_T) :: dset_id       ! Dataset identifier
    INTEGER(HID_T) :: null_dset     ! Null dataset identifier
    INTEGER(HID_T) :: dspace_id     ! Dataspace identifier
    INTEGER(HID_T) :: null_dspace   ! Null dataspace identifier
    INTEGER(HID_T) :: dtype_id      ! Datatype identifier

    INTEGER(HSIZE_T), DIMENSION(2) :: dims = (/4,6/) ! Dataset dimensions
    INTEGER     ::   rank = 2                        ! Dataset rank

    INTEGER, DIMENSION(4,6) :: dset_data, data_out ! Data buffers
    INTEGER     ::   error ! Error flag

    INTEGER     :: i, j    !general purpose integers
    INTEGER(HSIZE_T), DIMENSION(2) :: data_dims
    INTEGER(HSIZE_T), DIMENSION(1) :: null_data_dim
    INTEGER     ::   null_dset_data = 1              ! null data
    INTEGER :: flag ! Space allocation status

    !
    ! Initialize the dset_data array.
    !
    DO i = 1, 4
       DO j = 1, 6
          dset_data(i,j) = (i-1)*6 + j;
       END DO
    END DO
    !
    ! Create a new file using default properties.
    !
    CALL h5_fixname_f(filename, fix_filename, H5P_DEFAULT_F, error)
    IF (error .NE. 0) THEN
       WRITE(*,*) "Cannot modify filename"
       STOP
    ENDIF
    CALL h5fcreate_f(fix_filename, H5F_ACC_TRUNC_F, file_id, error)
    CALL check("h5fcreate_f", error, total_error)
    !
    ! Create the dataspace.
    !
    CALL h5screate_simple_f(rank, dims, dspace_id, error)
    CALL check("h5screate_simple_f", error, total_error)
    !
    ! Create null dataspace.
    !
    CALL h5screate_f(H5S_NULL_F, null_dspace, error)
    CALL check("h5screate_simple_f", error, total_error)
    !
    ! Create the dataset with default properties.
    !
    CALL h5dcreate_f(file_id, dsetname, H5T_NATIVE_INTEGER, dspace_id, &
         dset_id, error)
    CALL check("h5dcreate_f", error, total_error)
    !
    ! Create the null dataset.
    !
    CALL h5dcreate_f(file_id, null_dsetname, H5T_NATIVE_INTEGER, null_dspace, null_dset, error)
    CALL check("h5dcreate_f", error, total_error)
    !
    ! Write the dataset.
    !
    data_dims(1) = 4
    data_dims(2) = 6
    CALL h5dwrite_f(dset_id, H5T_NATIVE_INTEGER, dset_data, data_dims, error)
    CALL check("h5dwrite_f", error, total_error)
    !
    ! Write null dataset.  Nothing can be written.
    !
    null_data_dim(1) = 1
    CALL h5dwrite_f(null_dset, H5T_NATIVE_INTEGER, null_dset_data, null_data_dim, error)
    CALL check("h5dwrite_f", error, total_error)
    !
    ! End access to the dataset and release resources used by it.
    !
    CALL h5dclose_f(dset_id, error)
    CALL check("h5dclose_f", error, total_error)
    CALL h5dclose_f(null_dset, error)
    CALL check("h5dclose_f", error, total_error)
    !
    ! Terminate access to the data space.
    !
    CALL h5sclose_f(dspace_id, error)
    CALL check("h5sclose_f", error, total_error)
    CALL h5sclose_f(null_dspace, error)
    CALL check("h5sclose_f", error, total_error)
    !
    ! Close the file.
    !
    CALL h5fclose_f(file_id, error)
    CALL check("h5fclose_f", error, total_error)
    !
    ! Open the existing file.
    !
    CALL h5fopen_f (fix_filename, H5F_ACC_RDWR_F, file_id, error)
    CALL check("h5fopen_f", error, total_error)
    !
    ! Open the existing dataset.
    !
    CALL h5dopen_f(file_id, dsetname, dset_id, error)
    CALL check("h5dopen_f", error, total_error)
    CALL h5dopen_f(file_id, null_dsetname, null_dset, error)
    CALL check("h5dopen_f", error, total_error)

    ! Test whether space has been allocated for a dataset
    CALL h5dget_space_status_f(dset_id, flag, error)
    CALL check("h5dget_space_status_f",error, total_error)
    CALL VERIFY("h5dget_space_status_f", flag, H5D_SPACE_STS_ALLOCATED_F, total_error)

    CALL h5dget_space_status_f(null_dset, flag, error)
    CALL check("h5dget_space_status_f",error, total_error)
    CALL VERIFY("h5dget_space_status_f", flag, H5D_SPACE_STS_NOT_ALLOCATED_F, total_error)
    !
    ! Get the dataset type.
    !
    CALL h5dget_type_f(dset_id, dtype_id, error)
    CALL check("h5dget_type_f", error, total_error)
    !
    ! Get the data space.
    !
    CALL h5dget_space_f(dset_id, dspace_id, error)
    CALL check("h5dget_space_f", error, total_error)
    !
    ! Read the dataset.
    !
    CALL h5dread_f(dset_id, H5T_NATIVE_INTEGER, data_out, data_dims, error)
    CALL check("h5dread_f", error, total_error)
    !
    ! Read the null dataset.  Nothing should be read.
    !
    CALL h5dread_f(null_dset, H5T_NATIVE_INTEGER, null_dset_data, null_data_dim, error)
    CALL check("h5dread_f", error, total_error)
    !
    !Compare the data.
    !
    DO i = 1, 4
       DO j = 1, 6
          IF (data_out(i,j) .NE. dset_data(i, j)) THEN
             WRITE(*, *) "dataset test error occured"
             WRITE(*,*) "data read is not the same as the data written"
          END IF
       END DO
    END DO
    !
    ! Check if no change to null_dset_data
    !
    IF (null_dset_data .NE. 1) THEN
       WRITE(*, *) "null dataset test error occured"
    END IF
    !
    ! End access to the dataset and release resources used by it.
    !
    CALL h5dclose_f(dset_id, error)
    CALL check("h5dclose_f", error, total_error)
    CALL h5dclose_f(null_dset, error)
    CALL check("h5dclose_f", error, total_error)
    !
    ! Terminate access to the data space.
    !
    CALL h5sclose_f(dspace_id, error)
    CALL check("h5sclose_f", error, total_error)

    !
    ! Terminate access to the data type.
    !
    CALL h5tclose_f(dtype_id, error)
    CALL check("h5tclose_f", error, total_error)
    !
    ! Close the file.
    !
    CALL h5fclose_f(file_id, error)
    CALL check("h5fclose_f", error, total_error)
    IF(cleanup) CALL h5_cleanup_f(filename, H5P_DEFAULT_F, error)
    CALL check("h5_cleanup_f", error, total_error)

    RETURN
  END SUBROUTINE datasettest

!
!the following subroutine tests h5dextend_f functionality
!

  SUBROUTINE extenddsettest(cleanup, total_error)

    IMPLICIT NONE

    LOGICAL, INTENT(IN)  :: cleanup
    INTEGER, INTENT(OUT) :: total_error

    !
    !the dataset is stored in file "extf.h5"
    !
    CHARACTER(LEN=4), PARAMETER :: filename = "extf"
    CHARACTER(LEN=80) :: fix_filename

    !
    !dataset name is "ExtendibleArray"
    !
    CHARACTER(LEN=15), PARAMETER :: dsetname = "ExtendibleArray"

    !
    !dataset rank is 2
    !
    INTEGER :: RANK = 2

    INTEGER(HID_T) :: file_id       ! File identifier
    INTEGER(HID_T) :: dset_id       ! Dataset identifier
    INTEGER(HID_T) :: dataspace     ! Dataspace identifier
    INTEGER(HID_T) :: memspace      ! memory Dataspace identifier
    INTEGER(HID_T) :: crp_list        ! dataset creatation property identifier

    !
    !dataset dimensions at creation time
    !
    INTEGER(HSIZE_T), DIMENSION(2) :: dims = (/3,3/)

    !
    !data dimensions
    !
    INTEGER(HSIZE_T), DIMENSION(2) :: dims1 = (/10,3/)

    !
    !Maximum dimensions
    !
    INTEGER(HSIZE_T), DIMENSION(2) :: maxdims

    !
    !data arrays for reading and writing
    !
    INTEGER, DIMENSION(10,3) :: data_in, data_out

    !
    !Size of data in the file
    !
    INTEGER(HSIZE_T), DIMENSION(2) :: size

    !
    !general purpose integer
    !
    INTEGER :: i, j
    INTEGER(HSIZE_T) :: ih, jh

    !
    !flag to check operation success
    !
    INTEGER :: error

    !
    !Variables used in reading data back
    !
    INTEGER(HSIZE_T), DIMENSION(2) :: dimsr, maxdimsr
    INTEGER :: rankr
    INTEGER(HSIZE_T), DIMENSION(2) :: data_dims

    !
    !data initialization
    !
    DO i = 1, 10
       DO j = 1, 3
          data_in(i,j) = 2
       END DO
    END DO

    !
    !Initialize FORTRAN predifined datatypes
    !
!          CALL h5init_types_f(error)
!               CALL check("h5init_types_f",error,total_error)

    !
    !Create a new file using default properties.
    !
    CALL h5_fixname_f(filename, fix_filename, H5P_DEFAULT_F, error)
    IF (error .NE. 0) THEN
       WRITE(*,*) "Cannot modify filename"
       STOP
    ENDIF
    CALL h5fcreate_f(fix_filename, H5F_ACC_TRUNC_F, file_id, error)
    CALL check("h5fcreate_f",error,total_error)

    !
    !Create the data space with unlimited dimensions.
    !
    maxdims = (/H5S_UNLIMITED_F, H5S_UNLIMITED_F/)

    CALL h5screate_simple_f(RANK, dims, dataspace, error, maxdims)
    CALL check("h5screate_simple_f",error,total_error)

    !
    !Modify dataset creation properties, i.e. enable chunking
    !
    CALL h5pcreate_f(H5P_DATASET_CREATE_F, crp_list, error)
    CALL check("h5pcreate_f",error,total_error)

    CALL h5pset_chunk_f(crp_list, RANK, dims1, error)
    CALL check("h5pset_chunk_f",error,total_error)

    !
    !Create a dataset with 3X3 dimensions using cparms creation propertie .
    !
    CALL h5dcreate_f(file_id, dsetname, H5T_NATIVE_INTEGER, dataspace, dset_id, error, crp_list )
    CALL check("h5dcreate_f",error,total_error)

    !
    !Extend the dataset. This call assures that dataset is 3 x 3.
    !
    SIZE(1) = 3
    SIZE(2) = 3
    CALL h5dextend_f(dset_id, size, error)
    CALL check("h5dextend_f",error,total_error)

    !
    !Extend the dataset. Dataset becomes 10 x 3.
    !
    SIZE(1)   = 10;
    SIZE(2)   = 3;
    CALL h5dextend_f(dset_id, size, error)
    CALL check("h5dextend_f",error,total_error)

    !
    !Write the data of size 10X3 to the extended dataset.
    !
    data_dims(1) = 10
    data_dims(2) = 3
    CALL H5dwrite_f(dset_id, H5T_NATIVE_INTEGER, data_in, data_dims, error)
    CALL check("h5dwrite_f",error,total_error)

    !
    !Close the dataspace for the dataset.
    !
    CALL h5sclose_f(dataspace, error)
    CALL check("h5sclose_f",error,total_error)

    !
    !Close the property list.
    !
    CALL h5pclose_f(crp_list, error)
    CALL check("h5pclose_f",error,total_error)
    !
    !Close the dataset.
    !
    CALL h5dclose_f(dset_id, error)
    CALL check("h5dclose_f",error,total_error)

    !
    !Close the file.
    !
    CALL h5fclose_f(file_id, error)
    CALL check("h5fclose_f",error,total_error)

    !
    !read the data back
    !
    !Open the file.
    !
    CALL h5fopen_f (fix_filename, H5F_ACC_RDONLY_F, file_id, error)
    CALL check("hfopen_f",error,total_error)

    !
    !Open the  dataset.
    !
    CALL h5dopen_f(file_id, dsetname, dset_id, error)
    CALL check("h5dopen_f",error,total_error)

    !
    !Get dataset's dataspace handle.
    !
    CALL h5dget_space_f(dset_id, dataspace, error)
    CALL check("h5dget_space_f",error,total_error)

    !
    !Get dataspace's rank.
    !
    CALL h5sget_simple_extent_ndims_f(dataspace, rankr, error)
    CALL check("h5sget_simple_extent_ndims_f",error,total_error)
    IF (rankr .NE. RANK) THEN
       WRITE(*,*) "dataset rank error occured"
       STOP
    END IF

    !
    !Get dataspace's dimensinons.
    !
    CALL h5sget_simple_extent_dims_f(dataspace, dimsr, maxdimsr, error)
    CALL check("h5sget_simple_extent_dims_f",error,total_error)
    IF ((dimsr(1) .NE. dims1(1)) .OR. (dimsr(2) .NE. dims1(2))) THEN
       WRITE(*,*) "dataset dimensions error occured"
       STOP
    END IF

    !
    !Get creation property list.
    !
    CALL h5dget_create_plist_f(dset_id, crp_list, error)
    CALL check("h5dget_create_plist_f",error,total_error)

    !
    !create memory dataspace.
    !
    CALL h5screate_simple_f(rankr, dimsr, memspace, error)
    CALL check("h5screate_simple_f",error,total_error)

    !
    !Read data
    !
    CALL H5dread_f(dset_id, H5T_NATIVE_INTEGER, data_out, data_dims, error, memspace, dataspace)
    CALL check("h5dread_f",error,total_error)

    !
    !Compare the data.
    !
    DO ih = 1, dims1(1)
       DO jh = 1, dims1(2)
          IF (data_out(ih,jh) .NE. data_in(ih, jh)) THEN
             WRITE(*, *) "extend dataset test error occured"
             WRITE(*, *) "read value is not the same as the written values"
          END IF
       END DO
    END DO

    !
    !Close the dataspace for the dataset.
    !
    CALL h5sclose_f(dataspace, error)
    CALL check("h5sclose_f",error,total_error)

    !
    !Close the memspace for the dataset.
    !
    CALL h5sclose_f(memspace, error)
    CALL check("h5sclose_f",error,total_error)

    !
    !Close the property list.
    !
    CALL h5pclose_f(crp_list, error)
    CALL check("h5pclose_f",error,total_error)

    !
    !Close the dataset.
    !
    CALL h5dclose_f(dset_id, error)
    CALL check("h5dclose_f",error,total_error)

    !
    !Close the file.
    !
    CALL h5fclose_f(file_id, error)
    CALL check("h5fclose_f",error,total_error)
    IF(cleanup) CALL h5_cleanup_f(filename, H5P_DEFAULT_F, error)
    CALL check("h5_cleanup_f", error, total_error)

    RETURN
  END SUBROUTINE extenddsettest

!
! The following subroutine tests h5dget_offset_f functionality
!

  SUBROUTINE test_userblock_offset(cleanup, total_error)

    USE ISO_C_BINDING

    IMPLICIT NONE
    LOGICAL, INTENT(IN) :: cleanup
    INTEGER, INTENT(OUT) :: total_error
    !
    !the dataset is stored in file "offset.h5"
    !
    INTEGER, PARAMETER :: dset_dim1=2, dset_dim2=10
    CHARACTER(LEN=6), PARAMETER :: filename = "offset"
    CHARACTER(LEN=80) :: fix_filename

    INTEGER(hid_t) :: file, fcpl, dataset, space
    INTEGER :: i, j, n, ios
    INTEGER(hsize_t), DIMENSION(1:2) :: dims
    INTEGER(haddr_t) :: offset
    INTEGER, DIMENSION(1:dset_dim1,1:dset_dim2), TARGET :: rdata, data_in
    INTEGER :: error
    TYPE(C_PTR) :: f_ptr
    !
    !Create a new file using default properties.
    !
    CALL h5_fixname_f(filename, fix_filename, H5P_DEFAULT_F, error)
    IF (error .NE. 0) THEN
       WRITE(*,*) "Cannot modify filename"
       STOP
    ENDIF

    CALL h5pcreate_f(H5P_FILE_CREATE_F, fcpl, error)
    CALL check("h5pcreate_f",error,total_error)

    ! Initialize the dataset
    n = 0
    DO i = 1, dset_dim1
       DO j = 1, dset_dim2
          n = n + 1
          data_in(i,j) = n
       END DO
    END DO
    CALL h5fcreate_f(fix_filename, H5F_ACC_TRUNC_F, file, error, fcpl)
    CALL check("h5fcreate_f",error,total_error)

    ! Create the data space
    dims(1:2) = (/dset_dim1,dset_dim2/)

    CALL h5screate_simple_f(2, dims, space, error)
    CALL check("h5screate_simple_f",error,total_error)

    ! Create the dataset
    CALL h5dcreate_f(file, "dset1", H5T_NATIVE_INTEGER, space, dataset, error)
    CALL check("h5dcreate_f", error, total_error)

    ! Test dataset address.  Should be undefined.
    CALL h5dget_offset_f(dataset, offset, error)
    CALL VERIFY("h5dget_offset_f",offset, HADDR_UNDEF_F, total_error)

    ! Write the data to the dataset
    f_ptr = C_LOC(data_in(1,1))
    CALL h5dwrite_f(dataset, H5T_NATIVE_INTEGER, f_ptr, error)
    CALL check("h5dwrite_f", error, total_error)

    ! Test dataset address in file. Open the same file as a C file, seek
    ! the data position as H5Dget_offset points to, read the dataset, and
    ! compare it with the data written in.
    CALL h5dget_offset_f(dataset, offset, error)
    CALL check("h5dget_offset_f", error, total_error)
    IF(offset.EQ.HADDR_UNDEF_F)THEN
       total_error = total_error + 1
    ENDIF

    CALL h5dclose_f(dataset, error)
    CALL check("h5dclose_f", error, total_error)
    CALL h5fclose_f(file, error)
    CALL check("h5fclose_f", error, total_error)

    IF(total_error.NE.0) RETURN

    OPEN(10,FILE=fix_filename, ACCESS="STREAM", IOSTAT=ios)
    IF(ios.NE.0)THEN
       WRITE(*,'(A)') "Failed to open file "//TRIM(fix_filename)
       total_error = total_error + 1
       RETURN
    ENDIF
    ! The pos= specifier illustrates that positions are in bytes,
    ! starting from byte 1 (as opposed to C, where they start from byte 0)
    READ(10, POS=offset+1, IOSTAT=ios) rdata
    IF(ios.NE.0)THEN
       WRITE(*,'(A)') "Failed to read data from stream I/O "
       total_error = total_error + 1
       CLOSE(10)
       RETURN
    ENDIF

    ! Check that the values read are the same as the values written
    DO i = 1, dset_dim1
       DO j = 1, dset_dim2
          CALL VERIFY("h5dget_offset_f",rdata(i,j), data_in(i,j), total_error)
          IF(total_error.NE.0)THEN
             WRITE(*,'(A)') "    Read different values than written."
             WRITE(*,'(2(A,I0))') "    At index ",i,",",j
             CLOSE(10)
             RETURN
          ENDIF
       END DO
    END DO

    CLOSE(10)

    IF(cleanup) CALL h5_cleanup_f(fix_filename, H5P_DEFAULT_F, error)
    CALL check("h5_cleanup_f", error, total_error)
    IF(cleanup) CALL h5_cleanup_f(filename, H5P_DEFAULT_F, error)
    CALL check("h5_cleanup_f", error, total_error)

  END SUBROUTINE test_userblock_offset

END MODULE TH5D


      subroutine nut_pminrl

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine computes p flux between the labile, active mineral
!!    and stable mineral p pools.
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    Intrinsic: Min

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use basin_module
      use organic_mineral_mass_module
      use hru_module, only : ihru
      use soil_module
      use output_landscape_module, only : hnb_d

      implicit none

      real, parameter :: bk = .0006     !              |
      integer :: j                      !none          |HRU number
      integer :: l                      !none          |counter
      real :: rto                       !              |
      real :: rmp1                      !kg P/ha       |amount of phosphorus moving from the solution
                                        !              |mineral to the active mineral pool in the
                                        !              |soil layer
      real :: roc                       !kg P/ha       |amount of phosphorus moving from the active
                                        !              |mineral to the stable mineral pool in the
                                        !              |soil layer

      j = ihru

      hnb_d(j)%lab_min_p = 0.
      hnb_d(j)%act_sta_p = 0.

      rto = bsn_prm%psp / (1. - bsn_prm%psp)

      do l = 1, soil(j)%nly
        rmp1 = (soil1(j)%mp(l)%lab - soil1(j)%mp(l)%act * rto)
        !! mike changed/added per isabelle beaudin"s email from 01/21/09
        if (rmp1 > 0.) rmp1 = rmp1 * 0.1
        if (rmp1 < 0.) rmp1 = rmp1 * 0.6
        !! mike changed/added per isabelle beaudin"s email from 01/21/09
        rmp1 = Min(rmp1, soil1(j)%mp(l)%lab)

        roc = bk * (4. * soil1(j)%mp(l)%act - soil1(j)%mp(l)%sta)
        if (roc < 0.) roc = roc * .1
        roc = Min(roc, soil1(j)%mp(l)%act)

        soil1(j)%mp(l)%sta = soil1(j)%mp(l)%sta + roc
        if (soil1(j)%mp(l)%sta < 0.) soil1(j)%mp(l)%sta = 0.

        soil1(j)%mp(l)%act = soil1(j)%mp(l)%act - roc + rmp1
        if (soil1(j)%mp(l)%act < 0.) soil1(j)%mp(l)%act = 0.

        soil1(j)%mp(l)%lab = soil1(j)%mp(l)%lab - rmp1
        if (soil1(j)%mp(l)%lab < 0.) soil1(j)%mp(l)%lab = 0.

        hnb_d(j)%lab_min_p = hnb_d(j)%lab_min_p + rmp1
        hnb_d(j)%act_sta_p = hnb_d(j)%act_sta_p + roc
      end do

      return
      end subroutine nut_pminrl


text:
! RUN: bbc -emit-fir %s -o - | FileCheck %s
! RUN: %flang_fc1 -emit-fir %s -o - | FileCheck %s

! CHECK-LABEL: exp_testr
! CHECK-SAME: (%[[AREF:.*]]: !fir.ref<f32> {{.*}}, %[[BREF:.*]]: !fir.ref<f32> {{.*}})
subroutine exp_testr(a, b)
  real :: a, b
! CHECK:  %[[A:.*]] = fir.load %[[AREF:.*]] : !fir.ref<f32>
! CHECK:  %[[RES:.*]] = fir.call @fir.exp.f32.f32(%[[A]]) : (f32) -> f32
! CHECK:  fir.store %[[RES]] to %[[BREF]] : !fir.ref<f32>
  b = exp(a)
end subroutine

! CHECK-LABEL: exp_testd
! CHECK-SAME: (%[[AREF:.*]]: !fir.ref<f64> {{.*}}, %[[BREF:.*]]: !fir.ref<f64> {{.*}})
subroutine exp_testd(a, b)
  real(kind=8) :: a, b
! CHECK:  %[[A:.*]] = fir.load %[[AREF:.*]] : !fir.ref<f64>
! CHECK:  %[[RES:.*]] = fir.call @fir.exp.f64.f64(%[[A]]) : (f64) -> f64
! CHECK:  fir.store %[[RES]] to %[[BREF]] : !fir.ref<f64>
  b = exp(a)
end subroutine

! CHECK-LABEL: exp_testc
! CHECK-SAME: (%[[AREF:.*]]: !fir.ref<!fir.complex<4>> {{.*}}, %[[BREF:.*]]: !fir.ref<!fir.complex<4>> {{.*}})
subroutine exp_testc(a, b)
  complex :: a, b
! CHECK:  %[[A:.*]] = fir.load %[[AREF:.*]] : !fir.ref<!fir.complex<4>>
! CHECK:  %[[RES:.*]] = fir.call @fir.exp.z4.z4(%[[A]]) : (!fir.complex<4>) -> !fir.complex<4>
! CHECK:  fir.store %[[RES]] to %[[BREF]] : !fir.ref<!fir.complex<4>>
  b = exp(a)
end subroutine

! CHECK-LABEL: exp_testcd
! CHECK-SAME: (%[[AREF:.*]]: !fir.ref<!fir.complex<8>> {{.*}}, %[[BREF:.*]]: !fir.ref<!fir.complex<8>> {{.*}})
subroutine exp_testcd(a, b)
  complex(kind=8) :: a, b
! CHECK:  %[[A:.*]] = fir.load %[[AREF:.*]] : !fir.ref<!fir.complex<8>>
! CHECK:  %[[RES:.*]] = fir.call @fir.exp.z8.z8(%[[A]]) : (!fir.complex<8>) -> !fir.complex<8>
! CHECK:  fir.store %[[RES]] to %[[BREF]] : !fir.ref<!fir.complex<8>>
  b = exp(a)
end subroutine

! CHECK-LABEL: private @fir.exp.f32.f32
! CHECK-SAME: (%[[ARG32_OUTLINE:.*]]: f32) -> f32
! CHECK: %[[RESULT32_OUTLINE:.*]] = fir.call @__fs_exp_1(%[[ARG32_OUTLINE]]) : (f32) -> f32
! CHECK: return %[[RESULT32_OUTLINE]] : f32

! CHECK-LABEL: private @fir.exp.f64.f64
! CHECK-SAME: (%[[ARG64_OUTLINE:.*]]: f64) -> f64
! CHECK: %[[RESULT64_OUTLINE:.*]] = fir.call @__fd_exp_1(%[[ARG64_OUTLINE]]) : (f64) -> f64
! CHECK: return %[[RESULT64_OUTLINE]] : f64

! CHECK-LABEL: private @fir.exp.z4.z4
! CHECK-SAME: (%[[ARG32_OUTLINE]]: !fir.complex<4>) -> !fir.complex<4>
! CHECK: %[[RESULT32_OUTLINE]] = fir.call @__fc_exp_1(%[[ARG32_OUTLINE]]) : (!fir.complex<4>) -> !fir.complex<4>
! CHECK: return %[[RESULT32_OUTLINE]] : !fir.complex<4>

! CHECK-LABEL: private @fir.exp.z8.z8
! CHECK-SAME: (%[[ARG64_OUTLINE]]: !fir.complex<8>) -> !fir.complex<8>
! CHECK: %[[RESULT64_OUTLINE]] = fir.call @__fz_exp_1(%[[ARG64_OUTLINE]]) : (!fir.complex<8>) -> !fir.complex<8>
! CHECK: return %[[RESULT64_OUTLINE]] : !fir.complex<8>

max_stars_repo_name: akashhansda/llvm-project

max_issues_repo_name: akashhansda/llvm-project

max_forks_repo_name: akashhansda/llvm-project


text:
C->>> -----------------------------------------------> ems_g_emsol_p <<<
c     Get the EMSOL pointers
c
      subroutine ems_g_emsol_p(rt_cod, is)
      implicit none
      integer rt_cod, is(0:*)
c      print*, 'Calling TRUE ems_g_emsol_p'
      call ems_g_ml_p(rt_cod, is)
      if (rt_cod .ne. 0) goto 7000
      call ems_g_rsmi_p(rt_cod, is)
      if (rt_cod .ne. 0) goto 7000
      call ems_g_inv_p(rt_cod, is)
      if (rt_cod .ne. 0) goto 7000
 7000 continue
      return
      end

C->>> -----------------------------------------> ems_se_emsol_p_no_p <<<
c     Get the EMSOL pointers
c
      subroutine ems_se_emsol_p_no_p(rt_cod, is)
      implicit none
      integer rt_cod, is(0:*)
c      print*, 'Calling TRUE ems_se_emsol_p_no_p'
      call ems_se_ml_p_no_p(rt_cod, is)
      if (rt_cod .ne. 0) goto 7000
      call ems_se_rsmi_p_no_p(rt_cod)
      if (rt_cod .ne. 0) goto 7000
      call ems_se_inv_p_no_p(rt_cod)
      if (rt_cod .ne. 0) goto 7000
 7000 continue
      return
      end

ext: f

max_stars_repo_name: ERGO-Code/EMSOL

max_issues_repo_name: ERGO-Code/EMSOL

max_forks_repo_name: ERGO-Code/EMSOL


text:
      program simple_iri_driver
          implicit none

          logical jf(50),print_output
          real alati,along,dhour,heibeg,heiend,heistp
          real oarr(100)
          real outf(20,1000)
          integer jmag,k
          integer mmdd,iyyyy

          !call read_ig_rz
          !call readapf107

          print_output = .true.

          do k=1,50
              jf(k) = .true.
          enddo

          jf(5)  = .false.
          jf(6)  = .false.
          jf(21) = .false.
          jf(22) = .false.
          jf(23) = .false.
          jf(29) = .false.
          jf(30) = .false.

          jmag = 0
          alati = 45.
          along = -80.
          iyyyy = 2002
          mmdd = -10
          dhour = 12.
          heibeg = 250.
          heiend = 251.
          heistp = 1.

          call read_ig_rz
          call readapf107

          call iri_sub(jf,jmag,alati,along,iyyyy,mmdd,dhour,
     &               heibeg,heiend,heistp,outf,oarr)

          if (print_output) then
              print *,'=============='
              print *,'  Outputs 1'
              print *,'=============='

              print *,'electron density, outf(1,*)     = ',outf(1,1)
              print *,'neutral temperature, outf(2,*)  = ',outf(2,1)
              print *,'ion temperature, outf(3,*)      = ',outf(3,1)
              print *,'electron temperature, outf(4,*) = ',outf(4,1)
              print *,'O+ ion density, outf(5,*)       = ',outf(5,1)
              print *,'H+ ion density, outf(6,*)       = ',outf(6,1)
              print *,'HE+ ion density, outf(7,*)      = ',outf(7,1)
              print *,'O2+ ion density, outf(8,*)      = ',outf(8,1)
              print *,'NO+ ion density, outf(9,*)      = ',outf(9,1)
          endif

          call iri_sub(jf,jmag,alati,along,iyyyy,mmdd,dhour,
     &               heibeg,heiend,heistp,outf,oarr)

          if (print_output) then
              print *,'=============='
              print *,'  Outputs 2'
              print *,'=============='

              print *,'electron density, outf(1,*)     = ',outf(1,1)
              print *,'neutral temperature, outf(2,*)  = ',outf(2,1)
              print *,'ion temperature, outf(3,*)      = ',outf(3,1)
              print *,'electron temperature, outf(4,*) = ',outf(4,1)
              print *,'O+ ion density, outf(5,*)       = ',outf(5,1)
              print *,'H+ ion density, outf(6,*)       = ',outf(6,1)
              print *,'HE+ ion density, outf(7,*)      = ',outf(7,1)
              print *,'O2+ ion density, outf(8,*)      = ',outf(8,1)
              print *,'NO+ ion density, outf(9,*)      = ',outf(9,1)
          endif

      end program simple_iri_driver

ext: for


text: #redirect V.I.P. Studios Photography

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Exclude DN ranges from histogram, as specified.
C
      SUBROUTINE EXCLUDE_HISTOGRAM(FORMAT,hist,npts)
      IMPLICIT NONE
      CHARACTER*8 FORMAT
      INTEGER*4 HIST(-32768:32767)
      INTEGER*4 NPTS			!Number of pixels in histogram

      COMMON/C3/INMIN,INMAX,DNMIN,DNMAX
      INTEGER*4 INMIN,INMAX,DNMIN,DNMAX

      INTEGER*4 I,IDN,NPAIR,IXL,IXH,ICOUNT,IDEF
      REAL*4 TMEAN,TSIGMA,CUT

      REAL*4 RPARM(1000)
      INTEGER*4 IPARM(1000)
      EQUIVALENCE (RPARM,IPARM)

      CHARACTER*80 PRT
      LOGICAL*4 XVPTST
   20 FORMAT('DN values',I7,' thru',I7,' excluded')
   25 FORMAT('DN values',I7,' and',I7,' excluded')
   30 FORMAT(I7,' Pixels Outside Range',I7,' to',I7)
   40 FORMAT('Histogram after exclusion: Mean=',F11.4,' Sigma=',F11.4)
   50 FORMAT('Histogram after cutting: Mean=',F11.4,' Sigma=',F11.4)

C     ....Default is to exclude 0 and the maximum dn.
      IF (.NOT. XVPTST('INCLUDE'))  THEN
         WRITE (PRT,25) DNMIN,DNMAX
         CALL XVMESSAGE(PRT,' ')
         HIST(DNMIN) = 0
         HIST(DNMAX) = 0
      ENDIF

      CALL XVPARM('REXCLUDE',IPARM,ICOUNT,IDEF,200)
      IF (ICOUNT.GT.0) THEN
         NPAIR=ICOUNT/2
         IF (2*NPAIR.NE.ICOUNT) GOTO 900
         CALL XVMESSAGE(' ',' ')
         DO I=1,ICOUNT,2
            IXL=IPARM(I)
            IXH=IPARM(I+1)
            IF (IXL.LT.INMIN) IXL=INMIN
            IF (IXH.GT.INMAX) IXH=INMAX
            IF (IXL.GT.IXH) GOTO 902
            WRITE (PRT,20) IXL,IXH
            CALL XVMESSAGE(PRT,' ')
            CALL ZIA(hist(ixl),IXH-IXL+1)	!zero out exluded entries
         END DO
      END IF

      CALL XVPARM('IEXCLUDE',IPARM,ICOUNT,IDEF,100)
      IF (ICOUNT.GT.0) THEN
         CALL XVMESSAGE(' ',' ')
         DO I=1,ICOUNT
            IDN=IPARM(I)
            IF (IDN.LT.INMIN) IDN=INMIN
            IF (IDN.GT.INMAX) IDN=INMAX
            HIST(IDN)=0				!zero out excluded entries
         END DO
      END IF

C     ....Calculate new mean and standard deviation
      CALL STATI(HIST,tmean,tsigma,npts)
      IF (NPTS.LE.0) GOTO 904
      CALL XVMESSAGE(' ',' ')
      WRITE (PRT,40) TMEAN,TSIGMA
      CALL XVMESSAGE(PRT,' ')

C     ....Exclude DNs with less than a specified percent of max DN level
      CALL XVPARM('CUT',RPARM,ICOUNT,IDEF,1)
      IF (ICOUNT.GT.0) THEN
         CUT=RPARM(1)
         IF (CUT.LT.0.0 .OR. CUT.GT.100.0) THEN
            CALL XVMESSAGE('INVALID CUT VALUE - RESET TO 5.0',' ')
            CUT=5.0
         END IF
         CALL FCLIP(HIST,INMIN,INMAX,CUT)
         CALL STATI(HIST,tmean,tsigma,npts)
         IF (NPTS.LE.0) GOTO 908
         CALL XVMESSAGE(' ',' ')
         WRITE (PRT,40) TMEAN,TSIGMA
         CALL XVMESSAGE(PRT,' ')
      END IF

C        Print number OF pixels outside expected range.
      RETURN

  900 CALL XVMESSAGE('***Invalid count for parameter REXCLUDE',' ')
      CALL ABEND
  902 CALL XVMESSAGE('***Error in rexclude range ***',' ')
      CALL ABEND
  904 CALL XVMESSAGE('***No pixels remain after exclusion',' ')
      CALL ABEND
  908 CALL XVMESSAGE('***No pixels remain after cutting',' ')
      CALL ABEND
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Exclude all histogram entries with a frequency less than FCLIPC percent
C from the maximum frequency.
C
      SUBROUTINE FCLIP(HIST,INMIN,INMAX,FCLIPC)
      INTEGER*4 INMIN,INMAX
      INTEGER*4 HIST(-32768:32767)
      REAL      RTOL

C     ....Find max frequency
      MAX=0
      DO I=INMIN,INMAX
         IF (HIST(I).GT.MAX) MAX=HIST(I)
      END DO

C     ....Zero out any frequencies less than FCLIPC percent of max.
      RTOL = FLOAT(MAX)*(FCLIPC/100.0)
      DO I=INMIN,INMAX
         IF (HIST(I).LT.RTOL) HIST(I)=0
      END DO
      RETURN
      END

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
      SUBROUTINE F (W,WI,TP,XKZ,BZ,F0R,F0I,F1R,F1I,F2R,F2I)
C      DIV=2.6
C
C	THIS ROUTINE RETURNS ZETA*F0, ETC., FROM STIX FIRST BOOK,
C	iF0 is Fried and Conte Z
C
c****      DOUBLE PRECISION Z2R,Z2I,T,TR,TI,SR,SI,XN5,TTR,TTI
      DOUBLE PRECISION Z2R,Z2I,T,TR,TI,SR,SI,XN5,TTR,TTI
      PARAMETER(DIV=4.33)
C
      AR=W-XKZ*BZ
      B=XKZ*SQRT(2.*TP)
C****        IF (ABS(AR*AR-WI*WI)-34.*B*B) 100,100,105
      IF (ABS(AR*AR-WI*WI)-80.*B*B) 100,100,105
  100 ZR=AR/B
      ZI=WI/B
      TE=-ZI*ZI+ZR*ZR
      U=0.
      IF (TE.LT.80.) U=EXP(-TE)
C	print*,'f,100,tp,u',tp,u
      UI=-U*SIN(2.*ZR*ZI)*1.772453851
      UR=U*COS(2.*ZR*ZI)*1.772453851
      U=PR(UR,UI,ZR,ZI)
      UI=PI(UR,UI,ZR,ZI)
      UR=U
      U2R=TPR(W,WI,W,WI,UR,UI)/XKZ/XKZ
      U2I=TPI(W,WI,W,WI,UR,UI)/XKZ/XKZ
      GO TO 115
  105 IF ((AR*AR-WI*WI).LT.0..AND.WI.LT.0.) GO TO 110
      UR=0.
      UI=0.
      U2R=0.
      U2I=0.
C	print*,'f,105,tp,ur',tp,ur
      GO TO 115
C*****  110 U=EXP(28.)
  110 U=EXP(80.)
C      WRITE (7,180) AR,WI,W,TP,B
C      print*,'ar,wi,w,b', AR,WI,W,B
      UR=WI*U
      UI=UR
      U2R=UR
      U2I=UR
  115 IF (AR*AR+WI*WI-(DIV*B)**2) 120,120,135
C
C	ENTRY FOR POWER SERIES
C
  120 Z2R=ZR*ZR-ZI*ZI
      Z2I=2.*ZR*ZI
C	print*,'f,120,tp,z2r',tp,z2r
      TR=Z2R
      TI=Z2I
      SR=TR+SIGN(.5,XKZ)*UI
      SI=TI-SIGN(.5,XKZ)*UR
      DO 125 N=1,100
        XN5=DFLOAT(N)+.5D00
        T=-(TR*Z2R-TI*Z2I)/XN5
        TI=-(TR*Z2I+TI*Z2R)/XN5
        TR=T
        SR=SR+TR
        SI=SI+TI
C	PRINT*,'N,UR,UI',N,UR,UI
C	PRINT*,'SR,SI,TR,TI,',SR,SI,TR,TI
        IF (DABS(SR)+DABS(SI)-1.0D9*(DABS(TR)+DABS(TI))) 125,125,130
  125   CONTINUE
  130 F0R=-2.*SI
      F0I=2.*SR
      F1R=(PR(W,WI,F0R,F0I)+WI)/XKZ
      F1I=(PI(W,WI,F0R,F0I)-AR)/XKZ
      F2R=(PR(W,WI,F1R,F1I)+BZ*WI)/XKZ
      F2I=(PI(W,WI,F1R,F1I)-BZ*AR)/XKZ
C	WRITE(7,*) 's',w,wi,tp,F0R,F0I,F1R,F1I,F2R,F2I
      RETURN
C
C	ENTRY FOR ASYMPTOTIC SERIES
C
  135 ZR=QR(B,0.,AR,WI)
      ZI=QI(B,0.,AR,WI)
C	print*,'f,135,tp',tp
      YR=ZR*ZR-ZI*ZI
      YI=2.*ZR*ZI
      SR=QR(W,WI,AR,WI)
      SI=QI(W,WI,AR,WI)
      TR=-PI(SR,SI,SR,SI)*TP
      TI=PR(SR,SI,SR,SI)*TP
      IF (WI) 140,155,140
  140 IF (XKZ*WI) 145,145,150
  145 C=SIGN(2.,XKZ)
      GO TO 160
  150 C=0.
      GO TO 160
  155 C=SIGN(1.,XKZ)
  160 SR=TR+C*U2R
      SI=TI+C*U2I+BZ*BZ
C	write(7,*) 'xkz,wi,c,u',xkz,wi,c,u
      DO 170 N=2,30
        XN5=DFLOAT(N)-.5D00
        TTR=XN5*(TR*YR-TI*YI)
        TTI=XN5*(TR*YI+TI*YR)
C		EXIT WHEN TERMS START TO GET LARGER
       IF (DABS(TR)+DABS(TI)-DABS(TTR)-DABS(TTI)) 175,165,165
  165   TR=TTR
        TI=TTI
        SR=SR+TR
        SI=SI+TI
C		EXIT IF TERMS ARE SMALL ENOUGH
        IF (DABS(SR)+DABS(SI)-1.0D+10*(DABS(TR)+DABS(TI))) 170,175,175
  170   CONTINUE
C      WRITE (7,180) SR,SI
  175 F2R=SR
      F2I=SI
      F1R=QR(XKZ*F2R-BZ*WI,XKZ*F2I+BZ*AR,W,WI)
      F1I=QI(XKZ*F2R-BZ*WI,XKZ*F2I+BZ*AR,W,WI)
      F0R=QR(XKZ*F1R-WI,XKZ*F1I+AR,W,WI)
      F0I=QI(XKZ*F1R-WI,XKZ*F1I+AR,W,WI)
C	WRITE(7,*) 'a',w,wi,tp,F0R,F0I,F1R,F1I,F2R,F2I
C      WRITE(7,*) F0R,F0I,F1R,F1I,F2R,F2I
	RETURN
C
  180 FORMAT (1X,6E16.8)
      END
      SUBROUTINE Fold(W,WI,TP,XKZ,BZ,F0R,F0I,F1R,F1I,F2R,F2I)
      DIV=2.6
      AR=W-XKZ*BZ
      B=XKZ*SQRT(2.*TP)
        IF (ABS(AR*AR-WI*WI)-34.*B*B) 100,100,105
  100 ZR=AR/B
      ZI=WI/B
      TE=-ZI*ZI+ZR*ZR
      U=0.
      IF (TE.LT.34.) U=EXP(-TE)
      UI=-U*SIN(2.*ZR*ZI)*1.772453851
      UR=U*COS(2.*ZR*ZI)*1.772453851
      U=PR(UR,UI,ZR,ZI)
      UI=PI(UR,UI,ZR,ZI)
      UR=U
      U2R=TPR(W,WI,W,WI,UR,UI)/XKZ/XKZ
      U2I=TPI(W,WI,W,WI,UR,UI)/XKZ/XKZ
      GO TO 115
  105 IF ((AR*AR-WI*WI).LT.0..AND.WI.LT.0.) GO TO 110
      UR=0.
      UI=0.
      U2R=0.
      U2I=0.
      GO TO 115
  110 U=EXP(28.)
      WRITE (7,180) AR,WI,W,TP,B
      UR=WI*U
      UI=UR
      U2R=UR
      U2I=UR
  115 IF (AR*AR+WI*WI-(DIV*B)**2) 120,120,135
  120 Z2R=ZR*ZR-ZI*ZI
      Z2I=2.*ZR*ZI
      TR=Z2R
      TI=Z2I
      SR=TR+SIGN(.5,XKZ)*UI
      SI=TI-SIGN(.5,XKZ)*UR
      DO 125 N=1,40
        XN=N
        T=-(TR*Z2R-TI*Z2I)/(XN+.5)
        TI=-(TR*Z2I+TI*Z2R)/(XN+.5)
        TR=T
        SR=SR+TR
        SI=SI+TI
        IF (ABS(SR)+ABS(SI)-1.0E+6*(ABS(TR)+ABS(TI))) 125,125,130
  125   CONTINUE
  130 F0R=-2.*SI
      F0I=2.*SR
      F1R=(PR(W,WI,F0R,F0I)+WI)/XKZ
      F1I=(PI(W,WI,F0R,F0I)-AR)/XKZ
      F2R=(PR(W,WI,F1R,F1I)+BZ*WI)/XKZ
      F2I=(PI(W,WI,F1R,F1I)-BZ*AR)/XKZ
C	WRITE(7,*) F0R,F0I,F1R,F1I,F2R,F2I
      RETURN
  135 ZR=QR(B,0.,AR,WI)
      ZI=QI(B,0.,AR,WI)
      YR=ZR*ZR-ZI*ZI
      YI=2.*ZR*ZI
      SR=QR(W,WI,AR,WI)
      SI=QI(W,WI,AR,WI)
      TR=-PI(SR,SI,SR,SI)*TP
      TI=PR(SR,SI,SR,SI)*TP
      IF (WI) 140,155,140
  140 IF (XKZ*WI) 145,145,150
  145 C=SIGN(2.,XKZ)
      GO TO 160
  150 C=0.
      GO TO 160
  155 C=SIGN(1.,XKZ)
  160 SR=TR+C*U2R
      SI=TI+C*U2I+BZ*BZ
      DO 170 N=2,30
        XN=N
        TTR=(XN-.5)*(TR*YR-TI*YI)
        TTI=(XN-.5)*(TR*YI+TI*YR)
        IF (ABS(TR)+ABS(TI)-ABS(TTR)-ABS(TTI)) 175,165,165
  165   TR=TTR
        TI=TTI
        SR=SR+TR
        SI=SI+TI
        IF (ABS(SR)+ABS(SI)-1.0E+6*(ABS(TR)+ABS(TI))) 170,175,175
  170   CONTINUE
C      WRITE (7,180) SR,SI
  175 F2R=SR
      F2I=SI
      F1R=QR(XKZ*F2R-BZ*WI,XKZ*F2I+BZ*AR,W,WI)
      F1I=QI(XKZ*F2R-BZ*WI,XKZ*F2I+BZ*AR,W,WI)
      F0R=QR(XKZ*F1R-WI,XKZ*F1I+AR,W,WI)
      F0I=QI(XKZ*F1R-WI,XKZ*F1I+AR,W,WI)
C      WRITE(7,*) F0R,F0I,F1R,F1I,F2R,F2I
	RETURN
C
  180 FORMAT (1X,6E16.8)
      END
      FUNCTION PR (A,B,C,D)
      PR=A*C-B*D
      RETURN
      END
      FUNCTION PI (A,B,C,D)
      PI=A*D+B*C
      RETURN
      END
      FUNCTION QR (A,B,C,D)
      CD = AMAX1(ABS(C),ABS(D))
      AT = A/CD
      BT = B/CD
      CT = C/CD
      DT = D/CD
      QR = (AT*CT+BT*DT)/(CT*CT+DT*DT)
      RETURN
      END
      FUNCTION QI (A,B,C,D)
      CD = AMAX1(ABS(C),ABS(D))
      AT = A/CD
      BT = B/CD
      CT = C/CD
      DT = D/CD
      QI=(BT*CT-AT*DT)/(CT*CT+DT*DT)
      RETURN
      END
      FUNCTION TPR (A,B,C,D,E,F)
      TPR=E*(A*C-B*D)-F*(B*C+A*D)
      RETURN
      END
      FUNCTION TPI (A,B,C,D,E,F)
      TPI=F*(A*C-B*D)+E*(B*C+A*D)
      RETURN
      END


text:
      program main

c*********************************************************************72
c
cc MAIN runs an example of Dijkstra's minimum distance algorithm.
c
c  Discussion:
c
c    Given the distance matrix that defines a graph, we seek a list
c    of the minimum distances between node 0 and all other nodes.
c
c    This program sets up a small example problem and solves it.
c
c    The correct minimum distances are:
c
c      0   35   15   45   49   41
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    BY Sourangshu Ghosh
c
c  Author:
c
c  Sourangshu Ghosh
c
      implicit none

      include 'omp_lib.h'

      integer nv
      parameter ( nv = 6 )

      integer i
      integer i4_huge
      parameter ( i4_huge = 2147483647 )
      integer j
      integer mind(nv)
      integer ohd(nv,nv)

      call timestamp ( )
      write ( *, '(a)' )  ' '
      write ( *, '(a)' )  'DIJKSTRA_OPENMP:'
      write ( *, '(a)' )  '  FORTRAN77 version'
      write ( *, '(a)' )
     &  '  Use Dijkstra''s algorithm to determine the minimum'
      write ( *, '(a)' )
     &  '  distance from node 1 to each node in a graph,'
      write ( *, '(a)' )
     &  '  given the distances between each pair of nodes.'
      write ( *, '(a)' )  ' '
      write ( *, '(a)' )
     &  '  Although a very small example is considered, we'
      write ( *, '(a)' )  '  demonstrate the use of OpenMP directives'
      write ( *, '(a)' )  '  for parallel execution.'
c
c  Initialize the problem data.
c
      call init ( nv, ohd )
c
c  Print the distance matrix.
c
      write ( *, '(a)' )  ' '
      write ( *, '(a)' )  '  Distance matrix:'
      write ( *, '(a)' )  ' '
      do i = 1, nv
        write ( *, '(6(2x,i3))' ) ( ohd(i,j), j = 1, nv )
      end do
c
c  Carry out the algorithm.
c
      call dijkstra_distance ( nv, ohd, mind )
c
c  Print the results.
c
      write ( *, '(a)' )  ' '
      write ( *, '(a)' )  '  Minimum distances from node 1:'
      write ( *, '(a)' )  ' '
      do i = 1, nv
        write ( *, '(2x,i2,2x,i2)' ) i, mind(i)
      end do
c
c  Terminate.
c
      write ( *, '(a)' )  ' '
      write ( *, '(a)' )  'DIJKSTRA_OPENMP:'
      write ( *, '(a)' )  '  Normal end of execution.'

      write ( *, '(a)' )  ' '
      call timestamp ( )

      stop
      end
      subroutine dijkstra_distance ( nv, ohd, mind )

c*********************************************************************72
c
cc DIJKSTRA_DISTANCE uses Dijkstra's minimum distance algorithm.
c
c  Discussion:
c
c    We essentially build a tree.  We start with only node 0 connected
c    to the tree, and this is indicated by setting CONNECTED(0) = TRUE.
c
c    We initialize MIND(I) to the one step distance from node 0 to node I.
c
c    Now we search among the unconnected nodes for the node MV whose minimum
c    distance is smallest, and connect it to the tree.  For each remaining
c    unconnected node I, we check to see whether the distance from 0 to MV
c    to I is less than that recorded in MIND(I), and if so, we can reduce
c    the distance.
c
c    After NV-1 steps, we have connected all the nodes to 0, and computed
c    the correct minimum distances.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    By Sourangshu Ghosh
c
c  Author:
c
c    Sourangshu Ghosh
c
c  Parameters:
c
c    Input, integer NV, the number of nodes.
c
c    Input, integer OHD(NV,NV), the distance of the direct
c    link between nodes I and J.
c
c    Output, integer MIND(NV), the minimum
c    distance from node 1 to each node.
c
      implicit none

      include 'omp_lib.h'

      integer nv

      logical connected(nv)
      integer i
      integer i4_huge
      parameter ( i4_huge = 2147483647 )
      integer md
      integer mind(nv)
      integer mv
      integer my_first
      integer my_id
      integer my_last
      integer my_md
      integer my_mv
      integer my_step
      integer nth
      integer ohd(nv,nv)
c
c  Start out with only node 1 connected to the tree.
c
      connected(1) = .true.
      do i = 2, nv
        connected(i) = .false.
      end do
c
c  Initialize the minimum distance to the one-step distance.
c
      do i = 1, nv
        mind(i) = ohd(1,i)
      end do
c
c  Begin the parallel region.
c
c$omp parallel private ( my_first, my_id, my_last, my_md, my_mv, my_step )
c$omp&         shared ( connected, md, mind, mv, nth, ohd )

      my_id = omp_get_thread_num ( )
      nth = omp_get_num_threads ( )
      my_first =   (   my_id       * nv ) / nth + 1
      my_last  =   ( ( my_id + 1 ) * nv ) / nth
c
c  The SINGLE directive means that the block is to be executed by only
c  one thread, and that thread will be whichever one gets here first.
c
c$omp single
      write ( *, '(a)' ) ' '
      write ( *, * ) '  P', my_id,
     &  ': Parallel region begins with ', nth, ' threads.'
      write ( *, '(a)' ) ' '
c$omp end single
      write ( *, * ) '  P', my_id, ':  First=', my_first,
     &  '  Last=', my_last
c
c  Attach one more node on each iteration.
c
      do my_step = 2, nv
c
c  Before we compare the results of each thread, set the shared variable
c  MD to a big value.  Only one thread needs to do this.
c
c$omp single
        md = i4_huge
        mv = -1
c$omp end single
c
c  Each thread finds the nearest unconnected node in its part of the graph.
c  Some threads might have no unconnected nodes left.
c
        call find_nearest ( my_first, my_last, nv, mind, connected,
     &    my_md, my_mv )
c
c  In order to determine the minimum of all the MY_MD's, we must insist
c  that only one thread at a time execute this block!
c
c$omp critical
        if ( my_md .lt. md ) then
          md = my_md
          mv = my_mv
        end if
c$omp end critical
c
c  This barrier means that ALL threads have executed the critical
c  block, and therefore MD and MV have the correct value.  Only then
c  can we proceed.
c
c$omp barrier
c
c  If MV is -1, then NO thread found an unconnected node, so we're done early.
c  OpenMP does not like to BREAK out of a parallel region, so we'll just have
c  to let the iteration run to the end, while we avoid doing any more updates.
c
c  Otherwise, we connect the nearest node.
c
c$omp single
        if ( mv .ne. - 1 ) then
          connected(mv) = .true.
          write ( *, * ) '  P', my_id, ': Connecting node ', mv
        else
          write ( *, * )
     &      '  P', my_id, ': No connecting node on step ', my_step
        end if
c$omp end single
c
c  Again, we don't want any thread to proceed until the value of
c  CONNECTED is updated.
c
c$omp barrier
c
c  Now each thread should update its portion of the MIND vector,
c  by checking to see whether the trip from 0 to MV plus the step
c  from MV to a node is closer than the current record.
c
        if ( mv .ne. -1 ) then
          call update_mind ( my_first, my_last, nv, connected, ohd,
     &      mv, mind )
        end if
c
c  Before starting the next step of the iteration, we need all threads
c  to complete the updating, so we set a BARRIER here.
c
c$omp barrier

      end do
c
c  Once all the nodes have been connected, we can exit.
c
c$omp single
      write ( *, * ) ' '
      write ( *, * ) '  P', my_id, ': Exiting parallel region.'
c$omp end single

c$omp end parallel

      return
      end
      subroutine find_nearest ( s, e, nv, mind, connected, d, v )

c*********************************************************************72
c
cc FIND_NEAREST finds the nearest unconnected node.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c   By Sourangshu Ghosh
c
c  Author:
c
c    Sourangshu Ghosh
c
c  Parameters:
c
c    Input, integer S, E, the first and last nodes that
c    are to be checked.
c
c    Input, integer NV, the number of nodes.
c
c    Input, integer MIND(NV), the currently computed minimum
c    distance from node 1 to each node.
c
c    Input, logical CONNECTED(NV), is true for each connected node, whose
c    minimum distance to node 1 has been determined.
c
c    Output, integer D, the distance from node 1 to the nearest
c    unconnected node in the range S to E.
c
c    Output, integer V, the index of the nearest unconnected node
c    in the range S to E.
c
      implicit none

      integer nv

      logical connected(nv)
      integer d
      integer e
      integer i
      integer i4_huge
      parameter ( i4_huge = 2147483647 )
      integer mind(nv)
      integer s
      integer v

      d = i4_huge
      v = -1

      do i = s, e
        if ( .not. connected(i) .and. mind(i) .lt. d ) then
          d = mind(i)
          v = i
        end if
      end do

      return
      end
      subroutine init ( nv, ohd )

c*********************************************************************72
c
cc INIT initializes the problem data.
c
c  Discussion:
c
c    The graph uses 6 nodes, and has the following diagram and
c    distance matrix:
c
c    N0--15--N2-100--N3           0   40   15  Inf  Inf  Inf
c      \      |     /            40    0   20   10   25    6
c       \     |    /             15   20    0  100  Inf  Inf
c        40  20  10             Inf   10  100    0  Inf  Inf
c          \  |  /              Inf   25  Inf  Inf    0    8
c           \ | /               Inf    6  Inf  Inf    8    0
c            N1
c            / \
c           /   \
c          6    25
c         /       \
c        /         \
c      N5----8-----N4
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    By Sourangshu Ghosh
c
c  Author:
c
c   Sourangshu Ghosh
c
c  Parameters:
c
c    Input, integer NV, the number of nodes.
c
c    Output, integer OHD(NV,NV), the distance of the direct
c    link between nodes I and J.
c
      implicit none

      integer nv

      integer i
      integer i4_huge
      parameter ( i4_huge = 2147483647 )
      integer j
      integer ohd(nv,nv)

      do j = 1, nv
        do i = 1, nv
          ohd(i,j) = i4_huge
        end do
      end do

      do i = 1, nv
        ohd(i,i) = 0
      end do

      ohd(1,2) = 40
      ohd(1,3) = 15
      ohd(2,3) = 20
      ohd(2,4) = 10
      ohd(2,5) = 25
      ohd(3,4) = 100
      ohd(2,6) = 6
      ohd(5,6) = 8

      ohd(2,1) = 40
      ohd(3,1) = 15
      ohd(3,2) = 20
      ohd(4,2) = 10
      ohd(5,2) = 25
      ohd(4,3) = 100
      ohd(6,2) = 6
      ohd(6,5) = 8

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc TIMESTAMP prints out the current YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    By Sourangshu Ghosh
c
c  Author:
c
c    Sourangshu Ghosh
c
c  Parameters:
c
c    None
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ',
     &  'May      ', 'June     ', 'July     ', 'August   ',
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 0 ) then
            ampm = 'Midnight'
          else
            ampm = 'AM'
          end if
        end if
      end if

      write ( *,
     &  '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' )
     &  d, month(m), y, h, ':', n, ':', s, '.', mm, ampm

      return
      end
      subroutine update_mind ( s, e, nv, connected, ohd, mv, mind )

c*********************************************************************72
c
cc UPDATE_MIND updates the minimum distance vector.
c
c  Discussion:
c
c    We've just determined the minimum distance to node MV.
c
c    For each node I which is not connected yet,
c    check whether the route from node 0 to MV to I is shorter
c    than the currently known minimum distance.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    By Sourangshu
c
c  Author:
c
c    Sourangshu Ghosh
c
c  Parameters:
c
c    Input, integer S, E, the first and last nodes that are
c    to be checked.
c
c    Input, integer NV, the number of nodes.
c
c    Input, logical CONNECTED(NV), is true for each connected node, whose
c    minimum distance to node 0 has been determined.
c
c    Input, integer OHD(NV,NV), the distance of the direct link
c    between nodes I and J.
c
c    Input, integer MV, the node whose minimum distance to node 20
c    has just been determined.
c
c    Input/output, integer MIND(NV), the currently computed
c    minimum distances from node 1 to each node.  On output, the values for
c    nodes S through E have been updated.
c
      implicit none

      integer nv

      logical connected(nv)
      integer e
      integer i
      integer i4_huge
      parameter ( i4_huge = 2147483647 )
      integer mind(nv)
      integer mv
      integer ohd(nv,nv)
      integer s

      do i = s, e
        if ( .not. connected(i) ) then
          if ( ohd(mv,i) .lt. i4_huge ) then
            if ( mind(mv) + ohd(mv,i) .lt. mind(i) ) then
              mind(i) = mind(mv) + ohd(mv,i)
            end if
          end if
        end if
      end do

      return
      end


text:
Biotechnology involves the study of the use of living things to produce products or methods for human use. Biotech has been around for thousands of years with the use of yeast for cheeses, bread, and brewing booze. Modern biotechnology often uses genetic engineering to manipulate the genes within organisms to produce many different pharmaceuticals. Research in biotechnology has been used to produce everything from insulin to a fish that glows in the dark.  A local genetically engineered vegetable is the Flavr Savr Tomatoes Flavr Savr Tomato.

Classes/Major and Programs

There is a Biotechnology Major at UC Davis.  There are also related grad school disciplines such as the Biotechnology emphasis for the Microbiology graduate group, and Biochemistry, a specialization of the Chemistry doctoral program.

There is also a Biotechnology Program, created by Dean Hess in 1988.  In 1993, they decided that oversight would be a joint effort by: OVCR; CA&ES; College of Engineering; College of Letters & Sciences; Division of Biological Sciences; School of Medicine and School of Veterinary Medicine.  The program is a unit of the Office of Research.

http://www.news.ucdavis.edu/search/news_detail.lasso?id5940&titleUC%20Biotech%20Program%20Moves%20to%20UC%20Davis In 2001, the University of California Systemwide Biotechnology Research and Teaching Program moved to Davis from UC Berkeley, making Davis the home for lots of biotechnology research.  This helped get the $95 million needed to make the Genome Center.

Local Corporate Facilities

Quite a bit of Research & Development goes on in Davis.  Given the local agricultural emphasis, some of it is directed toward biotechnology in the ag sector.  For instance, Arcadia Biosciences develops agricultural products, while Novozymes, Inc. is a local company that is the the R&D arm of Danish company Novozymes.

Calgene LLC was founded in a local garage in 1980, and remained based out of Davis until they were bought by Monsanto in 1996.  They developed the genetically modified Flavr Savr Tomatoes Flavr Savr and ran subsidiaries across the country.

Zoogen specializes in sex identification of birds but also offers a range of genetic testing services for several species of mammals.

Seminis Seminis Vegetable Seed Company  the largest developer, grower and marketer of fruit and vegetable seeds in the world  has their second largest lab located in nearby Woodland.  Seminis was bought by http://en.wikipedia.org/wiki/Monsanto Monsanto in 2005, making Monsanto the largest GM and conventional seed company in the world.

Biotech Industry Funding at UC Davis

Third in a series of Seeds of Doubt 2004 Sacramento Bee articles on Biotechnology looks athttp://www.sacbee.com/static/live/news/projects/biotech/c3_1.html UC Davis research fundiing.

 Legal Matters

Family farmers are http://www.osgata.org/osgatapressreleases complaining against Monsantos gmos.

http://www.youtube.com/watch?vryiu5nd3Ih4&featureshare India sues Monsanto for piracy  youtube video


text:
SUBROUTINE elemnt(TASK, name, deltc, ttime, istop, igrav, ivect, lnod,   &
                  flag1, flag2)
!********************************************************************
! Standard ELEMENT routine
!********************************************************************
USE ctrl_db, ONLY: npoin, top, bottom
USE esets_db,ONLY: esets, eset, nelms
USE npo_db
IMPLICIT NONE

  !--- Dummy variables
  CHARACTER(len=*),INTENT(IN):: TASK
    !Optional parameters
  LOGICAL,OPTIONAL:: flag1, flag2
  CHARACTER(len=*),OPTIONAL:: name
  INTEGER(kind=4),OPTIONAL:: istop, igrav, ivect(:)
  INTEGER (kind=4), POINTER, OPTIONAL :: lnod(:,:)
  REAL(kind=8),OPTIONAL:: deltc, ttime
  !--- Local variables
  INTEGER(kind=4):: iset, i

  INTERFACE
    INCLUDE 'elemt1.h'
    INCLUDE 'elemt2.h'
    INCLUDE 'elem04.h'
    INCLUDE 'elem05.h'
    INCLUDE 'elemt6.h'
    INCLUDE 'elemt7.h'
    INCLUDE 'elemt8.h'
    INCLUDE 'elemt9.h'
    INCLUDE 'elem10.h'
    INCLUDE 'elem11.h'
    INCLUDE 'elem13.h'
    INCLUDE 'elem14.h'
    INCLUDE 'elem15.h'
    INCLUDE 'elem16.h'
    INCLUDE 'elem17.h'
    INCLUDE 'elem18.h'
    INCLUDE 'elem19.h'
    INCLUDE 'elem20.h'
    INCLUDE 'elem21.h'
    INCLUDE 'elem22.h'
    INCLUDE 'elem24.h'
    INCLUDE 'elem25.h'
  END INTERFACE

  IF (TRIM(task) == 'RESVPL') THEN
    IF (bottom) coorb=0d0
    IF (top) coort=0d0
    IF (top .OR. bottom) ifact=0
  END IF

  IF (PRESENT(flag2) .AND. TASK /= 'OUTDYN') flag2=.FALSE.

  DO iset=1,esets
    SELECT CASE(eset(iset))
    CASE (1)
      CALL elemt1(TASK, nelms, name, deltc, ttime, flag1, flag2)
    CASE (2)
      CALL elemt2(TASK, nelms, name, deltc, ttime, istop, flag1, flag2)
    CASE (4)
      CALL elem04(TASK, nelms, name, deltc, ttime, istop, &
                  lnod, flag1, flag2)
    CASE (5)
      CALL elem05(TASK, nelms, name, deltc, ttime, istop, &
                  lnod, flag1, flag2)
    CASE (6)
      CALL elemt6(TASK, nelms, name, deltc, ttime, istop, &
                  flag1, flag2)
    CASE (7)
      CALL elemt7(TASK, nelms, name, deltc, ttime, istop, &
                  ivect, flag1, flag2)
    CASE (8)
      CALL elemt8(TASK, nelms, name, deltc, ttime, istop, &
                  flag1, flag2)
    CASE (9)
      CALL elemt9(TASK, nelms, name, deltc, ttime, istop, &
                  flag1, flag2)
    CASE (10)
      CALL elem10(TASK, nelms, name, istop, flag2)
    CASE (11)
      CALL elem11(TASK, nelms, name, deltc, ttime, istop, flag1, flag2)
    CASE (13)
      CALL elem13(TASK, nelms, name, deltc, ttime, istop, ivect, flag1, flag2)
    CASE (14)
      CALL elem14(TASK, nelms, name, deltc, ttime, istop, ivect, flag1, flag2)
    CASE (15)
      CALL elem15(TASK, deltc, ttime, istop, flag1, flag2)
    CASE (16)
      CALL elem16(TASK, nelms, name, deltc, ttime, istop, &
                  lnod, flag1, flag2)
    CASE (17)
      CALL elem17(TASK, nelms, name, deltc, ttime, istop, &
                  lnod, flag1, flag2)
    CASE (18)
      CALL elem18(TASK, nelms, name, deltc, ttime, istop, &
                  ivect, lnod, flag1, flag2)
    CASE (19)
      CALL elem19(TASK, name, deltc, ttime, istop, &
                  lnod, flag1, flag2)
    CASE (20)
      CALL elem20(TASK, name, deltc, ttime, istop, ivect, lnod, flag1, flag2)
    CASE (21)
      CALL elem21(TASK, deltc, ttime, istop)
    CASE (22)
      CALL elem22(TASK, name, deltc, ttime, istop, flag2)
    CASE (24)
      CALL elem24(TASK, name, deltc, ttime, istop, ivect,   &
                  flag1, flag2)
    CASE (25)
      CALL elem25(TASK, nelms, name, deltc, ttime, istop, &
                  lnod, flag1, flag2)
    END SELECT

    IF (TRIM(task) == 'BOUNDA' .OR.         &
        TRIM(task) == 'DELETE' .OR.         &
        TRIM(task) == 'INIGAU' .OR.         &
        TRIM(task) == 'NODSET' .OR.         &
        TRIM(task) == 'RFSTBC' .OR.         &
        TRIM(task) == 'SEARCH' .OR.         &
        TRIM(task) == 'SMOOTH' .OR.         &
        TRIM(task) == 'SURFAC' ) THEN
      IF (flag2) THEN
        IF (TRIM(task) == 'SEARCH')  igrav=eset(iset)   !=etype (type of element is 'sname')
        EXIT        !set found, exit loop
      END IF
    END IF

  END DO
  IF (TRIM(task) == 'RESVPL') THEN
    IF (bottom) THEN
      DO i=1,npoin
        IF (ifact(i) > 0) THEN
          coorb(:,i) = coorb(:,i)/ifact(i) + coora(:,i)
        ELSE
          coorb(:,i) = coora(:,i)
        END IF
      END DO
    END IF
    IF (top) THEN
      DO i=1,npoin
        IF (ifact(i) > 0) THEN
          coort(:,i) = coort(:,i)/ifact(i) + coora(:,i)
        ELSE
          coort(:,i) = coora(:,i)
        END IF
      END DO
    END IF
  END IF

RETURN
END SUBROUTINE elemnt

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


program test_pic_3d

use mesh_fields_m
use poisson_3d_m
use output_m
use particles_m
use compute_rho_cic_m
use interpolation_cic_m

implicit none

type(mesh_t)       :: mesh
type(fields_3d_t)  :: fields
type(poisson_3d_t) :: poisson
type(particles_t)  :: particles

integer,    parameter :: nx = 64, ny = 64, nz = 4
integer(8), parameter :: nbpart = nx * ny * 100
real(8)               :: xmin, xmax, ymin, ymax, zmin, zmax
real(8)               :: pi
integer               :: i, j, k, m
real(8)               :: x, y, z
real(8)               :: dx, dy, dz
real(8)               :: err_x, err_y, err_z

pi = 4d0 * atan(1d0)

xmin = 0d0
xmax = 18d0
ymin = 0d0
ymax = 18d0
zmin = 0d0
zmax = 1d0

call init_mesh( mesh, xmin, xmax, nx, ymin, ymax, ny, zmin, zmax, nz )

dx = mesh%dx
dy = mesh%dy
dz = mesh%dz

call init_fields( fields, mesh )

call init_poisson( poisson, mesh)

print*, " number of particles = ", nbpart

call init_particles_3d( particles, nbpart, mesh )

call compute_rho_cic( fields, particles)

print*, sum(fields%rho) * dx * dy * dz, nbpart * particles%w

call solve_poisson( poisson, fields)

call write_data( 1, fields)

do k = 1, nz+1
    z = (k-1) * dz
    do j = 1, ny+1
        y = (j-1) * dy
        do i = 1, nx+1
            x = (i-1) * dx
            fields%e(1,i,j,k) = sin(x)
            fields%e(2,i,j,k) = sin(y)
            fields%e(3,i,j,k) = sin(z)
        end do
    end do
end do

call interpolate_eb_cic( particles, fields)

err_x = 0d0
err_y = 0d0
err_z = 0d0

do m = 1, nbpart
    x = (i-1) * mesh%dx
    err_x = err_x + abs(particles%e(1,m) - sin(particles%x(1,m)))
    err_y = err_y + abs(particles%e(2,m) - sin(particles%x(2,m)))
    err_z = err_z + abs(particles%e(3,m) - sin(particles%x(3,m)))
end do

print*, err_x/nbpart, err_y/nbpart, err_z/nbpart

end


text:
MODULE ewald_module

  use ewald_variables, only: eta,mg,mr,LG,LR,ipair,mpair
  use aa_module, only: aa, Va
  use bb_module, only: bb
  use atom_module, only: Natom, aa_atom, ki_atom
  use pseudopot_module, only: Zps
  use parallel_module
  use watch_module
  use hsort_module

  implicit none

  PRIVATE
  PUBLIC :: test_ewald, calc_ewald, cal_ewald

  real(8),allocatable :: zatom(:)
  real(8) :: Qtot,Qtot2,rrcut,ecut,Vcell
  integer,allocatable :: id(:),ir(:)
  real(8) :: ewldg_0=0.d0, g_0=0.d0
  real(8) :: ewldr_0=0.d0, r_0=0.d0

CONTAINS

  SUBROUTINE cal_ewald(Ewld)
    implicit none
    real(8),intent(OUT) :: Ewld
    real(8) :: ewldg,ewldr,pi
    pi=acos(-1.d0)
    call calc_ewldg(mg,LG,ewldg)
    ewldg=ewldg*4.d0*pi/Vcell
    ewldg=ewldg-pi/(eta*Vcell)*Qtot**2
    ewldg=ewldg-2.d0*sqrt(eta/pi)*Qtot2
    call calc_ewldr(mr,LR,ewldr)
    Ewld=0.5d0*(ewldg+ewldr)
  END SUBROUTINE cal_ewald

  SUBROUTINE test_ewald(Ewld)
    implicit none
    real(8),intent(OUT) :: Ewld
    integer,parameter :: maxloop=50,max_loop_r=10,max_loop_g=10
    integer :: loop_eta,loop_g,loop_r,mg0,mr0,i,m1,m2,m3,m,n,j
    integer :: mg_tot_0,mr_tot_0,mr_tot,mg_tot,mr_tmp,mr_dif,ierr
    integer,allocatable :: LR_tmp(:,:)
    real(8),parameter :: epsg=1.d-12, epsr=1.d-12, ep=1.d-10
    real(8),parameter :: factor=1.0d0
    real(8) :: const,const1,ewldr,ewldr0,ewldg,ewldg0,eta_small,eta_big
    real(8) :: r,r0,rr,alpha,g,g0,t,pi,max_bb,max_aa,eta_0,tt0,ttt,gg
    real(8) :: ewldr_tmp,gmax
    real(8) :: ct0,et0,ct1,et1,timg(2),timg0(2),timr(2),timr0(2)
    logical :: flag_conv
    integer,allocatable :: grd_chk(:,:,:),grd_chk0(:,:,:)

    pi=acos(-1.d0)

    Vcell = abs( aa(1,1)*aa(2,2)*aa(3,3)+aa(1,2)*aa(2,3)*aa(3,1) &
                +aa(1,3)*aa(2,1)*aa(3,2)-aa(1,3)*aa(2,2)*aa(3,1) &
                -aa(1,2)*aa(2,1)*aa(3,3)-aa(1,1)*aa(2,3)*aa(3,2) )

    max_aa=0.d0
    do i=1,3
       t=sqrt( sum(aa(1:3,i)**2) )
       max_aa=max(max_aa,t)
    end do
    max_bb=0.d0
    do i=1,3
       t=sqrt( sum(bb(1:3,i)**2) )
       max_bb=max(max_bb,t)
    end do

    if ( .not.allocated(zatom) ) then
       allocate( zatom(Natom) )
       do i=1,Natom
          zatom(i)=nint( Zps(ki_atom(i)) )
       end do
    end if

    if ( .not.allocated(id) ) then
       allocate( id(0:nprocs-1) )
       allocate( ir(0:nprocs-1) )
    end if

    if ( .not.allocated(ipair) ) then
       mpair = ( Natom*(Natom+1) )/2
       call prep_parallel(mpair,nprocs,ir,id)
       mpair = ir(myrank)
       allocate( ipair(2,mpair) )
       m=0
       n=0
       do j=1,Natom
       do i=j,Natom
          m=m+1
          if ( m < id(myrank)+1 .or. id(myrank)+ir(myrank) < m ) cycle
          n=n+1
          ipair(1,n)=i
          ipair(2,n)=j
       end do
       end do
    end if

    eta = pi/Vcell**(2.d0/3.d0)

    eta_small = 0.d0
    eta_big   = 0.d0
    eta_0     = 0.d0
    tt0       = 0.d0
    mg_tot_0  = 0
    mr_tot_0  = 0

    Qtot =0.d0
    Qtot2=0.d0
    do i=1,Natom
       Qtot =Qtot +zatom(i)
       Qtot2=Qtot2+zatom(i)**2
    end do

    const1 = 4.d0*Pi/Vcell

    do loop_eta=1,maxloop

!--- G

       const=1.d0/(4.d0*eta)
!       do i=1,10000
!          g=i
!          t=const1*exp(-const*g*g)/(g*g)
!          if ( t <= epsg ) exit
!       end do
       g=12.d0*sqrt(eta)

       ewldg0 = 0.d0
       flag_conv = .false.
       do loop_g=1,max_loop_g

          ecut=g*g

          call search_grid(bb,ecut,m1,m2,m3,mg_tot)
          call prep_parallel(mg_tot,nprocs,ir,id)
          mg=ir(myrank)
          if ( allocated(LG) ) deallocate( LG )
          allocate( LG(3,mg) )
          call construct_grid(bb,ecut,mg_tot,m1,m2,m3,mg,LG)

          call watch(ct0,et0)
          call calc_ewldg(mg,LG,ewldg)
          ewldg=ewldg*4.d0*pi/Vcell
          ewldg=ewldg-pi/(eta*Vcell)*Qtot**2
          ewldg=ewldg-2.d0*sqrt(eta/pi)*Qtot2
          call watch(ct1,et1) ; timg(1)=ct1-ct0 ; timg(2)=et1-et0

          ir(myrank)=mg
          call mpi_allgather &
               (mg,1,mpi_integer,ir,1,mpi_integer,mpi_comm_world,ierr)

          if ( disp_switch_parallel ) then
             write(*,'(1x,2i4,2i8,2g26.15,2g12.5)') &
                  loop_eta,loop_g,mg_tot,mg,ewldg,sqrt(ecut),timg(1:2)
          end if

          if ( abs((ewldg-ewldg0)/ewldg) < 1.d-15 ) then
             flag_conv=.true.
             ewldg=ewldg0
             mg=mg0
             g=g0
             timg(:)=timg0(:)
             exit
          else
             ewldg0=ewldg
             mg0=mg
             g0=g
             g=(g+max_bb)*factor
             timg0(:)=timg(:)
          end if

       end do ! loop_g

       if ( flag_conv ) then
          if ( mg <= 0 ) then
             gg=0.d0
          else
             gg=( bb(1,1)*LG(1,mg)+bb(1,2)*LG(2,mg)+bb(1,3)*LG(3,mg) )**2 &
               +( bb(2,1)*LG(1,mg)+bb(2,2)*LG(2,mg)+bb(2,3)*LG(3,mg) )**2 &
               +( bb(3,1)*LG(1,mg)+bb(3,2)*LG(2,mg)+bb(3,3)*LG(3,mg) )**2
          end if
          call mpi_allreduce(gg,ecut,1,mpi_real8,mpi_max,mpi_comm_world,i)
          call search_grid(bb,ecut,m1,m2,m3,mg_tot)
          call prep_parallel(mg_tot,nprocs,ir,id)
          mg=ir(myrank)
          if ( allocated(LG) ) deallocate( LG )
          allocate( LG(3,mg) )
          call construct_grid(bb,ecut,mg_tot,m1,m2,m3,mg,LG)
       else
          stop "ewldg is not converged"
       end if

       if ( disp_switch_parallel ) then
          write(*,'(1x,"ewald(G)  =",2g24.15,2i8)') ewldg*0.5d0,sqrt(ecut),mg,mg_tot
       end if

!--- R

       alpha = sqrt(eta)
!       do i=1,1000
!          r=i
!          t=erfc(alpha*r)/r
!          if ( t <= epsr ) exit
!       end do
       r=6.d0/sqrt(eta)

       ewldr0    = 0.d0
!       ewldr     = 0.d0
       ewldr     = ewldg ; ewldg=0.d0
       ewldr_tmp = ewldr
       flag_conv = .false.
       mr_tmp    = 0
       timr(:)   = 0.d0
       do loop_r=1,max_loop_r

          rrcut=r*r

          call watch(ct0,et0)

          call search_grid(aa,rrcut,m1,m2,m3,mr_tot)
!          call prep_parallel(mr_tot,nprocs,ir,id) ; mr=ir(myrank)
          mr=mr_tot ; id(:)=0 ; ir(:)=mr

          allocate( grd_chk(-m1:m1,-m2:m2,-m3:m3) ) ; grd_chk=0

          if ( loop_r == 1 ) then
          else
             allocate( grd_chk0(-m1:m1,-m2:m2,-m3:m3) ) ; grd_chk0=0
             mr_tmp = size(LR,2)
             do i=1,mr_tmp
                grd_chk0( LR(1,i),LR(2,i),LR(3,i) ) = 1
             end do
             m=size(grd_chk)
             call mpi_allreduce(grd_chk0,grd_chk,m &
                               ,mpi_integer,mpi_sum,mpi_comm_world,ierr)
             deallocate( grd_chk0 )
          end if

          if ( allocated(LR) ) deallocate(LR)
          allocate( LR(3,mr) )

          call construct_grid(aa,rrcut,mr_tot,m1,m2,m3,mr,LR)

          if ( allocated(LR_tmp) ) deallocate( LR_tmp )
          allocate( LR_tmp(3,mr) ) ; LR_tmp=0

          if ( loop_r == 1 ) then
             m=mr
             LR_tmp(1:3,1:mr) = LR(1:3,1:mr)
          else
             m=0
             do i=1,mr
                if ( grd_chk(LR(1,i),LR(2,i),LR(3,i)) == 0 ) then
                   m=m+1
                   LR_tmp(1:3,m) = LR(1:3,i)
                end if
             end do
          end if

          deallocate( grd_chk )

          call watch(ct1,et1)
!(1)
!          call watch(ct0,et0)
!          ewldr=ewldr_tmp
!!          ewldr=0.d0
!          call calc_ewldr(mr,LR,ewldr)
!          call watch(ct1,et1)
!          timr(1)=ct1-ct0 ; timr(2)=et1-et0
!(2)
          call watch(ct0,et0)
          call calc_ewldr(m,LR_tmp,ewldr)
          call watch(ct1,et1)
          timr(1)=timr(1)+ct1-ct0 ; timr(2)=timr(2)+et1-et0

          deallocate( LR_tmp )

          if ( disp_switch_parallel ) then
             write(*,'(1x,2i4,3i8,2g26.15,2g12.5)') &
                  loop_eta,loop_r,mr_tot,mr,m,ewldr,sqrt(rrcut),timr(1:2)
          end if

          if ( abs((ewldr-ewldr0)/ewldr) < 1.d-15 ) then
             flag_conv=.true.
             ewldr=ewldr0
             mr=mr0
             r=r0
             timr(:)=timr0(:)
             exit
          else
             ewldr0=ewldr
             mr0=mr
             r0=r
             r=(r+max_aa)*factor
             timr0(:)=timr(:)
          end if

       end do ! loop_r

       if ( flag_conv ) then
          rr=( aa(1,1)*LR(1,mr)+aa(1,2)*LR(2,mr)+aa(1,3)*LR(3,mr) )**2 &
            +( aa(2,1)*LR(1,mr)+aa(2,2)*LR(2,mr)+aa(2,3)*LR(3,mr) )**2 &
            +( aa(3,1)*LR(1,mr)+aa(3,2)*LR(2,mr)+aa(3,3)*LR(3,mr) )**2
          call mpi_allreduce(rr,rrcut,1,mpi_real8,mpi_max,mpi_comm_world,i)
          call search_grid(aa,rrcut,m1,m2,m3,mr_tot)
!          call prep_parallel(mr_tot,nprocs,ir,id) ; mr=ir(myrank)
          mr=mr_tot ; id(:)=0 ; ir(:)=mr
          if ( allocated(LR) ) deallocate( LR )
          allocate( LR(3,mr) )
          call construct_grid(aa,rrcut,mr_tot,m1,m2,m3,mr,LR)
       else
          stop "ewldr is not converged"
       end if

       if ( disp_switch_parallel ) then
          write(*,'(1x,"ewald(G+R)=",2g24.15,2i8)') ewldr*0.5d0,sqrt(rrcut),mr,mr_tot
       end if

       Ewld=0.5d0*(ewldr+ewldg)

       call mpi_allreduce(timg,timg0,2,mpi_real8,mpi_max,mpi_comm_world,i)
       call mpi_allreduce(timr,timr0,2,mpi_real8,mpi_max,mpi_comm_world,i)

       if ( disp_switch_parallel ) then
          write(*,'(1x,i6,4x,"Ewld=",3g22.15,2i8,2x,3g12.5)') &
               loop_eta,Ewld,eta,eta-eta_0,mg_tot,mr_tot &
               ,timg0(2),timr0(2),timg0(2)+timr0(2)
       end if

       if ( abs(eta-eta_0) < ep ) exit
       if ( mg_tot_0 == mg_tot .and. mr_tot_0 == mr_tot ) exit

       mg_tot_0 = mg_tot
       mr_tot_0 = mr_tot

       ttt=timg0(2)+timr0(2)
       if ( loop_eta == 1 ) then
          eta_0 = eta
          tt0 = ttt
          eta = eta*2.d0
       else
          if ( eta_big == 0.d0 ) then
             if ( ttt < tt0 ) then
                eta_small = eta_0
                eta_0 = eta
                tt0 = ttt
                eta = eta*2.d0
             else
                eta_big = eta
                if ( mod(loop_eta,2) == 0 ) then
                   eta=(eta_small+eta_0)/2
                else
                   eta=(eta_0+eta_big)/2
                end if
             end if
          else
             if ( mod(loop_eta,2) == 1 ) then
                if ( ttt > tt0 ) then
                   eta_small = eta
                   eta = (eta_0+eta_big)/2
                else
                   eta_big = eta_0
                   eta_0 = eta
                   tt0 = ttt
                   eta = (eta_0+eta_big)/2
                end if
             else
                if ( ttt > tt0 ) then
                   eta_big = eta
                   eta = (eta_small+eta_0)/2
                else
                   eta_small = eta_0
                   eta_0 = eta
                   tt0 = ttt
                   eta = (eta_small+eta_0)/2
                end if
             end if
          end if
       end if

       if ( loop_eta == maxloop ) eta=eta_0

    end do ! loop_eta

    return
  END SUBROUTINE test_ewald

  SUBROUTINE search_grid(v,vvcut,m1,m2,m3,mv)
    implicit none
    real(8),intent(IN)  :: v(3,3),vvcut
    integer,intent(OUT) :: m1,m2,m3,mv
    integer :: i,m,i1,i2,i3
    real(8) :: vx,vy,vz,vv
    m1=1
    m2=1
    m3=1
    mv=0
    do i=1,100000
       m=0
       do i3=-m3,m3
       do i2=-m2,m2
       do i1=-m1,m1
          vx=v(1,1)*i1+v(1,2)*i2+v(1,3)*i3
          vy=v(2,1)*i1+v(2,2)*i2+v(2,3)*i3
          vz=v(3,1)*i1+v(3,2)*i2+v(3,3)*i3
          vv=vx*vx+vy*vy+vz*vz
          if ( vv < vvcut + 1.d-8 ) then
             m=m+1
          end if
       end do
       end do
       end do
       if ( m == mv ) exit
       mv=m
       m1=m1+1
       m2=m2+1
       m3=m3+1
    end do
  END SUBROUTINE search_grid

  SUBROUTINE prep_parallel(mv,nprocs,ir,id)
    implicit none
    integer,intent(IN) :: mv,nprocs
    integer,intent(OUT) :: ir(0:nprocs-1),id(0:nprocs-1)
    integer :: i,j
    ir(:)=0
    do i=1,mv
       j=mod(i-1,nprocs)
       ir(j)=ir(j)+1
    end do
    do j=0,nprocs-1
       id(j)=sum(ir(0:j))-ir(j)
    end do
  END SUBROUTINE prep_parallel

  SUBROUTINE construct_grid(v,vvcut,mt,m1,m2,m3,mv,LV)
    implicit none
    real(8),intent(IN)  :: v(3,3),vvcut
    integer,intent(IN)  :: mt,m1,m2,m3,mv
    integer,intent(OUT) :: LV(3,mv)
    integer :: m,n,i1,i2,i3,i,j
    integer,allocatable :: indx(:),LVtmp(:,:)
    real(8) :: vx,vy,vz,vv
    real(8),allocatable :: ll(:)
    allocate( LVtmp(3,mt) ) ; LVtmp=0
    m=0
    n=0
    do i3=-m3,m3
    do i2=-m2,m2
    do i1=-m1,m1
       vx=v(1,1)*i1+v(1,2)*i2+v(1,3)*i3
       vy=v(2,1)*i1+v(2,2)*i2+v(2,3)*i3
       vz=v(3,1)*i1+v(3,2)*i2+v(3,3)*i3
       vv=vx*vx+vy*vy+vz*vz
       if ( vv < vvcut + 1.d-8 ) then
          m=m+1
          LVtmp(1,m)=i1
          LVtmp(2,m)=i2
          LVtmp(3,m)=i3
          if ( m < id(myrank)+1 .or. id(myrank)+ir(myrank) < m ) cycle
          n=n+1
          LV(1,n)=i1
          LV(2,n)=i2
          LV(3,n)=i3
       end if
    end do
    end do
    end do
    if ( n /= mv ) then
       write(*,*) "n,m,mv,mt=",n,m,mv,mt,myrank
       stop "stop at construct_grid(ewald1)"
    end if
    if ( m /= mt ) then
       write(*,*) "n,m,mv,mt=",n,m,mv,mt,myrank
       stop "stop at construct_grid(ewald2)"
    end if

    allocate( ll(m),indx(m) )
    do i=1,m
       i1=LVtmp(1,i)
       i2=LVtmp(2,i)
       i3=LVtmp(3,i)
       vx=v(1,1)*i1+v(1,2)*i2+v(1,3)*i3
       vy=v(2,1)*i1+v(2,2)*i2+v(2,3)*i3
       vz=v(3,1)*i1+v(3,2)*i2+v(3,3)*i3
       ll(i)=vx*vx+vy*vy+vz*vz
    end do
    call indexx(m,ll,indx)
    if ( m == n ) then
       do i=1,m
          LV(1:3,i)=LVtmp(1:3,indx(i))
       end do
    else
    n=0
    do i=1,m
       j=mod(i-1,nprocs)
       if ( j == myrank ) then
          n=n+1
          LV(1:3,n)=LVtmp(1:3,indx(i))
       end if
    end do
    end if
    deallocate( indx,ll,LVtmp )

  END SUBROUTINE construct_grid

  SUBROUTINE calc_ewldg(mg,LG,ewldg)
    implicit none
    integer,intent(INOUT) :: mg
    integer,intent(IN)    :: LG(3,mg)
    real(8),intent(OUT) :: ewldg
    integer :: i,a,m1
    real(8) :: x,y,z,t,pi2,gx,gy,gz,gg,e0,e1,g1,const
    complex(8) :: sg
    pi2=2.d0*acos(-1.d0)
    const=1.d0/(4.d0*eta)
    e0=0.d0
    e1=0.d0
    g1=0.d0
    m1=0
    do i=1,mg
       x=pi2*LG(1,i)
       y=pi2*LG(2,i)
       z=pi2*LG(3,i)
       if ( all( LG(1:3,i) == 0 ) ) cycle
       sg=(0.d0,0.d0)
!$OMP parallel do private( t ) reduction(+:sg)
       do a=1,Natom
          t=x*aa_atom(1,a)+y*aa_atom(2,a)+z*aa_atom(3,a)
          sg=sg+zatom(a)*dcmplx(cos(t),-sin(t))
       end do
!$OMP end parallel do
       gx=bb(1,1)*LG(1,i)+bb(1,2)*LG(2,i)+bb(1,3)*LG(3,i)
       gy=bb(2,1)*LG(1,i)+bb(2,2)*LG(2,i)+bb(2,3)*LG(3,i)
       gz=bb(3,1)*LG(1,i)+bb(3,2)*LG(2,i)+bb(3,3)*LG(3,i)
       gg=gx*gx+gy*gy+gz*gz
       t = exp(-const*gg)/gg
       e0 = e0 + abs(sg)**2*t
!       write(*,'(1x,i6,4g16.8)') i,gg,t,abs(sg)**2,e0
       if ( e0 /= e1 ) then
          e1=e0
          g1=gg
          m1=i
       end if
    end do
    call mpi_allreduce(e0,ewldg,1,mpi_real8,mpi_sum,mpi_comm_world,i)
    mg=m1
  END SUBROUTINE calc_ewldg

  SUBROUTINE calc_ewldr(mr,LR,ewldr)
    implicit none
    integer,intent(IN) :: mr,LR(3,mr)
    real(8),intent(INOUT) :: ewldr
    integer :: i,i1,i2,i1i2
    real(8) :: a1_2,a2_2,a3_2,a1_1,a2_1,a3_1,z1,z2
    real(8) :: sum0,sum1,sum2,a1,a2,a3,x,y,z,r,t,alpha
    alpha=sqrt(eta)
!    ewldr=0.d0
    sum2=0.d0
    do i1i2=1,mpair
       i1=ipair(1,i1i2)
       i2=ipair(2,i1i2)
       a1_2=aa_atom(1,i2)
       a2_2=aa_atom(2,i2)
       a3_2=aa_atom(3,i2)
       z2=zatom(i2)
       a1_1=aa_atom(1,i1)
       a2_1=aa_atom(2,i1)
       a3_1=aa_atom(3,i1)
       z1=zatom(i1)
       sum0=0.d0
!$OMP parallel do private( a1,a2,a3,x,y,z,r,t ) reduction(+:sum0)
       do i=1,mr
          a1=LR(1,i)+a1_1-a1_2
          a2=LR(2,i)+a2_1-a2_2
          a3=LR(3,i)+a3_1-a3_2
          x=aa(1,1)*a1+aa(1,2)*a2+aa(1,3)*a3
          y=aa(2,1)*a1+aa(2,2)*a2+aa(2,3)*a3
          z=aa(3,1)*a1+aa(3,2)*a2+aa(3,3)*a3
          r=sqrt(x*x+y*y+z*z)
          if ( r == 0.d0 ) cycle
          t=erfc(alpha*r)/r
          sum0=sum0+t
       end do
!$OMP end parallel do
       if ( i1 /= i2 ) sum0=sum0*2.d0
       sum2=sum2+z1*z2*sum0
    end do
    call mpi_allreduce(sum2,sum1,1,mpi_real8,mpi_sum,mpi_comm_world,i)
    ewldr=ewldr+sum1
  END SUBROUTINE calc_ewldr

  SUBROUTINE calc_ewald(Ewld)
    implicit none
    real(8),intent(OUT) :: Ewld
    integer,parameter :: maxloop=50,max_loop_r=10,max_loop_g=10
    integer :: i,m1,m2,m3,m,n,j,loop_g,loop_r
    integer :: mg_tot_0,mr_tot_0,mr_tot,mg_tot,mr_tmp,mr_dif,ierr
    integer,allocatable :: LR_tmp(:,:)
    real(8),parameter :: epsg=1.d-12, epsr=1.d-12, ep=1.d-10
    real(8),parameter :: factor=1.0d0
    real(8) :: const,const1,ewldr,ewldg
    real(8) :: r,rr,alpha,g,t,pi,max_bb,max_aa,tt0,ttt,gg
    real(8) :: ewldr_tmp,gmax
    real(8) :: ct0,et0,ct1,et1,timg(2),timg0(2),timr(2),timr0(2)
    integer,allocatable :: grd_chk(:,:,:),grd_chk0(:,:,:)
    logical :: disp_sw

    call check_disp_switch( disp_sw, 0 )
    call  write_border( 0, "calc_ewald(start)" )

    pi=acos(-1.d0)

    Vcell = abs( aa(1,1)*aa(2,2)*aa(3,3)+aa(1,2)*aa(2,3)*aa(3,1) &
                +aa(1,3)*aa(2,1)*aa(3,2)-aa(1,3)*aa(2,2)*aa(3,1) &
                -aa(1,2)*aa(2,1)*aa(3,3)-aa(1,1)*aa(2,3)*aa(3,2) )

    max_aa=0.d0
    do i=1,3
       t=sqrt( sum(aa(1:3,i)**2) )
       max_aa=max(max_aa,t)
    end do
    max_bb=0.d0
    do i=1,3
       t=sqrt( sum(bb(1:3,i)**2) )
       max_bb=max(max_bb,t)
    end do

    if ( .not.allocated(zatom) ) then
       allocate( zatom(Natom) )
       do i=1,Natom
          zatom(i)=nint( Zps(ki_atom(i)) )
       end do
    end if

    if ( .not.allocated(id) ) then
       allocate( id(0:nprocs-1) )
       allocate( ir(0:nprocs-1) )
    end if

    if ( .not.allocated(ipair) ) then
       mpair = ( Natom*(Natom+1) )/2
       call prep_parallel(mpair,nprocs,ir,id)
       mpair = ir(myrank)
       allocate( ipair(2,mpair) )
       m=0
       n=0
       do j=1,Natom
       do i=j,Natom
          m=m+1
          if ( m < id(myrank)+1 .or. id(myrank)+ir(myrank) < m ) cycle
          n=n+1
          ipair(1,n)=i
          ipair(2,n)=j
       end do
       end do
    end if

    Qtot =0.d0
    Qtot2=0.d0
    do i=1,Natom
       Qtot =Qtot +zatom(i)
       Qtot2=Qtot2+zatom(i)**2
    end do

    const1 = 4.d0*Pi/Vcell

    if ( eta == 0.d0 ) eta = pi/Vcell**(2.d0/3.d0)

!--- G

    const=1.d0/(4.d0*eta)
    if ( g_0 == 0.d0 ) g_0=12.d0*sqrt(eta)
    g=g_0

    do loop_g=1,max_loop_g

       ecut=g*g

       call search_grid(bb,ecut,m1,m2,m3,mg_tot)
       call prep_parallel(mg_tot,nprocs,ir,id)
       mg=ir(myrank)
       if ( allocated(LG) ) deallocate( LG )
       allocate( LG(3,mg) )
       call construct_grid(bb,ecut,mg_tot,m1,m2,m3,mg,LG)

       call watch(ct0,et0)
       call calc_ewldg(mg,LG,ewldg)
       ewldg=ewldg*4.d0*pi/Vcell
       ewldg=ewldg-pi/(eta*Vcell)*Qtot**2
       ewldg=ewldg-2.d0*sqrt(eta/pi)*Qtot2
       call watch(ct1,et1) ; timg(1)=ct1-ct0 ; timg(2)=et1-et0

       ir(myrank)=mg
       call mpi_allgather &
            (mg,1,mpi_integer,ir,1,mpi_integer,mpi_comm_world,ierr)

!       if ( disp_sw ) then
!          write(*,'(1x,i4,2i8,2g26.15,2g12.5)') &
!               loop_g,mg_tot,mg,ewldg,sqrt(ecut),timg(1:2)
!       end if

       if ( abs((ewldg-ewldg_0)/ewldg) < 1.d-15 ) then
          ewldg=ewldg_0
          exit
       else
          ewldg_0=ewldg
          g_0=g
          g=(g+max_bb)*factor
       end if

    end do ! loop_g

!    if ( disp_sw ) then
!       write(*,'(1x,"ewald(G)  =",2g24.15,2i8)') ewldg*0.5d0,sqrt(ecut),mg,mg_tot
!    end if

!--- R

    alpha = sqrt(eta)
    if ( r_0 == 0.d0 ) r_0=6.d0/sqrt(eta)
    r=r_0

    ewldr     = ewldg ; ewldg=0.d0
    ewldr_tmp = ewldr
    mr_tmp    = 0
    timr(:)   = 0.d0
    do loop_r=1,max_loop_r

       rrcut=r*r

       call watch(ct0,et0)

       call search_grid(aa,rrcut,m1,m2,m3,mr_tot)
       mr=mr_tot ; id(:)=0 ; ir(:)=mr

       allocate( grd_chk(-m1:m1,-m2:m2,-m3:m3) ) ; grd_chk=0

       if ( loop_r == 1 ) then
       else
          allocate( grd_chk0(-m1:m1,-m2:m2,-m3:m3) ) ; grd_chk0=0
          mr_tmp = size(LR,2)
          do i=1,mr_tmp
             grd_chk0( LR(1,i),LR(2,i),LR(3,i) ) = 1
          end do
          m=size(grd_chk)
          call mpi_allreduce(grd_chk0,grd_chk,m &
               ,mpi_integer,mpi_sum,mpi_comm_world,ierr)
          deallocate( grd_chk0 )
       end if

       if ( allocated(LR) ) deallocate(LR)
       allocate( LR(3,mr) )

       call construct_grid(aa,rrcut,mr_tot,m1,m2,m3,mr,LR)

       if ( allocated(LR_tmp) ) deallocate( LR_tmp )
       allocate( LR_tmp(3,mr) ) ; LR_tmp=0

       if ( loop_r == 1 ) then
          m=mr
          LR_tmp(1:3,1:mr) = LR(1:3,1:mr)
       else
          m=0
          do i=1,mr
             if ( grd_chk(LR(1,i),LR(2,i),LR(3,i)) == 0 ) then
                m=m+1
                LR_tmp(1:3,m) = LR(1:3,i)
             end if
          end do
       end if

       deallocate( grd_chk )

       call watch(ct1,et1)

       call watch(ct0,et0)
       call calc_ewldr(m,LR_tmp,ewldr)
       call watch(ct1,et1)
       timr(1)=timr(1)+ct1-ct0 ; timr(2)=timr(2)+et1-et0

       deallocate( LR_tmp )

!       if ( disp_sw ) then
!          write(*,'(1x,i4,3i8,2g26.15,2g12.5)') &
!               loop_r,mr_tot,mr,m,ewldr,sqrt(rrcut),timr(1:2)
!       end if

       if ( abs((ewldr-ewldr_0)/ewldr) < 1.d-15 ) then
          ewldr=ewldr_0
          exit
       else
          ewldr_0=ewldr
          r_0=r
          r=(r+max_aa)*factor
       end if

    end do ! loop_r

!    if ( disp_sw ) then
!       write(*,'(1x,"ewald(G+R)=",2g24.15,2i8)') ewldr*0.5d0,sqrt(rrcut),mr,mr_tot
!    end if

    Ewld=0.5d0*(ewldr+ewldg)

    call write_border( 0, " calc_ewald(end)" )

    return
  END SUBROUTINE calc_ewald

END MODULE ewald_module

max_stars_repo_name: j-iwata/RSDFT

max_issues_repo_name: j-iwata/RSDFT

max_forks_repo_name: j-iwata/RSDFT