text:
!     SPH: INTEGER(iprec) -> INTEGER
      MODULE bfield_grid_mod
!--------------------------------------------------------------------------
!
! FUNCTION: Module containing routines necessary to compute the line-integrated
!           Faraday rotation angle and the line-integrated polarization phase shift.
!
! modified by JS 11/8/06 to use VMEC instead of AJAX for the flux coord transformations
!------------------------------------------------------------------------------

      USE stel_kinds
      USE stel_constants
      USE ajax_mod
      USE safe_open_mod
      USE v3_utilities     ! contains subroutine assert_eq
      USE v3read_wout      ! reads equilibrium from a wout file into AJAX (includes READEQ_FROM_WOUT)
      USE math_utilities   ! module containing vector manipulation routines
      USE b_transform_mod  ! module containing B field transformation routines
      USE b_transform_vmec  ! module containing VMEC B field transformation routines
      IMPLICIT NONE

      REAL(rprec), ALLOCATABLE, DIMENSION(:,:,:) :: flux_surface     ! array of pts along each flux surface

      INTEGER, PARAMETER :: USE_VMEC = 0
      INTEGER, PARAMETER :: USE_FS = 1

      INTEGER :: transform_select = USE_VMEC
!      INTEGER :: transform_select = USE_FS

      CONTAINS
!===================================================================================
      SUBROUTINE B_GRID_DRIVER()
!--------------------------------------------------------------------------------
!
! Function:  Top-level "driver" subroutine for module "bfield_grid_mod"
!----------------------------------------------------------------------------------
      USE faraday_mod   ! contains conversion routine TRACK_TRANS2FLX
      USE b_transform_vmec, ONLY: VMEC_B_INIT, VMEC_CYL2FLX
      USE safe_open_mod
      IMPLICIT NONE

      INTEGER :: i, j, iflag, ierr_vmec
      INTEGER :: istat
      INTEGER :: m, n
      INTEGER :: ngridpts_R      ! number of grid POINTS along R
      INTEGER :: ngridpts_Z      ! number of grid POINTS along Z
      INTEGER :: m_R             ! R dimension for Bcar array
      INTEGER :: n_Z             ! Z dimension for Bcar array
      INTEGER       :: ns_rho        ! number of radial grid pts in VMEC arrays
!      REAL(rprec)       :: rho
      REAL(rprec)       :: stepsize_R
      REAL(rprec)       :: stepsize_Z
      REAL(rprec)       :: Rmin_temp, Rmax_temp
      REAL(rprec)       :: Zmin_temp, Zmax_temp

      REAL(rprec), DIMENSION(3) :: dims
      REAL(rprec), DIMENSION(3) :: cyl_vec
      REAL(rprec), DIMENSION(3) :: x_cyl_temp
      REAL(rprec), DIMENSION(3) :: xflux_ajax         ! position vector in AJAX definition of flux coords
      REAL(rprec), DIMENSION(3) :: xflux_vmec         ! position vector from VMEC calc of flux coords
      REAL(rprec), DIMENSION(3) :: xflux, xflux2
      REAL(rprec), DIMENSION(6) :: g_cyl
      REAL(rprec), ALLOCATABLE, DIMENSION(:,:,:) :: B_car
      REAL(rprec), ALLOCATABLE, DIMENSION(:) :: rho_rz      ! VMEC radial flux coord array
      REAL(rprec), DIMENSION(3) :: B_car_temp
      REAL(rprec), DIMENSION(3) :: B_temp, B_temp2
      REAL(rprec) :: testarray(5)
      CHARACTER(len=80)  :: nml_infile
      CHARACTER(len=80)  :: default_nml_infile
      CHARACTER(len=80)  :: message

      CHARACTER(len=80) :: mychar1
      CHARACTER(len=*), PARAMETER :: subname = 'B_GRID_DRIVER: '
      LOGICAL        :: QUERY_USER_FOR_NML_FILE
      INTEGER :: iou_wout

!............NML variables...................................................!
      INTEGER :: nsegments_R       ! number of grid SEGMENTS along R
      INTEGER :: nsegments_Z       ! number of grid SEGMENTS along Z
      INTEGER :: nsegments_theta   ! number of grid SEGMENTS along poloidal flux coord theta
      REAL(rprec)    :: Rmin, Rmax
      REAL(rprec)    :: Zmin, Zmax
      REAL(rprec)    :: phi
      CHARACTER(len=80) :: wout_infile
      CHARACTER(len=80) :: cdf_outfile

      NAMELIST /mylist/ Rmin, Rmax, Zmin, Zmax,                                &
     &                  nsegments_R, nsegments_Z,  phi,                        &
     &                  nsegments_theta, wout_infile, cdf_outfile

      write(*,*) '===================================================='
      write(*,*) ' '
!      write(*,*) 'Hi.  Now in ', subname

!      QUERY_USER_FOR_NML_FILE = .false.
      QUERY_USER_FOR_NML_FILE = .true.

!..........Create a name for a default nml file that can be used in place of user input........!
      default_nml_infile = 'bgrid_in.nml'

      if ( QUERY_USER_FOR_NML_FILE) then
        message = "Enter nml filename or select 'd' for default file"
        write(*,*) message
        READ(*,*) nml_infile

        if( nml_infile == 'd') then
          write(*,*) 'Using default nml file: ', default_nml_infile
          nml_infile = default_nml_infile
        end if
      else
!......."hard-wire" in the name of the default input nml file containing the other values.........!
        nml_infile = default_nml_infile
      end if

!..........open the input nml file..........................!
      CALL safe_open(iou_wout, istat, TRIM(nml_infile),                 &
     &               'old', 'formatted')
      CALL assert_eq(0,istat,subname //                                       &
     &   ' Safe_open of input NML file failed')

      READ( UNIT = iou_wout, NML = mylist)
      CLOSE(iou_wout)

      write(*,*) subname, 'Rmin = ', Rmin
      write(*,*) subname, 'Rmax = ', Rmax
      write(*,*) subname, 'Zmin = ', Zmin
      write(*,*) subname, 'Zmax = ', Zmax
      write(*,*) subname, 'm = ', nsegments_R
      write(*,*) subname, 'n = ', nsegments_Z
      write(*,*) subname, 'phi = ', phi
      write(*,*) subname, 'nsegments_theta = ', nsegments_theta
      write(*,*) subname, 'input wout = ', wout_infile
      write(*,*) subname, 'cdf_outfile = ', cdf_outfile

!...........fill AJAX stuff in memory from the input wout file.................!
      call DEAL_WITH_WOUT(wout_infile, ns_rho)

!........At this point , AJAX stuff for the wout should be in memory................!
!......also note that wout variables are now globally available via a USE call to read_wout_mod...!

!........initialize the VMEC flux coord transformation routines..........!
      call VMEC_B_INIT

!.......fill the array of (R,Z) points that describe each flux surface...........!
      call FILL_FLUX_SURFACE_ARRAY(ns_rho, phi, nsegments_theta)

!!.......write out the flux_surface array to std output to verify that it filled correctly....!
!      do i = 1, ns_rho
!        write(*,*) 'flux suface number = ', i
!        do j = 1, nsegments_theta
!          write(*,*) flux_surface(i,j,:)
!         end do
!      end do

!.......convert the # of grid SEGMENTS into the # of grid POINTS (note that it is.....
!.......the grid point values (m_R, n_Z) that define the array size)................!
      ngridpts_R = nsegments_R + 1
      ngridpts_Z = nsegments_Z + 1

!.......allocate the B_car array.......................................!
      m_R = nsegments_R
      n_Z = nsegments_Z
      ALLOCATE( B_car(0:m_R, 0:n_Z, 3) )

!!.....DEBUGGING: Hard-wire cyl_vec values that give results inside the plasma. JS 12/16/05
!    (1.602,2.114, -0.117)   and (1.606, 2.104, -0.0738)
!      cyl_vec(1) = 1.602
!      cyl_vec(2) = 2.114
!      cyl_vec(3) = -0.117
!      Rmin_temp = 1.45
!      Rmax_temp = 1.606
!      Zmin_temp = -0.117
!      Zmax_temp = -0.0738

      stepsize_R = ABS(Rmax - Rmin) / REAL(nsegments_R)
      stepsize_Z = ABS(Zmax - Zmin) / REAL(nsegments_Z)

      do m = 0, m_R
        do n = 0, n_Z

!........Note that the x_cyl order in AJAX is the conventional (r,phi,z)...........!
          x_cyl_temp(1) = Rmin + m * stepsize_R
          x_cyl_temp(2) = phi
          x_cyl_temp(3) = Zmin + n * stepsize_Z
          call TRACK_TRANS2FLX(x_cyl_temp, xflux_ajax, iflag, message)

          if (iflag .ne. 0 ) then
            write(*,*) 'ERROR in TRACK_TRANS2FLX! MESSAGE IS: ', message
          end if

!...............use VMEC to compute flux coord position vector for debugging.  JS 12/4/06
          call VMEC_CYL2FLX(x_cyl_temp, xflux_vmec, ierr_vmec)

!           write(*,*) 'm = ', m, 'n = ', n
!           write(*,*) subname, 'x_cyl_temp = ', x_cyl_temp
!           write(*,*) subname, 'xflux_ajax = ', xflux_ajax

!!................verify that the current grid point is INSIDE the plasma................!
!!.........(and set upper limit to 0.9801 instead of 1.0 to avoid weird plasma "edge effects")......!
!          if ( xflux_ajax(1) .gt. 0.9801 ) then

!...............replace AJAX plasma in/out check with FS and/or VMEC check.  JS 12/2/06....!
!          if ( INVALID_RHO(x_cyl_temp) ) then
          if ( INVALID_RHO_VMEC(x_cyl_temp) ) then
            write(*,*) subname, 'WARNING: POINT IS OUTSIDE THE PLASMA!',       &
     &                 'xflux_ajax(1) = ',  xflux_ajax(1)

            B_car(m,n,:) = 0.0
          else
!            write(*,*) subname, 'xflux_ajax(1) = ', xflux_ajax(1)

!...............find B field at point x_cyl_temp & save to B_car array.............!
!             write(*,*) subname, 'now calling COMPUTE_B_CAR'
              call COMPUTE_B_CAR( x_cyl_temp, B_temp)
              B_car(m,n,:) = B_temp
          end if
        end do
      end do

!!..........write B_car arrays to std output for debugging (grouping by m [R] values)...............!
!      do i = 0,m_R
!          write(*,*) '----------------------------------------------'
!          write(*,*) subname, 'm value = ', i
!        do j = 0, n_Z
!          write(*,*) subname, 'B_car = ', B_car(i,j,:)
!        end do
!      end do
!
!!..........write B_car arrays to std output for debugging.(grouping by n [Z] values)...............!
!      do j = 0, n_Z
!          write(*,*) '----------------------------------------------'
!          write(*,*) subname, 'n value = ', j
!        do i = 0,m_R
!          write(*,*) subname, 'B_car = ', B_car(i,j,:)
!        end do
!      end do

!......write the results to an output netCDF file.........................!
      call OUTPUT_TO_CDF(Rmin, Rmax, Zmin, Zmax, nsegments_R,                       &
     &                   nsegments_Z, phi, B_car, ns_rho,                           &
     &                   nsegments_theta, cdf_outfile)

!........now read file back in to verify that the data was output to the .nc correctly....!
!      call READ_TEST_FOR_CDF(cdf_outfile)

      RETURN
      END SUBROUTINE B_GRID_DRIVER

!===================================================================================
      SUBROUTINE DEAL_WITH_WOUT(wout_infile, nr_rz)
!--------------------------------------------------------------------------------
!
! Function:  Retrieve necessary info from VMEC wout file
!
!----------------------------------------------------------------------------------
      USE read_wout_mod  ! module that holds VMEC variables in memory
      IMPLICIT NONE

!............dummy variables......................................!
      CHARACTER(len=*), INTENT(IN)  :: wout_infile
      INTEGER, INTENT(OUT) :: nr_rz ! number of radial grid pts in VMEC arrays

!............local variables......................................!
      INTEGER :: i, j, iflag
      INTEGER :: istat
      REAL(rprec), DIMENSION(3) :: dims
      REAL(rprec), DIMENSION(3) :: xflux
      REAL(rprec), DIMENSION(3) :: x_cyl
      CHARACTER(len=80)  :: message
      CHARACTER(len=*), PARAMETER :: subname = 'DEAL_WITH_WOUT: '

!      write(*,*) '===================================================='
!      write(*,*) ' '
!      write(*,*) 'Hi.  Now in ', subname

!...........fill AJAX stuff in memory from the input wout file.................!
      call READEQ_FROM_WOUT(wout_infile)

!........At this point , AJAX stuff for the wout should be in memory.  Also note that wout variables..........!
!......are now globally available to this subroutine due the call to USE read_wout_mod ...!

!.......extract the number of radial grid points from global VMEC variables.............!
      dims(1:2) = shape(rmnc)
      nr_rz = dims(2)

      write(*,*) subname, 'nr_rz = ', nr_rz, 'ns = ', ns

      if ( nr_rz .ne. ns ) then
        message = 'WARNING! MISMATCH IN # OF VMEC RADIAL GRIDPTS!!'
        write(*,*) subname, message
      end if

      RETURN
      END SUBROUTINE DEAL_WITH_WOUT

!===================================================================================
      SUBROUTINE FILL_FLUX_SURFACE_ARRAY(ns_rho, input_phi,                     &
     &                                   nsegments_theta)
!--------------------------------------------------------------------------------
!
! Function:  Fill an array with cylinderical coordinate points along the flux surfaces
!
!----------------------------------------------------------------------------------
      USE ajax_mod          ! contains conversion routine AJAX_FLX2CYL
      USE b_transform_vmec, ONLY: VMEC_FLX2CYL_JS
      IMPLICIT NONE

!...........dummy variables...............................................!
      REAL(rprec), INTENT(IN)       :: input_phi           ! phi value specified in input NML file
      INTEGER, INTENT(IN)    :: ns_rho              ! number of radial grid pts in VMEC arrays
      INTEGER, INTENT(IN)    :: nsegments_theta     ! # of grid segments along the theta path

!...........local variables...............................................!
      INTEGER    :: i, j, iflag
      INTEGER    :: istat
      REAL(rprec)       :: theta        ! poloidal flux coordinate
      REAL(rprec)       :: delta_theta  ! grid step size in the poloidal direction

      REAL(rprec), DIMENSION(3) :: xflux
      REAL(rprec), DIMENSION(3) :: x_cyl
      REAL(rprec), DIMENSION(3) :: x_cyl_temp, xflx_temp
      REAL(rprec), DIMENSION(3) :: B_con_temp
      REAL(rprec), ALLOCATABLE, DIMENSION(:) :: rho_rz      ! VMEC radial flux coord array
      CHARACTER(len=80)  :: message

      CHARACTER(len=80) :: mychar1
      CHARACTER(len=*), PARAMETER :: subname='FILL_FLUX_SURFACE_ARRAY: '
      INTEGER :: iou_wout

      write(*,*) '===================================================='
      write(*,*) ' '
!      write(*,*) 'Hi.  Now in ', subname

      ALLOCATE( rho_rz(1:ns_rho) )

!...........allocate output flux surface array........................!
      ALLOCATE( flux_surface(ns_rho, nsegments_theta, 2) )

!......fill rho_rz as an evenly spaced "full-grid", following Houlberg.........!
      do i = 1, ns_rho
        rho_rz(i) = SQRT( REAL(i-1) / REAL(ns_rho - 1) )
      end do

      if ( nsegments_theta .ne. 0 ) then
      delta_theta = (2.0 * PI) / REAL(nsegments_theta)
      else
        delta_theta = 0.0
      end if

!..........convert each point along the flux surface contour into cylinderical coords..!

      do i = 1, ns_rho
        do j = 1, nsegments_theta

          xflux(1) = rho_rz(i)
          xflux(2) = (j-1) * delta_theta
          xflux(3) = input_phi

!!...........transform flux vector xflux() into cylin coord vector x_cyl() via AJAX
!          call AJAX_FLX2CYL(xflux, x_cyl, iflag, message)

!...........transform flux vector xflux() into cylin coord vector x_cyl() via VMEC
          call VMEC_FLX2CYL_JS(xflux, x_cyl, iflag)

!...........fill the array of (R,Z) points along each flux surface (note: R=1, Z=2)................!
          flux_surface(i,j,1) = x_cyl(1)
          flux_surface(i,j,2) = x_cyl(3)
         end do
      end do

      RETURN
      END SUBROUTINE FILL_FLUX_SURFACE_ARRAY

!===================================================================================
      SUBROUTINE OUTPUT_TO_CDF(Rmin, Rmax, Zmin, Zmax, nsegments_R,               &
     &                         nsegments_Z, phi, B_Car, ns_rho,                   &
     &                         nsegments_theta, cdf_outfile)
!--------------------------------------------------------------------------------
!
! Function:  Writes filled variables to an output netCDF file.
!
! Global inputs: flux_surface
! Global outputs: none
!
! CAVEAT:   There is no mechanism analogous to "safe_open" for the cdf_open subroutine
!         (at least none that I'm aware of...).  Therefore, it is *theoretically* possible
!         that the i/o unit number for the file could conflict with a previously opened
!          i/o unit.
!----------------------------------------------------------------------------------
      USE ezcdf
      IMPLICIT NONE

!.........passed dummy variables...............................!
      INTEGER :: nsegments_R
      INTEGER :: nsegments_Z
      INTEGER :: ns_rho
      INTEGER :: nsegments_theta
      REAL(rprec) :: Rmin, Rmax
      REAL(rprec) :: Zmin, Zmax
      REAL(rprec) :: phi
      REAL(rprec), DIMENSION(:,:,:) :: B_car
      CHARACTER(len=*) :: cdf_outfile

!...........local variables..............................!
      INTEGER           :: counter = 0
      INTEGER           :: i, j, nbeams, n_edge, iflag
      INTEGER           :: n_ip_units, istat
      REAL(rprec) :: testarray(5)

      CHARACTER(len=80)  :: nml_infile

      CHARACTER(len=80) :: mychar1
      CHARACTER(len=80) :: mychar2, message
      CHARACTER(len=*), PARAMETER :: subname = 'OUTPUT_TO_CDF: '

      INTEGER :: iou_cdf = 132

!      write(*,*) 'Hi.  Now in ', subname

!........enter values for testarray.....!
      testarray(1) = 5.0
      testarray(2) = 6.0
      testarray(3) = 7.0
      testarray(4) = 8.0
      testarray(5) = 9.0

!      write(*,*) subname, 'testarray = ', testarray

!..........open the output netCDF file..........................!

      call cdf_open(iou_cdf, cdf_outfile, 'w', istat)
      CALL cdf_define(iou_cdf, 'Rmin', Rmin)
      CALL cdf_define(iou_cdf, 'Rmax', Rmax)
      CALL cdf_define(iou_cdf, 'Zmin', Zmin)
      CALL cdf_define(iou_cdf, 'Zmax', Zmax)
      CALL cdf_define(iou_cdf, 'phi', phi)
      CALL cdf_define(iou_cdf, 'nsegments_R', nsegments_R)
      CALL cdf_define(iou_cdf, 'nsegments_Z', nsegments_Z)
      CALL cdf_define(iou_cdf, 'ns_rho', ns_rho)
      CALL cdf_define(iou_cdf, 'nsegments_theta', nsegments_theta)
!      CALL cdf_define(iou_cdf, 'testarray', testarray)
      CALL cdf_define(iou_cdf, 'B_car', B_car)
      CALL cdf_define(iou_cdf, 'flux_surface', flux_surface)

      CALL cdf_write(iou_cdf, 'Rmin', Rmin)
      CALL cdf_write(iou_cdf, 'Rmax', Rmax)
      CALL cdf_write(iou_cdf, 'Zmin', Zmin)
      CALL cdf_write(iou_cdf, 'Zmax', Zmax)
      CALL cdf_write(iou_cdf, 'phi', phi)
      CALL cdf_write(iou_cdf, 'nsegments_R', nsegments_R)
      CALL cdf_write(iou_cdf, 'nsegments_Z', nsegments_Z)
      CALL cdf_write(iou_cdf, 'ns_rho', ns_rho)
      CALL cdf_write(iou_cdf, 'nsegments_theta', nsegments_theta)
!      CALL cdf_write(iou_cdf, 'testarray', testarray)
      CALL cdf_write(iou_cdf, 'B_car', B_car)
      CALL cdf_write(iou_cdf, 'flux_surface', flux_surface)
      call cdf_close(iou_cdf)

!      write(*,*) subname, 'Rmin = ', Rmin
!      write(*,*) subname, 'Rmax = ', Rmax
!      write(*,*) subname, 'Zmin = ', Zmin
!      write(*,*) subname, 'Zmax = ', Zmax
!      write(*,*) subname, 'm = ', nsegments_R
!      write(*,*) subname, 'n = ', nsegments_Z
!      write(*,*) subname, 'phi = ', phi
!      write(*,*) subname, 'cdf_outfile = ', cdf_outfile

      RETURN
      END SUBROUTINE OUTPUT_TO_CDF

!===================================================================================
      SUBROUTINE READ_TEST_FOR_CDF(cdf_infile)
!--------------------------------------------------------------------------------
!
! Function:  Test routine to verify that data was written to netCDF correctly
!
! CAVEAT:   There is no mechanism analogous to "safe_open" for the cdf_open subroutine
!         (at least none that I'm aware of...).  Therefore, it is *theoretically* possible
!         that the i/o unit number for the file could conflict with a previously opened
!          i/o unit.
!----------------------------------------------------------------------------------
      USE ezcdf
      IMPLICIT NONE

!.........passed dummy variables...............................!

      CHARACTER(len=*) :: cdf_infile

!...........local variables..............................!
      INTEGER           :: counter = 0
      INTEGER           :: i, j, nbeams, n_edge, iflag
      INTEGER           :: n_ip_units, istat

      INTEGER :: nsegments_R
      INTEGER :: nsegments_Z
      INTEGER :: ns_rho
      INTEGER :: nsegments_theta
      REAL(rprec) :: Rmin, Rmax
      REAL(rprec) :: Zmin, Zmax
      REAL(rprec) :: phi
      REAL(rprec), ALLOCATABLE, DIMENSION(:) :: testarray
      REAL(rprec), ALLOCATABLE, DIMENSION(:,:,:) :: B_car
      REAL(rprec), ALLOCATABLE, DIMENSION(:,:,:) :: my_fsurface

      CHARACTER(len=80) :: mychar1
      CHARACTER(len=80) :: mychar2, message
      CHARACTER(len=*), PARAMETER :: subname = 'READ_TEST_FOR_CDF: '

      INTEGER :: iou_cdf = 132
      INTEGER ::  arr_size
      INTEGER ::  size1, size2, size3

      INTEGER, DIMENSION(3) :: dim

      write(*,*) 'Hi.  Now in ', subname

!........open the input .nc file and read in the variables..............!
      call cdf_open(iou_cdf, cdf_infile, 'r', istat)

!      call cdf_inquire(iou_cdf, 'testarray', dim)
      call cdf_inquire(iou_cdf, 'B_car', dim)

!      write(*,*) subname, 'dim for testarray = ', dim
      write(*,*) subname, 'dim for B_car = ', dim

!!........allocate testarray size based on cdf_inquire results.......!
!      arr_size = dim(1)
!      write(*,*) subname, 'dim for testarray = ', arr_size
!      ALLOCATE( testarray(arr_size) )
!
!!........read in values for testarray...........................!
!      call cdf_read(iou_cdf, 'testarray', testarray)

!!........allocate B_car size based on cdf_inquire results.......!
      size1 = dim(1)
      size2 = dim(2)
      size3 = dim(3)
      write(*,*) subname, 'dim for B_car = ', dim
      ALLOCATE( B_car(0:(size1 -1),0:(size2 -1),size3) )

!........read in values for B_car...........................!
      call cdf_read(iou_cdf, 'B_car', B_car)

!.......repeat allocation process for flux_surface array...............!
      call cdf_inquire(iou_cdf, 'flux_surface', dim)
       size1 = dim(1)
       size2 = dim(2)
       size3 = dim(3)
      write(*,*) subname, 'dim for B_car = ', dim
      ALLOCATE( my_fsurface(size1, size2, size3) )

!........read in flux_surface values into local array "my_fsurface"...........................!
      call cdf_read(iou_cdf, 'flux_surface', my_fsurface)

!..........read in non-array variables.....................!
      CALL cdf_read(iou_cdf, 'Rmin', Rmin)
      CALL cdf_read(iou_cdf, 'Rmax', Rmax)
      CALL cdf_read(iou_cdf, 'Zmin', Zmin)
      CALL cdf_read(iou_cdf, 'Zmax', Zmax)
      CALL cdf_read(iou_cdf, 'phi', phi)
      CALL cdf_read(iou_cdf, 'nsegments_R', nsegments_R)
      CALL cdf_read(iou_cdf, 'nsegments_Z', nsegments_Z)
      CALL cdf_read(iou_cdf, 'ns_rho', ns_rho)
      CALL cdf_read(iou_cdf, 'nsegments_theta', nsegments_theta)
      call cdf_close(iou_cdf)

      write(*,*) subname, 'Rmin = ', Rmin
      write(*,*) subname, 'Rmax = ', Rmax
      write(*,*) subname, 'Zmin = ', Zmin
      write(*,*) subname, 'Zmax = ', Zmax
      write(*,*) subname, 'm = ', nsegments_R
      write(*,*) subname, 'n = ', nsegments_Z
      write(*,*) subname, 'phi = ', phi
      write(*,*) subname, 'ns_rho = ', ns_rho
      write(*,*) subname, 'nsegments_theta = ', nsegments_theta

!!.........test read-in of B_car array............!
!      do i = 0,nsegments_R
!        do j = 0, nsegments_Z
!          write(*,*) subname, 'B_car = ', B_car(i,j,:)
!        end do
!      end do

!.........test read-in of flux_surface array............!
!      do i = 1, ns_rho
      do i = 30, ns_rho
        write(*,*) subname, 'flux suface number = ', i
        do j = 1, nsegments_theta
          write(*,*) subname, my_fsurface(i,j,:)
         end do
      end do

      RETURN
      END SUBROUTINE READ_TEST_FOR_CDF

      end module bfield_grid_mod


text:
!
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!
!           c_sizeof() can be applied to elements of a pointer, allocatable,
!           assumed-shape, and assumed-size array
module bi04
   use iso_c_binding
   integer result(4)
   integer expect(4)
   data expect/4,8,4,8/
   contains
   subroutine sub(aa, pp, yy, zz)
     real, allocatable :: aa(:)
     real, pointer :: pp(:)
     real(8)       :: yy(:)
     real(8)       :: zz(*)
     result(1) = c_sizeof(pp(1))
     result(2) = c_sizeof(yy(1))
     result(3) = c_sizeof(aa(1))
     result(4) = c_sizeof(zz(1))
!     print 99, result
!     print 99, expect
!  99 format(6i4)
     call check(result, expect, 4)
    endsubroutine
end
use bi04
real, allocatable :: aa(:)
real, pointer :: pp(:)
real(8)       :: yy(7)
allocate(aa(0:4))
allocate(pp(3:9))
call sub(aa, pp, yy, yy)
end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
Angelique Donovan, aka Angel Donovan is part of the Management Team for Academy Lane Apartments.

Hello, I am the Property Manager at Academy Lane Apartments, previously Americana Arms.  I was moved up from the Los Angeles area with my husband Jason, by Virtu Property Management when they purchased the property June 2012. My husband is also part of the Management Team, as the Head of Maintenance. Since we took over the property, we have been working on improvements  not only what you see from the street, but also structural repairs, insulation in all the interior and exterior walls, and just recently all the new paint/trim and carpet in the interior hallways.

Come check out all our hard work!

Besides managing the property, we also have a total of 4 kids. His daughters are 15, one of which lives with us and the other lives and goes to school in LA with her mother. We also have 2 kids together, a son and a daughter, 8 and 11 who attend North Davis Elementary... We love this town, and have enjoyed the changes so far.

20130319 16:11:20 nbsp Welcome to the Wiki/Business Owner Welcome to the Wiki! Nice page and pics. Users/PeteB

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*----------------------------------------------------------------------*
      subroutine store_def_intm(fl_item,
     &     label,occ,rst,nj,nblk,
     &     parent1,parent2,tra_new,tra1,tra2,
     &     orb_info)
*----------------------------------------------------------------------*
*     write info for creating new intermediate on formula list item
*     (fl_item)
*     the operator labels parent1 and parent2 are needed for
*     obtaining the list-attributes when evaluating the formula
*     (e.g. the IRREP of the new intermediate is the direct product
*      of the IRREPs of parent1 and parent2; parent2 may be empty)
*----------------------------------------------------------------------*

      implicit none

      include 'stdunit.h'
      include 'opdim.h'
      include 'def_contraction.h'
      include 'def_operator.h'
      include 'def_formula_item.h'
      include 'def_orbinf.h'

      integer, parameter ::
     &     ntest = 00

      type(formula_item), intent(inout) ::
     &     fl_item
      character(len=*), intent(in) ::
     &     label, parent1, parent2
      type(orbinf), intent(in) ::
     &     orb_info
      integer, intent(in) ::
     &     nj, nblk,
     &     occ(ngastp,2,nj*nblk),
     &     rst(2,orb_info%ngas,2,2,orb_info%nspin,nj*nblk)
      logical, intent(in) ::
     &     tra_new, tra1, tra2

      if (ntest.ge.100) write(lulog,*) 'storing def of intermediate'

      fl_item%interm%n_occ_cls = nblk
      fl_item%interm%njoined = nj
      fl_item%interm%name = trim(label)
      call init_operator(fl_item%interm,orb_info)

      fl_item%interm%ihpvca_occ = occ
      fl_item%interm%igasca_restr = rst

      call occ2caocc(fl_item%interm%ica_occ,occ,nj,nblk)

      fl_item%interm%formal_blk = .false.

      fl_item%parent1 = trim(parent1)
      fl_item%parent2 = trim(parent2)

      fl_item%tra = tra_new
      fl_item%tra1 = tra1
      fl_item%tra2 = tra2

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
!> Backward auto diff
program main

    use auto_diff, only: tree_t, operator(*), operator(+), exp, rk, assignment(=)
    implicit none
    type(tree_t) :: a, b, c
    type(tree_t) :: y

    a = 2.0_rk
    b = 1.0_rk
    c = 0.0_rk

    print *, "demo1: y = (a + b*b)*exp(c)"

    y = (a + b*b)*exp(c)

    print *, "y     = ", y%get_value()
    call y%backward()

    print *, "dy/da = ", a%get_grad()
    print *, "dy/db = ", b%get_grad()
    print *, "dy/dc = ", c%get_grad()

    ! - - -

    a = 2.0_rk
    b = 1.0_rk

    print *, ""
    print *, "demo2: y = (a + b)*(b + 1.0)"

    y = (a + b)*(b + 1.0_rk)

    print *, "y     = ", y%get_value()

    call y%backward()

    print *, "dy/da = ", a%get_grad()
    print *, "dy/db = ", b%get_grad()

end program main

!>  demo1: y = (a + b*b)*exp(c)
!>  y     =    3.0000000000000000
!>  dy/da =    1.0000000000000000
!>  dy/db =    2.0000000000000000
!>  dy/dc =    3.0000000000000000
!>
!>  demo2: y = (a + b)*(b + 1.0)
!>  y     =    6.0000000000000000
!>  dy/da =    2.0000000000000000
!>  dy/db =    5.0000000000000000

max_stars_repo_name: zoziha/Auto-Diff

max_issues_repo_name: zoziha/Auto-Diff

max_forks_repo_name: zoziha/Auto-Diff


text:
      PROGRAM xchebft
C     driver for routine chebft
      INTEGER NVAL
      REAL PIO2,EPS
      PARAMETER(NVAL=40, PIO2=1.5707963, EPS=1E-6)
      INTEGER i,j,mval
      REAL a,b,dum,f,func,t0,t1,term,x,y,c(NVAL)
      EXTERNAL func
      a=-PIO2
      b=PIO2
      call chebft(a,b,c,NVAL,func)
C     test result
10    write(*,*) 'How many terms in Chebyshev evaluation?'
      write(*,'(1x,a,i2,a)') 'Enter n between 6 and ',NVAL,
     *     '. Enter n=0 to END.'
      read(*,*) mval
      if ((mval.le.0).or.(mval.gt.NVAL)) goto 20
      write(*,'(1x,t10,a,t19,a,t28,a)') 'X','Actual','Chebyshev fit'
      do 12 i=-8,8,1
        x=i*PIO2/10.0
        y=(x-0.5*(b+a))/(0.5*(b-a))
C     evaluate Chebyshev polynomial without using routine CHEBEV
        t0=1.0
        t1=y
        f=c(2)*t1+c(1)*0.5
        do 11 j=3,mval
          dum=t1
          t1=2.0*y*t1-t0
          t0=dum
          term=c(j)*t1
          f=f+term
11      continue
        write(*,'(1x,3f12.6)') x,func(x),f
12    continue
      goto 10
20    END
      REAL FUNCTION func(x)
      REAL x
      func=(x**2)*(x**2-2.0)*sin(x)
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Examples/xchebft.for


text:
Dr. Kurt Meyers is a former Sport Psychology Postdoctoral Fellow at USC. Dr. Kurt Meyers has worked with individuals, teams, and coaches at the amateur, olympic, and professional levels. Kurt Meyers, Psy.D. is a licensed clinical psychologist with a Doctorate in Psychology and a Master’s Degree in Sport Psychology. License number PSY23156.

For additional medical and psychological services, see our Medical Services page.


text:
PAUL

Im new to the Wiki! Ive been reading for awhile now and will finally edit some pages. :)

Im a UC Davis alumnus currently living/working in the Bay Area Martinez/El Sobrante area.

20080528 18:28:35 nbsp Welcome to the Wiki. You may want to read through some of the Wiki Community articles to understand why there is a vested interest in retaining history here on the wiki. Users/JasonAller

20080528 18:36:32 nbsp The concerns you raise continue to be discussed on the wiki. See Reviews and Anonymous Criticism. Please review the Importance of using your RealName as it goes a long way when advocating for a change in the way the wiki operates. Users/JasonAller

20080528 19:19:07 nbsp Aggie alumnus as well? Users/PhilipNeustrom

20080528 19:56:50 nbsp Hey, Paul. I did revert/change a couple of your entries so I wanted to give you an explanation. It seemed like you removed quite a bit of the negative information. While it is negative, it is also true. I was at the meetings and know some of the folks involved. I was wondering why you chose to delete it? Is it a question about it being true?

Feel free to send me a message back...

andrew Users/AndrewPeake

20080528 22:21:22 nbsp Hey, I agree that in the scheme of life these things are not all that important. However, as phillip said, there is no reason to delete it. Sure it is a little negative, but people have to live with what they did. Check out my page, being an ASUCD senator isnt much, but by having that info out there it provides a history, and lessons on what not to do for future students. I was involved via mediaboard (I was an designee at the time) Users/AndrewPeake

20080528 22:39:13 nbsp it means I went to meetings on behalf of someone who was on the board. I later was a rep. Media Board oversees the aggie and kdvs. Did you work at the paper? Users/AndrewPeake

20080529 01:06:12 nbsp Hey thanks for all the new info about Aggie stuff. My uncle used to live in Martinez... what do you do there? Users/JeremyOgul

20080529 09:04:42 nbsp Hey, I read what you http://daviswiki.org/Users/PhilipNeustrom?actiondiff&at_date1212035730.2655101 wrote on Philips page. I do think its interesting and important to have all the details; this is a wiki, not an abridged dictionary for Davis. Maybe its not a resource you could quote in a paper but its helpful in giving visitors and residents a relevant history and background of Davis. Unfortunately sometimes negative things happen but that doesnt mean we should censor them.  Users/ElleWeber

20080529 14:42:46 nbsp Things on this website have appeared in articles about me over the years. While the Aggie Ed board revolves continuously, the wiki provides a certain amount of institutional memory. For better or for worse info here stay forever  which I think is a great resource so we can all learn from past mistakes or resurrect old ideas that never came to fruition. Id like to take back a few things Ive said over the years because I have since either learned better, or my intent was lost, etc. But live and let live. And let the wiki admins worry about bandwidth usage and storage capacity. Users/RobRoy


text:
      PROGRAM W_CURL
      USE MOD_ZA  ! HYCOM array I/O interface
      IMPLICIT NONE
C
C     WIND ARRAYS.
C
      INTEGER, ALLOCATABLE :: MSK(:,:)
      REAL*4,  ALLOCATABLE :: RSCX(:,:),RSCY(:,:),PANG(:,:)
      REAL*4,  ALLOCATABLE :: TXM(:,:),TYM(:,:),CURL(:,:),QCURL(:,:),
     &                        TXP(:),TYP(:)
C
      REAL*4,  PARAMETER   :: ONEM=1.0
C
      CHARACTER PREAMBL(5)*79
C
C**********
C*
C 1)  CALCULATE THE CURL OF EXISTING HYCOM WIND STRESS FILES.
C     ALSO WORKS ON 10M WIND VELOCITY FILES.
C
C 2)  INPUT
C        ON UNIT 10:    .a/.b FORMAT MODEL TAUEWD  FILE, SEE (3).
C        ON UNIT 11:    .a/.b FORMAT MODEL TAUNWD  FILE, SEE (3).
C        ON UNIT 21:    .a/.b FORMAT regional.grid FILE.
C     OUTPUT:
C        ON UNIT 12:    .a/.b FORMAT MODEL CURL    FILE, SEE (4).
C
C 3)  THE INPUT WIND STRESSES HAVE THEIR COMPONENTS ON EITHER
C      EVERY POINT OF THE MODEL'S 'U' AND 'V' GRIDS RESPECTIVELY,
C      OR BOTH ON EVERY POINT OF THE MODEL'S 'P' GRID.  IF ON THE
C      P-GRID, IT IS FIRST INTERPOLATED TO THE U AND V GRIDS AND
C      THEN THE CURL CALCULATED FROM THE STRESSES ON THESE GRIDS.
C
C 4)  THE CURL IS CALCULATED ON THE STREAMFUNCTION-GRID, AS IS
C      NATURAL FOR STRESSES ON THE U AND V GRIDS, AND INTERPOLATED
C      TO THE P-GRID FOR OUTPUT.  THIS MAY REDUCE ACCURACY.
C      ARRAY SIZE IS 'IDM' BY 'JDM', AND THE DATA IS OUTPUT .a/.b
C      FORMAT TOGETHER WITH EITHER (A) THE MONTH, OR (B) THE DAY THE
C      WIND REPRESENTS AND THE INCREMENT IN DAYS TO THE NEXT WIND
C      RECORD.
C
C 5)  ALAN J. WALLCRAFT,  NRL,  NOVEMBER 2002.
C*
C**********
C
      LOGICAL      PGRID,GLOBAL,ARCTIC,ROTATE
      CHARACTER*80 CLINE
      INTEGER      I,II,IOS,IM1,IP1,J,JJ,JM1,JP1,KREC
      REAL*4       TXMIJ,TYMIJ,COSPANG,SINPANG
      REAL*4       HMINA,HMINB,HMAXA,HMAXB
      REAL*4       XMIN,XMAX
C
C --- MODEL ARRAYS.
C
      CALL XCSPMD  !define idm,jdm
      ALLOCATE(   MSK(IDM,JDM) )
      ALLOCATE(  RSCX(IDM,JDM) )
      ALLOCATE(  RSCY(IDM,JDM) )
      ALLOCATE(  PANG(IDM,JDM) )
      ALLOCATE(   TXM(IDM,JDM) )
      ALLOCATE(   TYM(IDM,JDM) )
      ALLOCATE(  CURL(IDM,JDM) )
      ALLOCATE( QCURL(IDM,JDM) )
      ALLOCATE(   TXP(0:IDM) )
      ALLOCATE(   TYP(0:JDM) )
C
C     GRID INPUT.
C
      CALL ZAIOST
C
      CALL ZHOPNC(21, 'regional.grid.b', 'FORMATTED', 'OLD', 0)
      CALL ZAIOPF('regional.grid.a', 'OLD', 21)
C
      READ(21,*) ! skip idm
      READ(21,*) ! skip jdm
      READ(21,*) ! skip mapflg
      READ(21,'(A)') CLINE  !read plon range
      WRITE(6,'(A)') CLINE
      I = INDEX(CLINE,'=')
      READ(CLINE(I+1:),*) HMINB,HMAXB
      GLOBAL = HMAXB-HMINB .GT. 350.0
      CALL ZAIORD(RSCX,MSK,.FALSE., HMINA,HMAXA, 21)  !temporarily store plon
      READ(21,*) ! skip plat
      CALL ZAIORD(RSCY,MSK,.FALSE., HMINA,HMAXA, 21)  !temporarily store plat
      ARCTIC = RSCX(3*IDM/4+1,JDM-1).EQ.RSCX(IDM/4,JDM) .AND.
     &         RSCY(3*IDM/4+1,JDM-1).EQ.RSCY(IDM/4,JDM)
      READ(21,*) ! skip qlon
      CALL ZAIOSK(21)
      READ(21,*) ! skip qlat
      CALL ZAIOSK(21)
      READ(21,*) ! skip ulon
      CALL ZAIOSK(21)
      READ(21,*) ! skip ulat
      CALL ZAIOSK(21)
      READ(21,*) ! skip vlon
      CALL ZAIOSK(21)
      READ(21,*) ! skip vlat
      CALL ZAIOSK(21)
C
      READ(21,'(A)') CLINE
      WRITE(6,'(A)') CLINE
      I = INDEX(CLINE,'=')
      READ(CLINE(I+1:),*) HMINB,HMAXB
      CALL ZAIORD(PANG,MSK,.FALSE., HMINA,HMAXA, 21)
      IF     (ABS(HMINA-HMINB).GT.ABS(HMINB)*1.E-4 .OR.
     &        ABS(HMAXA-HMAXB).GT.ABS(HMAXB)*1.E-4     ) THEN
        WRITE(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (pang):',
     &    '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &    '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
        CALL ZHFLSH(6)
        STOP
      ENDIF
C
      READ(21,*) ! skip pscx
      CALL ZAIOSK(21)
      READ(21,*) ! skip pscy
      CALL ZAIOSK(21)
C
      READ(21,'(A)') CLINE
      WRITE(6,'(A)') CLINE
      I = INDEX(CLINE,'=')
      READ(CLINE(I+1:),*) HMINB,HMAXB
      CALL ZAIORD(RSCX,MSK,.FALSE., HMINA,HMAXA, 21)
      IF     (ABS(HMINA-HMINB).GT.ABS(HMINB)*1.E-4 .OR.
     &        ABS(HMAXA-HMAXB).GT.ABS(HMAXB)*1.E-4     ) THEN
        WRITE(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (qscx):',
     &    '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &    '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
        CALL ZHFLSH(6)
        STOP
      ENDIF
C
      READ(21,'(A)') CLINE
      WRITE(6,'(A)') CLINE
      I = INDEX(CLINE,'=')
      READ(CLINE(I+1:),*) HMINB,HMAXB
      CALL ZAIORD(RSCY,MSK,.FALSE., HMINA,HMAXA, 21)
      IF     (ABS(HMINA-HMINB).GT.ABS(HMINB)*1.E-4 .OR.
     &        ABS(HMAXA-HMAXB).GT.ABS(HMAXB)*1.E-4     ) THEN
        WRITE(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (qscy):',
     &    '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &    '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
        CALL ZHFLSH(6)
        STOP
      ENDIF
C
      CLOSE(UNIT=21)
      CALL ZAIOCL(21)
C
      DO I= 1,IDM
        DO J= 1,JDM
          RSCX(I,J) = 1.0/MAX(ONEM,RSCX(I,J)) !1/qscx
          RSCY(I,J) = 1.0/MAX(ONEM,RSCY(I,J)) !1/qscy
        ENDDO
      ENDDO
C
C     INITIALIZE HYCOM INPUT AND OUTPUT.
C
      CALL ZAIOPN('OLD', 10)
      CALL ZAIOPN('OLD', 11)
      CALL ZAIOPN('NEW', 12)
      CALL ZHOPEN(10, 'FORMATTED', 'OLD', 0)
      CALL ZHOPEN(11, 'FORMATTED', 'OLD', 0)
      CALL ZHOPEN(12, 'FORMATTED', 'NEW', 0)
C
      READ( 10,'(A79)') PREAMBL
      READ( 11,'(A79)') PREAMBL
      PGRID  = INDEX(PREAMBL(2),'p-grid') .NE.0
      ROTATE = INDEX(PREAMBL(2),'E-wards').NE.0
      IF     (INDEX(PREAMBL(2),'Velocities').EQ.0) THEN
        PREAMBL(2) = 'Wind Stress Curl on p-grid'
      ELSE
        PREAMBL(2) = 'Wind Velocity Curl on p-grid'
      ENDIF
      WRITE(12,'(A79)') PREAMBL
      WRITE( 6,'(A79)') PREAMBL
C
      IF     (PGRID) THEN
        WRITE(6,*)
        IF     (ROTATE) THEN
          WRITE(6,*)
     &      'Winds first rotated and interpolated to u and v grids.'
        ELSE
          WRITE(6,*)
     &      'Winds first interpolated to u and v grids.'
        ENDIF
        WRITE(6,*)
      ELSE
        IF     (ARCTIC) THEN
          WRITE(6,'(/ a /)')
     &    'error - wind stresses must be on p-grid for tri-pole regions'
          CALL ZHFLSH(6)
          STOP
        ELSEIF (ROTATE) THEN
          WRITE(6,'(/ a /)')
     &    'error - E-ward,N-ward wind stresses must be on p-grid'
          CALL ZHFLSH(6)
          STOP
        ENDIF
        WRITE(6,*)
        WRITE(6,*) 'Wind stresses input on u and v grids.'
        WRITE(6,*)
      ENDIF
C
C     PROCESS ALL THE WIND RECORDS.
C
      DO 810 KREC= 1,999999
C
C       READ THE INPUT WIND STRESSES.
C
        READ(10,'(A)',IOSTAT=IOS) CLINE
        IF     (IOS.NE.0) THEN
          EXIT !probably end of file
        ENDIF
        WRITE(6,'(A)') CLINE
        READ(11,'(A)') CLINE
        WRITE(6,'(A)') CLINE
        CALL ZAIORD(TXM,MSK,.FALSE., XMIN,XMAX, 10)
        CALL ZAIORD(TYM,MSK,.FALSE., XMIN,XMAX, 11)
        IF     (ROTATE) THEN
C         rotate from e-ward,n-ward to x-ward,y-ward
          DO J= 1,JDM
            DO I= 1,IDM
              TXMIJ    = TXM(I,J)
              TYMIJ    = TYM(I,J)
              COSPANG  = COS(PANG(I,J))
              SINPANG  = SIN(PANG(I,J))
              TXM(I,J) = COSPANG*TXMIJ + SINPANG*TYMIJ
              TYM(I,J) = COSPANG*TYMIJ - SINPANG*TXMIJ
            ENDDO
          ENDDO
        ENDIF !rotate
        IF     (PGRID) THEN
C         interpolate from p ro u and v grids
          DO J= 1,JDM
            DO I= 1,IDM
              TXP(I) = TXM(I,J)
            ENDDO
            IF     (GLOBAL) THEN
              TXP(0) = TXP(IDM)  ! assume     periodic boundary
            ELSE
              TXP(0) = TXP(  1)  ! assume non-periodic boundary
            ENDIF
            DO I= 1,IDM
              TXM(I,J) = 0.5*(TXP(I) + TXP(I-1))
            ENDDO
          ENDDO
          DO I= 1,IDM
            DO J= 1,JDM
              TYP(J) = TYM(I,J)
            ENDDO
            TYP(0) = TYP(1)  ! assume no change across southern boundary
            DO J= 1,JDM
              TYM(I,J) = 0.5*(TYP(J) + TYP(J-1))
            ENDDO
          ENDDO
        ENDIF !pgrid
C
C       CALCULATE CURL ON Q-GRID.
C
        DO J= 1,JDM
          IF     (J.NE.1) THEN
            JM1 = J-1
          ELSE
            JM1 = 1    ! assume no change across southern boundary
          ENDIF
          DO I= 1,IDM
            IF     (I.NE.1) THEN
              IM1 = I-1
            ELSEIF (GLOBAL) THEN
              IM1 = IDM  ! assume     periodic boundary
            ELSE
              IM1 =   1  ! assume non-periodic boundary
            ENDIF
            QCURL(I,J) = RSCX(I,J)*(TYM(I,J) - TYM(IM1,J)) -
     &                   RSCY(I,J)*(TXM(I,J) - TXM(I,JM1))
          ENDDO
        ENDDO
        IF     (ARCTIC) THEN !correct northern boundary, q-grid
          J  = JDM
          JJ = JDM-1
          DO I= 1,IDM
            II = MOD(IDM-(I-1),IDM)+1
            QCURL(I,J) = QCURL(II,JJ)
          ENDDO !i
        ENDIF
C
C       INTERPOLATE CURL TO P-GRID.
C
        DO J= 1,JDM
          IF     (J.NE.JDM) THEN
            JP1 = J+1
          ELSE
            JP1 = JDM  ! assume no change across northern boundary
          ENDIF
          DO I= 1,IDM
            IF     (I.NE.IDM) THEN
              IP1 = I+1
            ELSEIF (GLOBAL) THEN
              IP1 =   1  ! assume     periodic boundary
            ELSE
              IP1 = IDM  ! assume non-periodic boundary
            ENDIF
            CURL(I,J) = 0.25*(QCURL(I,  J)   +
     &                        QCURL(IP1,J)   +
     &                        QCURL(I,  JP1) +
     &                        QCURL(IP1,JP1)  )
          ENDDO
        ENDDO
        IF     (ARCTIC) THEN !correct northern boundary, p-grid
          J  = JDM
          JJ = JDM-1
          DO I= 1,IDM
            II = IDM-MOD(I-1,IDM)
            CURL(I,J) = CURL(II,JJ)
          ENDDO !i
        ENDIF
C
C       WRITE OUT HYCOM CURL.
C
        CALL ZAIOWR(CURL,MSK,.FALSE., XMIN,XMAX, 12, .FALSE.)
        IF     (INDEX(CLINE,'month').NE.0) THEN
          WRITE(12,'(2A,1P2E16.7)') ' wndcurl',CLINE(9:26),XMIN,XMAX
          WRITE( 6,'(2A,1P2E16.7)') ' wndcurl',CLINE(9:26),XMIN,XMAX
        ELSE
          WRITE(12,'(2A,1P2E16.7)') ' wndcurl',CLINE(9:45),XMIN,XMAX
          WRITE( 6,'(2A,1P2E16.7)') ' wndcurl',CLINE(9:45),XMIN,XMAX
        ENDIF
        CALL ZHFLSH(12)
        CALL ZHFLSH( 6)
  810 CONTINUE
C
      CALL ZAIOCL(10)
      CLOSE( UNIT=10)
      CALL ZAIOCL(11)
      CLOSE( UNIT=11)
      CALL ZAIOCL(12)
      CLOSE( UNIT=12)
      STOP
C
C     END OF PROGRAM WNDCUR.
      END

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


text:
The School of Veterinary Medicine 3A page has photos that have captured much of the construction of the Veterinary Medicine Instructional Facility.

Diners at Scrubs have watched the progress of both buildings. The concrete curved tiers that will form a lecture hall on the ground floor are visible.

The Memorials memorial for California Police Dogs entitled Faithful Partner is just north of the VMIF.

http://www.news.ucdavis.edu/search/news_detail.lasso?id7713 news article
 Todays announcement was made in Gladys Valley Hall, the first building at UC Davis (and only the second in the University of California system) designed and built to achieve certification under the U.S. Green Building Councils LEED (Leadership in Energy and Environmental Design) rating system. It is expected to use onethird less energy than a standard design, plus conserve water, provide better indoor air quality and incorporate natural materials. The building is a tangible example of the innovations that will be emphasized by the Energy Efficiency Center. Less energy is consumed for lighting by employing natural light, photo sensors and motion detectors. Natural ventilation, evaporative cooling and radiant floor slabs are used to reduce cooling energy.

 When construction is finished in June, Gladys Valley Hall will become the instructional heart of the UC Davis School of Veterinary Medicine campus. The building is named for the late Gladys Valley in recognition of the longstanding generosity of the Wayne and Gladys Valley Foundation.

The building was open and lectures began Fall of 06.
Theres also a computer lab of sorts on the 2nd floor. Restricted access to the vetschool.

each room is dedicated to someone, even the elevators.
lots of cushy comfortable seating inside and out, and a snack setup system for us students and workers in the area.


text:
C    Copyright(C) 1999-2020 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    See packages/seacas/LICENSE for details

      SUBROUTINE LTADD (MDIM, MDIMO, N, LINK, IHOLD)
C***********************************************************************

C  SUBROUTINE LTADD = LOOKUP TABLE REINSURTION FOR DATA POINTER ARRAYS

C***********************************************************************

C  VARIABLES USED:
C     MDIM   = DIMENSION OF LINK ARRAY, AND BASE FOR LOOKUP START
C     LINK   = LOOKUP TABLE ARRAY OF ID'S AND POINTERS
C              LINK(1,I) = ID VALUE STORED IN I'TH ROW (0 IF EMPTY)
C              LINK(2,I) = DATA POINTER ASSOCIATED W/THIS I'TH ID VALUE
C     IHOLD  = TEMPORARY LINK ARRAY TO BE USED FOR REINSURRTION

C***********************************************************************

      DIMENSION LINK(2,MDIM), IHOLD(2,MDIM)

      LOGICAL ADDLNK

      ADDLNK=.TRUE.

C  SORT THROUGH THE ORIGINAL LINK LIST AND LINK VALID ENTRIES IN THE
C  IHOLD LIST

      DO 100 I = 1, MDIMO
         IF ( (LINK(1,I) .NE. 0) .AND. (LINK(2,I) .LE. N) )
     &      CALL LTSORT (MDIM, IHOLD, IABS(LINK(1,I)), LINK(2,I),
     &      ADDLNK)
  100 CONTINUE

C  TRANSFER THE IHOLD LIST BACK TO THE LINK LIST

      DO 120 I = 1, MDIM
         DO 110 J = 1, 2
            LINK(J,I) = IHOLD (J, I)
  110    CONTINUE
  120 CONTINUE

      RETURN

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
!
! This file is released under terms of BSD license
! See LICENSE file for more information
!
!
! Simple program to test the loop-hoist directive with reshape
! Code is toally dummy and does not perform something useful. Structure is
! important for transformation test
!

PROGRAM claw_test
  CALL claw_hoist2
END PROGRAM claw_test

SUBROUTINE claw_hoist2
  IMPLICIT NONE
  INTEGER :: jt,i,j,ntrac,klev,kproma
  REAL :: var1,var2,var3,var4,var5,time_step_len,zsedtend
  REAL, DIMENSION(10,100) :: array2d_1, array2d_2, array2d_3, array2d_4, &
                             array2d_5
  REAL, DIMENSION(10,100,2) :: array3d_1, array3d_2
  TYPE dummy
    REAL :: value1
    REAL :: value2
  END TYPE dummy
  TYPE(dummy) :: tdum

  zsedtend=1.0
  var3=2.0

  ntrac=2
  klev=10
  kproma=100
  array2d_1(:,:) = 0.0
  array2d_2(:,:) = 0.0
  array2d_3(:,:) = 0.0
  array2d_4(:,:) = 0.0
  array2d_5(:,:) = 0.0

  !$acc parallel loop gang vector collapse(2)
  DO jt=1,ntrac
    !$claw loop-hoist(j,i) reshape(array2d_2(0),array2d_1(1,2)) interchange
    IF (.FALSE.) CYCLE
    var1=tdum%value1
    var2=var1**2.
    IF (jt == 2) THEN
      DO j=1,klev
        DO i=1,kproma
          array2d_2(i,j)=MIN(tdum%value2*2,50.E-6)
        END DO
      END DO
    END IF
    DO j=1,klev
      DO i=1,kproma
        array2d_1(i,j) = 0.0
        IF (array2d_2(i,j) > 0.0) THEN
          var3 = array3d_1(i,j,jt)+array3d_2(i,j,jt)*time_step_len
          var5 = MIN(2.0/9.0, array2d_1(i,j))
          zsedtend = MIN(MAX(0.0, var3 * var5 / array2d_3(i,j)), 2.0)
          array3d_2(i,j,jt) = array3d_2(i,j,jt) - var4
          array2d_1(i,j) = var4 * array2d_4(i,j)
        END IF
      END DO
    END DO
    DO j=2,klev
      DO i=1,kproma
        array3d_2(i,j,jt) = array3d_2(i,j,jt) + (array2d_1(i,j-1) / array2d_5(i,j))
      END DO
    END DO
    !$claw end loop-hoist
    DO j=2,klev
      DO i=1,kproma
        array3d_2(i,j,jt) = array3d_2(i,j,jt) + (5.0 / array2d_5(i,j))
      END DO
    END DO
  END DO
  !$acc end parallel

END SUBROUTINE claw_hoist2

max_stars_repo_name: peclatj/claw-compiler

max_issues_repo_name: peclatj/claw-compiler

max_forks_repo_name: peclatj/claw-compiler


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2021 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 distattach_1d(xig,pneigh,pnode,a,p,ratio,nterms)
!
!     calculates the distance between the node with coordinates
!     in "pnode" and the node with local coordinate xig
!     on a line described by "nterms" nodes with coordinates
!     in pneigh
!
      implicit none
!
      integer nterms,i,j
!
      real*8 ratio(3),pneigh(3,*),pnode(3),a,xig,p(3)
!
!
!
      if(nterms.eq.2) then
         ratio(1)=(1.d0-xig)/2.d0
         ratio(2)=(1.d0+xig)/2.d0
      elseif(nterms.eq.3) then
         ratio(1)=xig*(xig-1.d0)/2.d0
         ratio(2)=(1.d0-xig)*(1.d0+xig)
         ratio(3)=xig*(xig+1.d0)/2.d0
      else
         write(*,*) '*ERROR in distattach1d: case with ',nterms
         write(*,*) '       terms is not covered'
         call exit(201)
      endif
!
!     calculating the position in the face
!
      do i=1,3
         p(i)=0.d0
         do j=1,nterms
            p(i)=p(i)+ratio(j)*pneigh(i,j)
         enddo
      enddo
!
!     calculating the distance
!
      a=(pnode(1)-p(1))**2+(pnode(2)-p(2))**2+(pnode(3)-p(3))**2
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


! Objetivo: Usar pares ordenados a partir de um gerador aleatorio do
!           tipo rapido e "sujo" para obter coordenadas de pontos em
!	    	um espaço bidimensional.

PROGRAM aleatorio_sujo_xy
	IMPLICIT none
	INTEGER, PARAMETER :: im = 6075, ia = 106, ic = 1283
	REAL, PARAMETER :: im_real = REAL(im)
	REAL :: ran_x, ran_y
	INTEGER, PARAMETER :: dados_aleatorio_xy = 7
	INTEGER :: num_repeticoes, jran = 1
! Abrir arquivo.dat de saida de dados
	OPEN(dados_aleatorio_xy,FILE="dados_aleatorio_xy.dat")
! Essas constantes dao um periodo de comprimento máximo im. Por isso im eh a
! quantidade de numeros aleatorios gerados diferentes. O gerador produz 2*im numeros
! aleatorios pois im eh impar. O primeiro numero a repetir sera na coordenada y, pois
! im eh impar e estamos distribuindo aos pares. O primeiro numero a se repetir na
! coordenada x sera apos 2im numeros gerados.
	DO num_repeticoes = 1, im
		jran = MOD(jran*ia+ic,im)
		ran_x = REAL(jran)/im_real

		jran = MOD(jran*ia+ic,im)
		ran_y = REAL(jran)/im_real
! Escrever as coordenadas no arquivo de dados
		WRITE(dados_aleatorio_xy,*) ran_x,ran_y
	END DO
! Fechar arquivo.dat de saida de dados
	CLOSE(dados_aleatorio_xy)
END PROGRAM aleatorio_sujo_xy


text: #redirect Vickers Automotive Repair


text:
      SUBROUTINE SHCURD(IYR, IMON, IDAY)
	CALL CURTIM(JUL, ITIME)
	I = JLIYMD(JUL, IYR, IMON, IDAY)
	RETURN
	END

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib


text:
c  Part of R package KernSmooth
c  Copyright (C) 1995  M. P. Wand
c
c  Unlimited use and distribution (see LICENCE).

cccccccccc FORTRAN subroutine linbin2D.f cccccccccc

c Obtains bin counts for bivariate data
c via the linear binning strategy. In this version
c observations outside the mesh are ignored.

      subroutine lbtwod(X,n,a1,a2,b1,b2,M1,M2,gcnts)
      integer n,M1,M2,i,li1,li2,ind1,ind2,ind3,ind4
      double precision X(*),a1,a2,b1,b2,gcnts(*)
      double precision lxi1,lxi2,delta1,delta2,rem1,rem2

c     Initialize grid cnts to zero

      do 10 i = 1,(M1*M2)
         gcnts(i) = dble(0)
10    continue

      delta1 = (b1 - a1)/(M1 - 1)
      delta2 = (b2 - a2)/(M2 - 1)
      do 20 i = 1,n
         lxi1 = ((X(i) - a1)/delta1) + 1
         lxi2 = ((X(n+i) - a2)/delta2) + 1

c        Find the integer part of "lxi1" and "lxi2"

         li1 = int(lxi1)
         li2 = int(lxi2)
         rem1 = lxi1 - li1
         rem2 = lxi2 - li2
         if (li1.ge.1) then
            if (li2.ge.1) then
               if (li1.lt.M1) then
                  if (li2.lt.M2) then
                     ind1 = M1*(li2-1) + li1
                     ind2 = M1*(li2-1) + li1 + 1
                     ind3 = M1*li2 + li1
                     ind4 = M1*li2 + li1 + 1
                     gcnts(ind1) = gcnts(ind1)+(1-rem1)*(1-rem2)
                     gcnts(ind2) = gcnts(ind2)+rem1*(1-rem2)
                     gcnts(ind3) = gcnts(ind3)+(1-rem1)*rem2
                     gcnts(ind4) = gcnts(ind4)+rem1*rem2
                  endif
               endif
            endif
         endif
20    continue

      return
      end

cccccccccc End of linbin2D.f cccccccccc


text:
! ============================ Module Description ===========================
! etdrk4_mod Implementation of ETDRK4 by Cox & Mathews. Contains subroutines:
!   1. etdrk4      - ETDRK4 timestepping code
!   2. initETDC    - initializes ETD coefficients
! ===========================================================================

module etdrk4_mod

    use tools_mod,  only: dp, tic, toc
    use phi_mod,    only: phi
    implicit none

    contains

    ! =======================================================================
    ! ETDRK4   Implements ETDRK4 Numerical Time Integrator
    !
    ! Arguments
    !
    !   L       (input) COMPLEX*16 array, dimension(N)
    !           vector cooresponding to PDE linear operator
    !
    !   N       (input) SUBROUTINE:
    !           cooresponds to PDE nonlinear operator. Must be of the form:
    !               subroutine N(t,y_in,N_out)
    !                   real(dp),    intent(in)  :: t
    !                   complex(dp), intent(in)  :: y_in(:)
    !                   complex(dp), intent(out) :: N_out(:)
    !               end subroutine N
    !
    !   y_in    (input) COMPLEX*16 array, dimension (N)
    !           Initial condition.
    !
    !   y_out   (output) COMPLEX*16 array, dimension (N)
    !           Initial condition.
    !
    !   times   (output) DOUBLE array, dimension(4)
    !           times(1:2) = [cputime,clocktime] for timestepping loop
    !           times(3:4) = [cputime,clocktime] for initializing coefficients,
    ! =======================================================================

    subroutine etdrk4(L,N,tspan,y_in,Nt,y_out,times)

        ! Arguments
        complex(dp), intent(in)     :: L(:)
        real(dp),    intent(in)     :: tspan(2)
        complex(dp), intent(in)     :: y_in(:)
        integer,         intent(in) :: Nt
        complex(dp), intent(out)    :: y_out(size(y_in))
        real(dp),    intent(out)    :: times(4)
        ! define interface subroutine function N
        interface
            subroutine N(t,yh_in, N_out, thread_id)
            import :: dp
            complex(dp), intent(in)  :: t
            complex(dp), intent(in)  :: yh_in(:)
            complex(dp), intent(out) :: N_out(:)
            integer, intent(in), optional :: thread_id
            end subroutine N
        end interface

        ! Local Variables
        integer :: i,nL
        real(dp) :: h
        complex(dp) :: t
        complex(dp), allocatable :: y(:)
        complex(dp), allocatable, dimension(:) :: E, E2, A123, b1, b23, b4, Nv, Na, Nb, Nc, a, b, c

        ! Initialize Variables
        nL = size(L)
        allocate(E(nL), E2(nL), A123(nL), b1(nL), b23(nL), b4(nL), Nv(nL), Na(nL), Nb(nL), Nc(nL), a(nL), b(nL), c(nL))
        y = y_in
        ! Set stepsize
        h = (tspan(2)-tspan(1))/(real(Nt,8))
        t = tspan(1)
        ! initialize ETD coefficients
        call tic()
        call initETDC(L,h,E, E2, A123, b1, b23, b4)
        call toc(times(3:4))
        ! timestepping loop
        call tic
        do i=1,Nt
            call N(t,y,Nv);
            a  = E2*y + 0.5_dp*A123*Nv;
            call N(t+h/2.0_dp, a, Na);
            b  = E2*y + 0.5_dp*A123*Na;
            call N(t+h/2.0_dp, b, Nb);
            c  = E2*a + A123*(Nb - 0.5*Nv);
            call N(t+h, c, Nc);
            y  = E*y + b1*Nv + b23*(Na + Nb) + b4*Nc;
            t = t + h;
        enddo
        call toc(times(1:2))
        y_out = y

    end subroutine etdrk4

    ! =======================================================================
    ! INITETDC   Initializes ETD coefficients for ETDRK4 scheme.
    !
    ! Arguments
    !
    !   L   (input) COMPLEX*16 array, dimension(N)
    !       vector/scaler cooresponding to PDE linear operator
    !
    !   h   (input) DOUBLE
    !       timestep
    !
    !   E, E2, A, b1, b23, b4   (output) COMPLEX*16 array, dimension(N,N)
    !                           coefficients for ETDRK4
    ! =======================================================================

    subroutine initETDC(L,h,E, E2, A, b1, b23, b4)

        ! Arguments
        complex(dp), intent(in) :: L(:)
        real(dp),    intent(in) :: h
        complex(dp), dimension(:), intent(out) :: E, E2, A, b1, b23, b4

        ! Local Variables
        complex(dp), allocatable :: P(:,:)
        allocate(P(4,size(L)))

        ! Initialize Coefficients
        call phi(h/2.0_dp * L,1,P(1:2,:))
            E2  = p(1,:)
            A   = h*p(2,:)
        call phi(h * L,3,P)
            E   = P(1,:)
            b1  = h*(4.0_dp*P(4,:) - 3.0_dp*P(3,:) + P(2,:))
            b23 = h*(-4.0_dp*P(4,:) + 2.0_dp*P(3,:))
            b4  = h*(4.0_dp*P(4,:) - P(3,:))
        deallocate(P)

    end subroutine initETDC

end module etdrk4_mod

max_stars_repo_name: buvoli/epbm

max_issues_repo_name: buvoli/epbm

max_forks_repo_name: buvoli/epbm


        !COMPILER-GENERATED INTERFACE MODULE: Fri May 22 16:59:12 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE PREPARE_ASSEMBLE_MPI__genmod
          INTERFACE
            SUBROUTINE PREPARE_ASSEMBLE_MPI
              USE SPECFEM_PAR, ONLY :                                   &
     &          NSPEC,                                                  &
     &          IBOOL,                                                  &
     &          KNODS,                                                  &
     &          NGNOD,                                                  &
     &          NGLOB,                                                  &
     &          ELASTIC,                                                &
     &          POROELASTIC,                                            &
     &          NINTERFACE,                                             &
     &          MY_NELMNTS_NEIGHBOURS,                                  &
     &          MY_INTERFACES,                                          &
     &          IBOOL_INTERFACES_ACOUSTIC,                              &
     &          IBOOL_INTERFACES_ELASTIC,                               &
     &          IBOOL_INTERFACES_POROELASTIC,                           &
     &          NIBOOL_INTERFACES_ACOUSTIC,                             &
     &          NIBOOL_INTERFACES_ELASTIC,                              &
     &          NIBOOL_INTERFACES_POROELASTIC,                          &
     &          INUM_INTERFACES_ACOUSTIC,                               &
     &          INUM_INTERFACES_ELASTIC,                                &
     &          INUM_INTERFACES_POROELASTIC,                            &
     &          NINTERFACE_ACOUSTIC,                                    &
     &          NINTERFACE_ELASTIC,                                     &
     &          NINTERFACE_POROELASTIC,                                 &
     &          MASK_ISPEC_INNER_OUTER,                                 &
     &          NIBOOL_INTERFACES_EXT_MESH,                             &
     &          IBOOL_INTERFACES_EXT_MESH_INIT
            END SUBROUTINE PREPARE_ASSEMBLE_MPI
          END INTERFACE
        END MODULE PREPARE_ASSEMBLE_MPI__genmod


text: #redirect Davis Joint Unified School District


text:
! To perform mesh motion using RBFs, we need to supply the values of displacement and velocity
! at the RBF source nodes contained in an rbf_node_patch_mm_t. This might come from evaluating a
! specified prescribed mesh motion function, or in a practical FSI simulation, from an external
! structural solver.
!
! To accommodate these different scenarios, we start here with an abstract derived type that
! performs the "update" task - that is, to update the dnodes and vnodes at RBF source nodes.
! Extensions of this abstract derive type will implement procedures to accomplish this in
! different ways.
!
!
module type_rbf_mm_driver
#include <messenger.h>
    use mod_kinds,              only: ik, rk
    implicit none

    !>
    !!
    !!
    !!
    !---------------------------------------------------------------------
    type, public, abstract :: rbf_mm_driver_t
        character(len=:),       allocatable, private    :: name_
    contains

        procedure(init_interface),          deferred    :: init             !< Initialize function and register options
        procedure(compute_disp_interface),  deferred    :: compute_disp     !< Elemental function definition
        procedure(compute_vel_interface),   deferred    :: compute_vel      !< Elemental function definition

        procedure                                       :: set_name         !< Add function name. Only use for initialization
        procedure                                       :: get_name         !< Return the function name

    end type rbf_mm_driver_t
    !*********************************************************************

    abstract interface
        subroutine init_interface(self)
            import rbf_mm_driver_t

            class(rbf_mm_driver_t),  intent(inout)  :: self

        end subroutine
    end interface

    abstract interface
        function compute_disp_interface(self, time, node) result(val)
            use mod_kinds,  only: rk
            import rbf_mm_driver_t
            class(rbf_mm_driver_t),  intent(inout)  :: self
            real(rk),                intent(in)     :: time
            real(rk),                intent(in)     :: node(3)
            real(rk)                                :: val(3)
        end function
    end interface

    abstract interface
        function compute_vel_interface(self, time, node) result(val)
            use mod_kinds,  only: rk
            import rbf_mm_driver_t
            class(rbf_mm_driver_t),  intent(inout)  :: self
            real(rk),                intent(in)     :: time
            real(rk),                intent(in)     :: node(3)
            real(rk)                                :: val(3)
        end function
    end interface

contains

    !>  Add a name for the function
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/4/2016
    !!
    !!
    !--------------------------------------------------------------------------------------------
    subroutine set_name(self,fname)
        class(rbf_mm_driver_t),  intent(inout)   :: self
        character(*),       intent(in)      :: fname

        self%name_ = fname

    end subroutine set_name
    !********************************************************************************************

    !>  Return the name of the function.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/4/2016
    !!
    !-----------------------------------------------------------------------------------------
    function get_name(self) result(fname)
        class(rbf_mm_driver_t),  intent(inout)   :: self

        character(len=:),   allocatable :: fname

        fname = self%name_

    end function get_name
    !*****************************************************************************************

end module type_rbf_mm_driver

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
      INTEGER FUNCTION NPREROC( N, NB, IPROC, ISRCPROC, NPROCS )
*
*  -- ScaLAPACK tools routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     May 1, 1997
*
*     .. Scalar Arguments ..
      INTEGER              IPROC, ISRCPROC, N, NB, NPROCS
*     ..
*
*  Purpose
*  =======
*
*  NPREROC computes the Number of PREceeding Rows Or Columns of a
*  distributed matrix that are possessed by processes closer to
*  ISRCPROC than IPROC.  Therefore, if ISRCPROC=0 and IPROC=4, then
*  NPREROC returns the number of distributed matrix rows or columns
*  owned by processes 0, 1, 2, and 3.
*
*  Arguments
*  =========
*
*  N         (global input) INTEGER
*            The number of rows or columns in the distributed matrix.
*
*  NB        (global input) INTEGER
*            Block size, size of the blocks the distributed matrix is
*            split into.
*
*  IPROC     (local intput) INTEGER
*            The coordinate of the process whose preceeding distributed
*            matrix rows or columns are to be determined.
*
*  ISRCPROC  (global input) INTEGER
*            The coordinate of the process that possesses the first
*            row or column of the distributed matrix.
*
*  NPROCS    (global input) INTEGER
*            The total number processes over which the matrix is
*            distributed.
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER              EXTRABLKS, MYDIST, NBLOCKS
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC            MOD
*     ..
*     .. Executable Statements ..
*
*     Figure PROC's distance from source process
*
      MYDIST = MOD( NPROCS+IPROC-ISRCPROC, NPROCS )
*
*     Figure the total number of whole NB blocks N is split up into
*
      NBLOCKS = N / NB
*
*     Figure the minimum number of rows/cols previous processes could have
*
      NPREROC = (NBLOCKS/NPROCS) * NB * MYDIST
*
*     See if there are any extra blocks
*
      EXTRABLKS = MOD( NBLOCKS, NPROCS )
*
*     If I have an extra block, all processes in front of me got one too
*
      IF( MYDIST.LE.EXTRABLKS ) THEN
         NPREROC = NPREROC + NB*MYDIST
*
*     If I have don't have an extra block, add in extra blocks of
*     preceeding processes and the partial block, if it exists
*
      ELSE
         NPREROC = NPREROC + EXTRABLKS*NB + MOD( N, NB )
      END IF
*
      RETURN
*
*     End of NPREROC
*
      END

max_stars_repo_name: jschueller/scalapack

max_issues_repo_name: jschueller/scalapack

max_forks_repo_name: jschueller/scalapack


module mesh

    implicit none

    real(8),allocatable :: XYZB(:,:),DXYZB(:,:)
    !xyzb => (3,node_id)  node data
    !dxyzb => (3,nrml_id) derivative data
    integer,allocatable :: NCONB(:,:),NCONDB(:,:)
    ! nconb => node list body mesh
    ! ncondb => normal list body
    real(8),allocatable :: XYZE(:,:,:),DXYZE(:,:,:),TXYZE(:,:,:)
    ! xyze =>  before combine, new full mesh for node (3,8,elem_id)
    ! dxyze => before combine, new full mesh for normal (3,8,elem,id)
    real(8),allocatable :: XYZTP(:,:),DXYZTP(:,:)
    ! xyztp => convsb, combined node mesh, xyztp(3,node_id)
    ! dxyztp = > combined normal mesh, dxyztp(3,nrml_id)
    integer,allocatable :: NCN(:),NCON(:,:),NCOND(:,:),IETYPE(:)
    ! ncn => elem_type ,namely, the node number in a elem
    ! ncon => combined node list,ncon(8,elem_id)
    ! ncond => combined nrml list,ncond(8,elem_id)
    ! ietype => flag show if a elem is free surface mesh or a body mesh

    integer,allocatable :: NNORMN(:)
    real(8),allocatable :: xyz(:,:),dxyz(:,:)
    integer :: nsys,nelem,nnode,nnoded,isys
    ! nsys => about symmetr
    ! nelem => elem number in combined mesh
    ! nnode => node number in combined mesh
    ! nnoded => nrml nmber in combined mesh
    integer :: nelemb,nnb,nnbd,nelemf,nnf
    ! nelemb => elem num in body
    ! nnb => node num  in body
    ! nnbd => nrml num in body
    ! nelemf => elem num in fs
    ! nnf => node num in fs

    real(8),protected :: wk,h,beta,amp
    ! wave information

    integer,allocatable :: nodele(:,:),nodnoe(:),nodelj(:,:),nodqua(:)
contains

    subroutine read_mesh()
        implicit none

        INTEGER IFWKO,IPOL
      !INTEGER NTnum,IFLAG_T,IND,

        REAL*8  WL
      !REAL*8  FAMPR(6),FAMPI(6),FORCER(6),FORCEI(6)
      !REAL*8  PL_AMP(6),FORAMP
      !REAL*8  FCD_AMR,  FCD_AMI

        OPEN(1, FILE='INPUT/DATIN.txt',      STATUS='OLD')
        OPEN(2, FILE='INPUT/DATBDMS.txt',    STATUS='OLD')
        OPEN(3, FILE='INPUT/DATWPMS.txt',    STATUS='OLD')

        READ(1,*)      IFWKO
        READ(1,*)      H, AMP, WK, BETA
        !==================body mesh================================
        READ(2,*)   ISYS
        READ(2,*)   NELEMB, NNB, NNBD, IPOL

        IF(ISYS.EQ.0) NSYS=1
        IF(ISYS.EQ.1) NSYS=2
        IF(ISYS.EQ.2) NSYS=4

        READ(3,*)   NELEMF

        ! =========================================

        NELEM=NELEMB+NELEMF

        ALLOCATE (NCONB(NELEMB,8),NCONDB(NELEMB,8))!body node/normal list
        ALLOCATE (XYZB(3,NNB),DXYZB(3,NNBD))    !body node/normal

        ALLOCATE(NCN(NELEM),NCON(NELEM,8),NCOND(NELEM,8),IETYPE(NELEM))

        allocate(NNORMN(8*NELEM) )

        ALLOCATE( XYZE(3,8,NELEM),DXYZE(3,8,NELEM))
        !allocate(DAMPE(8,NELEM))
        ALLOCATE( TXYZE(3,8,NELEM))
        ALLOCATE( XYZTP(3,8*NELEM),DXYZTP(3,8*NELEM))!,DAMPTP(8*NELEM))

        !ALLOCATE(SAMB(NELEM,16,0:8),SAMBXY(NELEM,16,3),
        !&         DSAMB(NELEM,16,6))

        call MESHFS4()! Read in data on free surface mesh
        call MESHBD(IPOL) ! Read in data on body mesh

        close(2)
        close(1)
        close(3)

    end subroutine

    subroutine prepare_mesh()
        implicit none

        integer :: ind,node_max
        integer :: j,inode,l,ielem

        node_max = max(nnode,nnoded)

        allocate(xyz(3,node_max),dxyz(6,node_max))

        do ind = 1,nnode
            xyz(1:3,ind) = xyztp(1:3,ind)
        end do
       ! didn't assing damptp
        do ind = 1,nnoded
            dxyz(1:3,ind) = dxyztp(1:3,ind)
        end do

        deallocate(xyztp,dxyztp)

        !==============================================
        allocate(nodele(nnode,64),nodnoe(nnode),nodelj(nnode,64)&
                &,nodqua(nnode))

        DO 50 INODE=1, NNODE
        L=0
        DO 40 IELEM=1,  NELEM
        DO 30 J=1,      NCN(IELEM)
        IF(INODE.EQ.NCON(IELEM,J)) THEN
        L=L+1
        NODELE(INODE,L)=IELEM
        NODELJ(INODE,L)=J
        ENDIF
30      CONTINUE
40      CONTINUE
        NODNOE(INODE)=L
!
        NODQUA(INODE)=0
        IF( NSYS .GE. 2) THEN
          IF( DABS(XYZ(2,INODE)).LT.1.0E-06 ) THEN
          NODQUA(INODE)=2
          END IF
        END IF
!
        IF( NSYS .EQ. 4) THEN
          IF( DABS(XYZ(1,INODE)).LT.1.0E-06.AND.&
             & DABS(XYZ(2,INODE)).LT.1.0E-06) THEN
           NODQUA(INODE)=5
          ELSE IF( DABS(XYZ(1,INODE)).LT.1.0E-06 ) THEN
           NODQUA(INODE)=4
          ENDIF
        END IF
!
50      CONTINUE
    end subroutine
!       MESHFS4 + MESHBD
!
!C *******************************************************************
!C *                                                                 *
!C *  Read in the data file of the mesh on the free surface          *
!C *                                                                 *
!C *                                                                 *
!C *******************************************************************
!C
    subroutine meshfs4()

        implicit none

        integer :: IE,J,M

        XYZE(3,1:8, 1:NELEMF)=0.0d0

        DO 100 IE=1, NELEMF
            IETYPE(IE)=2
            READ(3, *)    M, NCN(IE)
            READ(3, *) (XYZE(1,J,IE), J=1, NCN(IE))
            READ(3, *) (XYZE(2,J,IE), J=1, NCN(IE))

      !dampe(:,:)=0.0d0

100    CONTINUE
        ! in dat file ,positive nrml is set as pointing into fulid field
      DXYZE(1, 1:8, 1:NELEMF)= 0.0d0
      DXYZE(2, 1:8, 1:NELEMF)= 0.0d0
      DXYZE(3, 1:8, 1:NELEMF)= 1.0d0
!     DXYZE(3, 1:8, 1:NELEMF)=-1.0d0

      RETURN
      END

!C ***************************************************************
!C *                                                             *
!C *  Generate nodal and element data on the body surface of     *
!C *  an arbitrary body                                          *
!C *                                                             *
!C ***************************************************************
!C
        SUBROUTINE MESHBD(NCOR)

!     USE MVAR_MOD
        IMPLICIT   NONE

!     INTEGER NELEM,NELEMB,NELEMF,NNODE,NNODED,NNB,NNBD,NNF,NNTCH

        INTEGER NCOR,IPL,IPOLAR(50)
        INTEGER I,IE,M,INODE,NCNN,K,KK
        REAL*8  A1,A2,R1,R2,Z1
        REAL*8  XOFSET(50),YOFSET(50),ZOFSET(50)
        real(8),parameter :: pi = 3.141592653589793
        !
        ! --------------------------------------------
        DO 5 I=1, NCOR
        READ(2,*)  M, IPOLAR(I), XOFSET(I),YOFSET(I),ZOFSET(I)
5       CONTINUE
!
        DO 10 INODE=1, NNB
        READ(2,*) M, IPL, R1, A1, Z1
        IF(IPOLAR(IPL).EQ.0) THEN
         XYZB(1,INODE)=R1+XOFSET(IPL)
         XYZB(2,INODE)=A1+YOFSET(IPL)
        ELSE IF(IPOLAR(IPL).EQ.1) THEN
         XYZB(1,INODE)=R1*DCOS(A1*PI/180.0D0)+XOFSET(IPL)
         XYZB(2,INODE)=R1*DSIN(A1*PI/180.0D0)+YOFSET(IPL)
        ENDIF
       XYZB(3,INODE)=Z1+ZOFSET(IPL)
10      CONTINUE
!
        DO 11 INODE=1, NNBD
        READ(2,*) M,IPL,R2,A2,DXYZB(3,INODE)

        IF(IPOLAR(IPL).EQ.0) THEN
          DXYZB(1,INODE)= R2
          DXYZB(2,INODE)= A2
        ELSE IF(IPOLAR(IPL).NE.0) THEN
          DXYZB(1,INODE)=R2*DCOS( A2*PI/180.0D0 )
          DXYZB(2,INODE)=R2*DSIN( A2*PI/180.0D0 )
        ENDIF
11      CONTINUE

        DO 20 I=1, NELEMB
      IE=I+NELEMF
      IETYPE(IE)=1
        READ(2,*) M, NCN(IE)
        READ(2,*) (NCONB(I,K), K=1, NCN(IE))
20    CONTINUE
!
     DO 30 I=1,  NELEMB
      IE=I+NELEMF
        READ(2,*)  M, NCNN
        READ(2,*) (NCONDB(I,K), K=1, NCN(IE))
30     CONTINUE
!
! ===============================================================
!    ×ª»»ÍøžñÊýŸÝžñÊœ£¬œ«Ë®ÏßÉÏœÚµãŽÓÎïÃæÉÏÏû³ý£¬¹éÈëµœË®ÃæÉÏ
!
       DO 100 I=1, NELEMB
      IE=I+NELEMF
      DO 80 K=1, NCN(IE)
         XYZE(1,K,IE)=XYZB(1, NCONB(I,K))
         XYZE(2,K,IE)=XYZB(2, NCONB(I,K))
         XYZE(3,K,IE)=XYZB(3, NCONB(I,K))
80      CONTINUE
!
100    CONTINUE
!
       DO 200 I=1, NELEMB
      IE=I+NELEMF
      DO 180 K=1, NCN(IE)
         DXYZE(1,K,IE)=DXYZB(1,NCONDB(I,K))
         DXYZE(2,K,IE)=DXYZB(2,NCONDB(I,K))
         DXYZE(3,K,IE)=DXYZB(3,NCONDB(I,K))
180     CONTINUE
!
200    CONTINUE
!
       END

    SUBROUTINE CONVSB

       !USE MVAR_mod
     IMPLICIT NONE

       INTEGER IEB,IEL,IND,I,L,M,NND0,kk
     REAL(8) X,Y,Z,DX,DY,DZ,DR2,tmp1,tmp2,tmp3,tol

        !C ------------------------------------------------------------
        !C  NELEM : total number of elements
        !C  NELEMF: total number of elements
        !C  NELEMB: total number of elements
        !C
        !C  NNODE : total number of nodes
        !C  NNF   : total number of nodes on the free surface
        !C  NNB   : total number of nodes on the body surface
        !C ------------------------------------------------------------
        !C
        !C ** check nodes' situation

        !WRITE(6,*) '  Inside CONVSB'
        !write(11,*) '  Inside CONVSB'
        !write(11,*) '  NELEMF=',NELEMF, '  NELEMB=',NELEMB

        !DO  IEL=1, NELEMF+NELEMB
         !   write(11,1000) IEL,NCN(IEL)
         !   write(11,1100) (XYZE(1,K,IEL),K=1,NCN(IEL))
         !   write(11,1100) (XYZE(2,K,IEL),K=1,NCN(IEL))
         !   write(11,1100) (XYZE(3,K,IEL),K=1,NCN(IEL))
        !ENDDO

1100    FORMAT(8E15.6)

      M=0

      IF (NELEMF .EQ. 0) THEN
        NNF=0
      ELSE

        DO 200 IEL=1,  NELEMF
            DO 180 I=1, NCN(IEL)

                X=XYZE(1,I,IEL)
                Y=XYZE(2,I,IEL)
                Z=XYZE(3,I,IEL)
                !========check coincident node
                DO 150 L=1, M
                    tmp1 = DABS( X-XYZTP(1,L) )
                    tmp2 = DABS( Y-XYZTP(2,L) )
                    tmp3 = DABS( Z-XYZTP(3,L) )
                    tol = 1.0E-04
               IF((tmp1.LT.tol).AND.(tmp2.LT.tol).AND.(tmp3.lt.tol)) THEN
                 NCON(IEL, I)=L
                 GOTO 180
                    END IF
 150           CONTINUE

!C ** if the point (XPM, YPM) does not concide with any other
!C    points had been checked,  give a new code
!C
             M=M+1
             NCON(IEL, I)=M
             ! a new node matrix,xyztp
             XYZTP(1,M)=X
             XYZTP(2,M)=Y
             XYZTP(3,M)=Z

180       CONTINUE
200     CONTINUE

       NNF=M
       END IF
!
! -----------------------------------------
!
         DO 300 IEB=1,  NELEMB
          IEL=IEB+NELEMF
           DO 280 I=1, NCN(IEL)

         X=XYZE(1,I,IEL)
         Y=XYZE(2,I,IEL)
         Z=XYZE(3,I,IEL)

             DO 250 L=1, M
               tmp1 = DABS( X-XYZTP(1,L) )
        tmp2 = DABS( Y-XYZTP(2,L) )
        tmp3 = DABS( Z-XYZTP(3,L) )
        tol = 1.0E-04
               IF((tmp1.LT.tol).AND.(tmp2.LT.tol).AND.(tmp3.lt.tol)) THEN
          NCON(IEL, I)=L
                 GOTO 280
               END IF
 250         CONTINUE
!C
!C ** if the point (XPM, YPM) does not concide with any other
!C    points had been checked,  give a new code
!C
             M=M+1
             NCON(IEL, I)=M
             XYZTP(1,M)=X
             XYZTP(2,M)=Y
             XYZTP(3,M)=Z

280       CONTINUE
300     CONTINUE

       NNODE=M
     NNB=NNODE-NNF
    ! write(11,*)  ' NNB=',NNB,' NNF=',NNF,' NNODE=',NNODE
     !write(11,*)
!
! ===========================================================
!
!
       NNODED=NNODE
!
!  对每个节点先找出一个方向向量
!
       DO 350 IND=1,  NNODE
         DO 320 IEL=1,  NELEM
           DO 320 I=1,  NCN(IEL)

          IF(NCON(IEL, I) .EQ. IND) THEN
               DX=DXYZE(1,I,IEL)
             DY=DXYZE(2,I,IEL)
               DZ=DXYZE(3,I,IEL)
             NCOND(IEL, I)=IND
             DXYZTP(1,IND)=DX
             DXYZTP(2,IND)=DY
             DXYZTP(3,IND)=DZ
             NNORMN(IND)=IND
             GOTO 350
         ENDIF
320   CONTINUE
350  CONTINUE
!
!  对每个节点查找是否还有其他方向向量
!
     DO 500 IND=1,  NNODE

      NND0=NNODED

         DO 480 IEL=1,  NELEM
           DO 450 I=1,  NCN(IEL)

          IF(NCON(IEL, I)  .EQ. IND) THEN
            DX=DXYZE(1,I,IEL)
            DY=DXYZE(2,I,IEL)
            DZ=DXYZE(3,I,IEL)

            DR2=((DX-DXYZTP(1,IND))**2+ &
                 (DY-DXYZTP(2,IND))**2+&
                 (DZ-DXYZTP(3,IND))**2)

           IF( DR2 .LT. 1.0E-6)  THEN
             NCOND(IEL, I)=IND
             GOTO 450
           ENDIF

           DO KK=NND0+1,  NNODED    ! 检查是否与已定义的法向量相同
             DR2=(( DX-DXYZTP(1,KK))**2+&
                  ( DY-DXYZTP(2,KK))**2+&
                  ( DZ-DXYZTP(3,KK))**2)

            IF( DR2 .LT. 1.0E-6) THEN
             NCOND(IEL, I)=KK
             GOTO 450
            ENDIF

           ENDDO
!

              NNODED=NNODED+1
            NCOND(IEL, I)=NNODED
            DXYZTP(1,NNODED)=DX
            DXYZTP(2,NNODED)=DY
            DXYZTP(3,NNODED)=DZ

            XYZTP(1,NNODED)=XYZTP(1,IND)
            XYZTP(2,NNODED)=XYZTP(2,IND)
            XYZTP(3,NNODED)=XYZTP(3,IND)
              NNORMN(NNODED)=IND

          ENDIF

!
450       CONTINUE
!

480      CONTINUE
500     CONTINUE

       END

end module


text:
      subroutine fem_update_stress(i, j, iUnload)
!***********************************************************************
! Copyright 2011 Los Alamos National Security, LLC  All rights reserved
! Unless otherwise indicated,  this information has been authored by an
! employee or employees of the Los Alamos National Security, LLC (LANS),
! operator of the  Los  Alamos National  Laboratory  under Contract  No.
! DE-AC52-06NA25396  with  the U. S. Department  of  Energy.  The  U. S.
! Government   has   rights  to  use,  reproduce,  and  distribute  this
! information.  The  public may copy  and  use this  information without
! charge, provided that this  Notice and any statement of authorship are
! reproduced on all copies.  Neither  the  Government nor LANS makes any
! warranty,   express   or   implied,   or   assumes  any  liability  or
! responsibility for the use of this information.
!***********************************************************************
!
! Computes stresses from strains in three-dimensional 'fem' computations
!
! Author : Sai Rapaka
!

      use comai, only: ns,iad,iout
      use comdi, only: t, phi, tini, phini
      use comsi, only: e1, e2, e3, du, dv, dw, alp, bulk, iPlastic
      use comsi, only: delta_u, delta_v, delta_w
      use comfem

      implicit none

      integer                      :: i,j
      integer                      :: k
      integer, dimension(8)        :: node
      real*8,  dimension(8)        :: alpha, beta, deltaT, deltaP
      real*8,  dimension(6, 24)    :: B
      real*8,  dimension(6, 6)     :: D
      real*8,  dimension(6)        :: gp_stress, gp_strain
      real*8,  dimension(24)       :: disp
      logical                      :: iUnload

      real*8                       :: e1bar, e2bar, e3bar
      real*8                       :: alphadeltaT, betadeltaP
      real*8                       :: sijeij

      ! Initialize unloading flag to false
      iUnload = .false.
      ! first compute the strain
      do k=1,ns
        node(k) = elnode(i,k)
        disp(3*k-2) = delta_u(node(k))
        disp(3*k-1) = delta_v(node(k))
        disp(3*k  ) = delta_w(node(k))
        alpha(k) = alp(node(k))
        deltaT(k) = t(node(k)) - tini(node(k))
        beta(k) = bulk(node(k))
        deltaP(k) = phi(node(k)) - phini(node(k))
      enddo

      B = 0.0d0
      do k=1,ns
        B(1,3*(k-1) + 1) = dPsidX(i, j, k)
        B(2,3*(k-1) + 2) = dPsidY(i, j, k)
        B(3,3*(k-1) + 3) = dPsidZ(i, j, k)
        B(4,3*(k-1) + 1) = dPsidY(i, j, k)
        B(4,3*(k-1) + 2) = dPsidX(i, j, k)
        B(5,3*(k-1) + 2) = dPsidZ(i, j, k)
        B(5,3*(k-1) + 3) = dPsidY(i, j, k)
        B(6,3*(k-1) + 1) = dPsidZ(i, j, k)
        B(6,3*(k-1) + 3) = dPsidX(i, j, k)
      enddo

      gp_strain = matmul(B, disp)

      CALL DoubleDot(fem_stress(i, j, :), gp_strain, sijeij)
      if(sijeij.lt.0.0d0) then
!        iUnload = .true.
      endif

      e1bar = 0.0d0
      e2bar = 0.0d0
      e3bar = 0.0d0

      alphadeltaT = 0.0d0
      betadeltaP = 0.0d0

      do k=1,ns
        e1bar = e1bar + Psi(i, j, k)*e1(node(k))
        e2bar = e2bar + Psi(i, j, k)*e2(node(k))
        e3bar = e3bar + Psi(i, j, k)*e3(node(k))

        alphadeltaT = alphadeltaT +
     &                 Psi(i,j,k)*alpha(k)*deltaT(k)
        betadeltaP = betadeltaP +
     &                 Psi(i,j,k)*beta(k)*deltaP(k)

      enddo

      if(iPlastic.eq.1) then
        call fem_material_stress_update(i, j, gp_stress, gp_strain,
     & iUnload)
      else

!        if(iad.eq.0) then
!          fem_strain(i,j,:) = gp_strain
!        else
          fem_strain(i,j,:) = conv_strain(i,j,:) + gp_strain
!        endif

        gp_strain = fem_strain(i,j,:)
        gp_strain(1) = gp_strain(1) - alphadeltaT - betadeltaP
        gp_strain(2) = gp_strain(2) - alphadeltaT - betadeltaP
        gp_strain(3) = gp_strain(3) - alphadeltaT - betadeltaP

        D = 0.0d0
        D(1,1) = e1bar
        D(1,2) = e2bar
        D(1,3) = e2bar
        D(2,1) = e2bar
        D(2,2) = e1bar
        D(2,3) = e2bar
        D(3,1) = e2bar
        D(3,2) = e2bar
        D(3,3) = e1bar
        D(4,4) = e3bar
        D(5,5) = e3bar
        D(6,6) = e3bar

        fem_stress(i,j,:) = matmul(D, gp_strain)
      endif

      ! fem_stress(i,j,:) = fem_stress(i,j,:) + gp_stress
      ! fem_strain(i,j,:) = fem_strain(i,j,:) + gp_strain
      ! fem_strain(i,j,:) = gp_strain + alphadeltaT + betadeltaP

      end subroutine fem_update_stress

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


! mesh.f90
! Assignment of charges to a 3-d mesh
PROGRAM mesh

  !------------------------------------------------------------------------------------------------!
  ! This software was written in 2016/17                                                           !
  ! by Michael P. Allen <m.p.allen@warwick.ac.uk>/<m.p.allen@bristol.ac.uk>                        !
  ! and Dominic J. Tildesley <d.tildesley7@gmail.com> ("the authors"),                             !
  ! to accompany the book "Computer Simulation of Liquids", second edition, 2017 ("the text"),     !
  ! published by Oxford University Press ("the publishers").                                       !
  !                                                                                !
  ! LICENCE                                                                                !
  ! Creative Commons CC0 Public Domain Dedication.                                                 !
  ! To the extent possible under law, the authors have dedicated all copyright and related         !
  ! and neighboring rights to this software to the PUBLIC domain worldwide.                        !
  ! This software is distributed without any warranty.                                             !
  ! You should have received a copy of the CC0 Public Domain Dedication along with this software.  !
  ! If not, see <http://creativecommons.org/publicdomain/zero/1.0/>.                               !
  !                                                                                !
  ! DISCLAIMER                                                                                !
  ! The authors and publishers make no warranties about the software, and disclaim liability       !
  ! for all uses of the software, to the fullest extent permitted by applicable law.               !
  ! The authors and publishers do not recommend use of this software for any purpose.              !
  ! It is made freely available, solely to clarify points made in the text. When using or citing   !
  ! the software, you should not imply endorsement by the authors or publishers.                   !
  !------------------------------------------------------------------------------------------------!

  ! This program assigns a set of charges to a cubic mesh using the
  ! triangular-shaped cloud distribution described by Hockney and Eastwood (1988)
  ! The charges are positioned in a box of unit length.
  ! The charge mesh is indexed from 0 to sc-1 in each coordinate direction

  USE, INTRINSIC :: iso_fortran_env, ONLY : input_unit, output_unit, error_unit, iostat_end, iostat_eor

  USE mesh_module, ONLY : mesh_function

  IMPLICIT NONE

  INTEGER                             :: n   ! Number of charges
  INTEGER                             :: sc  ! Dimension of mesh
  REAL, ALLOCATABLE, DIMENSION(:,:,:) :: rho ! Mesh cell charge density (0:sc-1,0:sc-1,0:sc-1)
  REAL, ALLOCATABLE, DIMENSION(:,:)   :: r   ! Charge positions (3,n)
  REAL, ALLOCATABLE, DIMENSION(:)     :: q   ! Charges (n)

  INTEGER :: n2, n3, ioerr
  REAL    :: h

  NAMELIST /nml/ n, sc

  WRITE ( unit=output_unit, fmt='(a)' ) 'mesh'
  WRITE ( unit=output_unit, fmt='(a)' ) '3-D mesh assignment of charges'
  WRITE ( unit=output_unit, fmt='(a)' ) 'Unit box length, coordinates in range (0,1)'

  ! Example values for illustration
  n  = 4  ! Number of charges
  sc = 10 ! Dimension of mesh

  ! Read parameters from namelist
  ! Comment out, or replace, this section if you don't like namelists
  READ ( unit=input_unit, nml=nml, iostat=ioerr )
  IF ( ioerr /= 0 ) THEN
     WRITE ( unit=error_unit, fmt='(a,i15)') 'Error reading namelist nml from standard input', ioerr
     IF ( ioerr == iostat_eor ) WRITE ( unit=error_unit, fmt='(a)') 'End of record'
     IF ( ioerr == iostat_end ) WRITE ( unit=error_unit, fmt='(a)') 'End of file'
     STOP 'Error in mesh'
  END IF

  h = 1.0 / REAL( sc ) ! Mesh spacing

  ! Write out run parameters
  WRITE ( unit=output_unit, fmt='(a,t40,i15)'   ) 'Number of charges', n
  WRITE ( unit=output_unit, fmt='(a,t40,i15)'   ) 'Dimension of mesh', sc
  WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Mesh spacing',      h

  ALLOCATE ( rho(0:sc-1,0:sc-1,0:sc-1) ) ! C-style indexing is convenient here
  ALLOCATE ( r(3,n), q(n) )

  ! For illustration we choose random charge positions with coordinates in range (0,1)
  ! In a real application, we would convert positions into this range
  CALL RANDOM_SEED() ! same random number sequence every time
  CALL RANDOM_NUMBER ( r )

  ! For illustration we choose +1 and -1 charges, alternately
  q(1::2) = 1.0
  q(2::2) = -1.0

  rho = mesh_function ( r, q, sc )

  ! Output mesh charge density
  DO n3 = 0, sc-1
     WRITE( unit=output_unit, fmt='(a,i5)' ) 'z-layer ', n3
     DO n2 = 0, sc-1
        WRITE( unit=output_unit, fmt='(*(f10.4))') rho(:,n2,n3)
     END DO
  END DO

  ! Finally check integrated charge density
  WRITE( unit=output_unit, fmt='(a,2f10.6)') 'Total charge = ', SUM ( q ), SUM ( rho ) * (h**3)

  DEALLOCATE ( r, q, rho )

END PROGRAM mesh