text:
! ###################################################################
! Copyright (c) 2013-2022, Marc De Graef Research Group/Carnegie Mellon University
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are
! permitted provided that the following conditions are met:
!
!     - Redistributions of source code must retain the above copyright notice, this list
!        of conditions and the following disclaimer.
!     - Redistributions in binary form must reproduce the above copyright notice, this
!        list of conditions and the following disclaimer in the documentation and/or
!        other materials provided with the distribution.
!     - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
!        of its contributors may be used to endorse or promote products derived from
!        this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ###################################################################
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
! EMsoft:ECCImod.f90
!--------------------------------------------------------------------------
!
! MODULE: ECCImod
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief EMECCI helper routines
!
!> @date 11/25/15 MDG 1.0 original
!---------------------------------------------------------------------------
module ECCImod

use local

IMPLICIT NONE

contains

!--------------------------------------------------------------------------
!
! SUBROUTINE: Calckvectorcone
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief compute a set of incident beam directions for single image ECCI mode
!
!> @param cell unit cell pointer
!> @param khead head of kvector linked list
!> @param k incident wave vector (zone axis)
!> @param ga principal g vector
!> @param ktxy tangential components
!> @param ktrad cone opening angle
!> @param ktstep number of steps along cone radius
!> @param numk resulting number of incident beam directions
!
!> @date 11/29/01 MDG 1.0 original
!> @date 12/05/13 MDG 2.0 adapted for ECCI simulations
!> @date 12/01/15 MDG 2.1 simplification of input parameters
!--------------------------------------------------------------------------
recursive subroutine Calckvectorcone(cell,khead,k,ga,ktxy,ktrad,ktstep,numk)
!DEC$ ATTRIBUTES DLLEXPORT :: Calckvectorcone

use io
use error
use diffraction
use crystal
use kvectors

IMPLICIT NONE

type(unitcell)                      :: cell
type(kvectorlist),pointer           :: khead
integer(kind=irg),INTENT(IN)        :: k(3)
integer(kind=irg),INTENT(IN)        :: ga(3)
real(kind=sgl),INTENT(IN)           :: ktxy(2)
real(kind=sgl),INTENT(IN)           :: ktrad
integer(kind=irg),INTENT(IN)        :: ktstep
integer(kind=irg),INTENT(OUT)       :: numk

type(kvectorlist),pointer           :: ktmp,ktail
integer                             :: istat,imin,imax,jmin,jmax,ijmax,i,j,ic,jc,ki
real                                :: kr(3),glen,delta,kstar(3),kt(3),gan(3),gperp(3),ktlen, dkt

! compute geometrical factors
 glen = CalcLength(cell,float(ga),'r')         ! length of ga
 gan = ga/glen                                 ! normalized ga
 delta = ktrad*glen/float(ktstep)              ! grid step size in nm-1
 dkt = ktrad/float(ktstep)
 call TransSpace(cell,float(k),kstar,'d','r')       ! transform incident direction to reciprocal space
 call CalcCross(cell,float(ga),kstar,gperp,'r','r',0)! compute g_perp = ga x k
 call NormVec(cell,gperp,'r')                       ! normalize g_perp
 call NormVec(cell,kstar,'r')                       ! normalize reciprocal beam vector

! deal only with the incident beam (parallel illumination)
if (ktstep.eq.0) then
 if (.not.associated(khead)) then     ! allocate the head and ktail of the linked list
   allocate(khead,stat=istat)         ! allocate new value
   if (istat.ne.0) call FatalError('Calckvectorcone: unable to allocate head pointer',' ')
   ktail => khead                      ! ktail points to new value
   nullify(ktail%next)                ! nullify next in new value
   numk = 1                          ! keep track of number of k-vectors so far
 ! this should be the center vector of the illumination cone !!!
   kt = - glen * (ktxy(1)*gan + ktxy(2) * gperp)
   ktail%kt = kt                           ! store tangential component of k
   ktlen = glen**2*(ktxy(1)**2+ktxy(2)**2)         ! squared length of tangential component

   kr = kt + sqrt(1.0/cell%mLambda**2 - ktlen)*kstar ! complete wave vector
   ktail%k = kr                            ! store in pointer list
   ktail%kn = CalcDot(cell,ktail%k,dble(kstar),'r')    ! normal component of k
 end if
else
! next, put the center of the cone in units of (i,j) (original ECP "screen" coordinates)
  ic = int(ktxy(1)*glen/delta)
  jc = int(ktxy(2)*glen/delta)
  ki = ktstep

 if (.not.associated(khead)) then     ! allocate the head and ktail of the linked list
   allocate(khead,stat=istat)         ! allocate new value
   if (istat.ne.0) call FatalError('Calckvectorcone: unable to allocate head pointer',' ')
   ktail => khead                      ! ktail points to new value
   nullify(ktail%next)                ! nullify next in new value
   numk = 1                          ! keep track of number of k-vectors so far
 ! this should be the center vector of the illumination cone !!!
   ktail%i = ic                            ! i-index of beam
   ktail%j = jc                            ! j-index of beam
   kt = -float(ktail%i)*delta*gan - float(ktail%j)*delta*gperp  ! tangential component of k
   ktail%kt = kt                           ! store tangential component of k
   ktlen = delta**2*(ktail%i**2+ktail%j**2)         ! squared length of tangential component

   kr = kt + sqrt(1.0/cell%mLambda**2 - ktlen)*kstar ! complete wave vector
   ktail%k = kr                            ! store in pointer list
   ktail%kn = CalcDot(cell,ktail%k,dble(kstar),'r')    ! normal component of k
 else
   call FatalError('Calckvectorcone: pointer head already allocated',' ')
 end if

! the following lines are quite different if symmetry is taken into account;
! check the MBsym.f90 program to determine how that can be done.
  imin =  -ki; imax = ki; jmin = -ki; jmax = ki;
  ijmax = ki**2
! now do the real work
  do i=imin,imax
   do j=jmin,jmax
    if (.not.((i.eq.0).and.(j.eq.0))) then  ! the point (0,0) has already been taken care of
     if ((i**2+j**2).le.ijmax) then   ! is point inside the incident cone ?
      allocate(ktail%next,stat=istat)  ! allocate new value
      if (istat.ne.0) call FatalError('Calckvectorcone: unable to allocate pointer',' ')
      ktail => ktail%next               ! ktail points to new value
      nullify(ktail%next)              ! nullify next in new value
      numk = numk + 1                 ! keep track of number of k-vectors so far
      ktail%i = ic+i                   ! i-index of beam
      ktail%j = jc+j                   ! j-index of beam
      kt = - float(ktail%i)*delta*gan - float(ktail%j)*delta*gperp  ! tangential component of k
      ktail%kt = kt                    ! store tangential component of k
      ktlen = delta**2*(ktail%i**2+ktail%j**2)         ! squared length of tangential component

      kr = kt + sqrt(1.0/cell%mLambda**2 - ktlen)*kstar ! complete wave vector
      ktail%k = kr                     ! store in pointer list
      ktail%kn = CalcDot(cell,ktail%k,dble(kstar),'r')    ! normal component of k
     end if
    end if
   end do
  end do
end if

end subroutine Calckvectorcone

!--------------------------------------------------------------------------
!
! SUBROUTINE: Calckvectortrace
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief compute a set of incident beam directions for line scan ECCI mode
!
!> @param cell unit cell pointer
!> @param khead head of k-vector list
!> @param k incident wave vector (zone axis)
!> @param ga principal g vector
!> @param ktxy tangential components of trace start point
!> @param ktxy2 tangential components of trace end point
!> @param ktrad cone opening angle
!> @param ktstep number of steps along cone radius
!> @param numk resulting number of incident beam directions
!
!> @date 12/08/13 MDG 1.0 original
!> @date 12/01/15 MDG 1.1 simplifcation of input variables
!--------------------------------------------------------------------------
recursive subroutine Calckvectortrace(cell,khead,k,ga,ktxy,ktxy2,ktrad,ktstep,numk)
!DEC$ ATTRIBUTES DLLEXPORT :: Calckvectortrace

use io
use error
use diffraction
use crystal
use kvectors

IMPLICIT NONE

type(unitcell)                      :: cell
type(kvectorlist),pointer           :: khead
integer(kind=irg),INTENT(IN)        :: k(3)
integer(kind=irg),INTENT(IN)        :: ga(3)
real(kind=sgl),INTENT(IN)           :: ktxy(2)
real(kind=sgl),INTENT(IN)           :: ktxy2(2)
real(kind=sgl),INTENT(IN)           :: ktrad
integer(kind=irg),INTENT(IN)        :: ktstep
integer(kind=irg),INTENT(OUT)       :: numk

type(kvectorlist),pointer           :: ktail
integer                             :: istat,j
real                                :: kr(3),glen,delta,kstar(3),kt(3),gan(3),gperp(3),ktlen, dktx, dkty

! compute geometrical factors
 glen = CalcLength(cell,float(ga),'r')              ! length of ga
 gan = ga/glen                                 ! normalized ga
 delta = 2.0*ktrad*glen/float(2*ktstep+1)      ! grid step size in nm-1
 call TransSpace(cell,float(k),kstar,'d','r')       ! transform incident direction to reciprocal space
 call CalcCross(cell,float(ga),kstar,gperp,'r','r',0)! compute g_perp = ga x k
 call NormVec(cell,gperp,'r')                       ! normalize g_perp
 call NormVec(cell,kstar,'r')                       ! normalize reciprocal beam vector

 j = 0
 if (.not.associated(khead)) then     ! allocate the head and ktail of the linked list
   allocate(khead,stat=istat)         ! allocate new value
   if (istat.ne.0) call FatalError('Calckvectortrace: unable to allocate head pointer',' ')
   ktail => khead                     ! ktail points to new value
   nullify(ktail%next)                ! nullify next in new value
   numk = 1                           ! keep track of number of k-vectors so far
! this should be the starting point of the line trace
!   kt = - glen * ( ktxy(1)*gan + ktxy(2) * gperp)
   kt = - glen * ( ktxy(1)*gan - ktxy(2) * gperp)
   ktail%kt = kt                           ! store tangential component of k
   ktlen = glen**2*(ktxy(1)**2+ktxy(2)**2)         ! squared length of tangential component
   kr = kt + sqrt(1.0/cell%mLambda**2 - ktlen)*kstar ! complete wave vector
   ktail%k = kr                            ! store in pointer list
   ktail%kn = CalcDot(cell,ktail%k,dble(kstar),'r')    ! normal component of k
 end if

 dktx = (ktxy2(1) - ktxy(1))/float(ktstep-1)
 dkty = (ktxy2(2) - ktxy(2))/float(ktstep-1)

 do j=1,ktstep-1
      allocate(ktail%next,stat=istat)  ! allocate new value
      if (istat.ne.0) call FatalError('Calckvectortrace: unable to allocate pointer',' ')
      ktail => ktail%next              ! ktail points to new value
      nullify(ktail%next)              ! nullify next in new value
      numk = numk + 1                  ! keep track of number of k-vectors so far
!     kt = - glen * ( (ktxy(1)+float(j)*dktx)*gan + (ktxy(2)+float(j)*dkty) * gperp) ! tangential component of k
      kt = - glen * ( (ktxy(1)+float(j)*dktx)*gan - (ktxy(2)+float(j)*dkty) * gperp) ! tangential component of k
      ktail%kt = kt                    ! store tangential component of k
      ktlen = delta**2*(ktail%i**2+ktail%j**2)         ! squared length of tangential component
      kr = kt + sqrt(1.0/cell%mLambda**2 - ktlen)*kstar ! complete wave vector
      ktail%k = kr                     ! store in pointer list
      ktail%kn = CalcDot(cell,ktail%k,dble(kstar),'r')    ! normal component of k
 end do

end subroutine Calckvectortrace

!--------------------------------------------------------------------------
!
! SUBROUTINE: ECCICalcSgh
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief compute structure factor-like array for ECCI and ECP simulations
!
!> @param nn dimension of array
!> @param Sgh output array
!> @param nat normalization array
!
!> @date 03/05/14  MDG 1.0 original (used to be in-line in ECP and ECCI programs)
!> @date 03/11/14  MDG 1.1 converted to diagonal Sgh array only
!> @date 10/22/14  MDG 1.2 corrected definition of Sgh and simplified loops; resolves issue #3
!--------------------------------------------------------------------------
recursive subroutine ECCICalcSgh(cell,nn,Sgh,nat)
!DEC$ ATTRIBUTES DLLEXPORT :: ECCICalcSgh

use local
use crystal
use gvectors
use constants
use symmetry
use typedefs

IMPLICIT NONE

type(unitcell)                          :: cell
integer(kind=irg),INTENT(IN)            :: nn
! corrected on 10/22/14 in response to issue #3 of release 2.0 raised by jfikar
! Sgh array changed to one dimensional vector instead of 2D array
complex(kind=dbl),INTENT(INOUT) :: Sgh(nn)
!f2py intent(in,out) ::  Sgh
integer(kind=irg),INTENT(INOUT) :: nat(maxpasym)
!f2py intent(in,out) ::  nat

integer(kind=irg)                       :: ip, ir, ic, kkk(3), ikk, n, numset
real(kind=sgl)                          :: Znsq, DBWF, kkl
complex(kind=dbl)                       :: carg
real(kind=dbl)                          :: ctmp(192,3),arg, tpi

  tpi = 2.D0 * cPi
  Sgh = cmplx(0.D0,0.D0)
  numset = cell % ATOM_ntype  ! number of special positions in the unit cell

! comment: this can likely be further simplified and we'll take a closer look at this for
! the next release.

! for each special position we need to compute its contribution to the Sgh vector
  do ip=1,numset
    call CalcOrbit(cell,ip,n,ctmp)
    nat(ip) = n
! get Zn-squared for this special position, and include the site occupation parameter as well
    Znsq = float(cell%ATOM_type(ip))**2 * cell%ATOM_pos(ip,4)
! loop over all contributing reflections
! ir is the row index
! we don't need to actually go through this list, so the lines are commented out
!   rltmpa => reflist%next    ! point to the front of the list
!   do ir=1,nn  ! we only need a single summation loop since Sgh is a column vector
! Debye-Waller exponential times Z^2
     Sgh = Sgh + cmplx(n * Znsq, 0.D0)
!    rltmpa => rltmpa%next  ! move to next row-entry
!   end do
  end do

! this is the older code which I'm leaving here commented [10/22/14]
!! for each special position we need to compute its contribution to the Sgh array
!  do ip=1,numset
!    call CalcOrbit(ip,n,ctmp)
!    nat(ip) = n
!! get Zn-squared for this special position, and include the site occupation parameter as well
!    Znsq = float(cell%ATOM_type(ip))**2 * cell%ATOM_pos(ip,4)
!! loop over all contributing reflections
!! ir is the row index
!    rltmpa => reflist%next    ! point to the front of the list
!    do ir=1,nn
!! ic is the column index
!      rltmpb => reflist%next    ! point to the front of the list
!      do ic=1,nn
!        kkk = rltmpb%hkl - rltmpa%hkl
!! We'll assume isotropic Debye-Waller factors for now ...
!! That means we need the square of the length of s=  kk^2/4
!        kkl = 0.25 * CalcLength(float(kkk),'r')**2
!! Debye-Waller exponential times Z^2
!        DBWF = Znsq * exp(-cell%ATOM_pos(ip,5)*kkl)
!! here is where we should insert the proper weight factor, Z^2 exp[-M_{h-g}]
!! and also the detector geometry...   For now, we do nothing with the detector
!! geometry; the Rossouw et al 1994 paper lists a factor A that does not depend
!! on anything in particular, so we assume it is 1.
!        do ikk=1,n
!! get the argument of the complex exponential
!          arg = tpi*sum(kkk(1:3)*ctmp(ikk,1:3))
!          carg = cmplx(dcos(arg),dsin(arg))
!! multiply with the prefactor and add
!          Sgh(ir,ic) = Sgh(ir,ic) + carg * cmplx(DBWF,0.D0)
!        end do
!        rltmpb => rltmpb%next  ! move to next column-entry
!      end do
!     rltmpa => rltmpa%next  ! move to next row-entry
!   end do
!  end do

end subroutine ECCICalcSgh

!C***********************************************************************
!C
!C                        naninfchk.f
!C
!C      *****************************************************************
!C      *                                                               *
!C	* 	Absoft Corporation 					*
!C 	*	2781 Bond Street					*
!C	*	Rochester Hills, MI  48309				*
!C	*								*
!C	*	This file contains example code for demonstration	*
!C	*	purposes only.  Absoft makes no warranty of the	*
!C	*	suitability of this code for any purpose.		*
!C	*								*
!C	*	In no event shall Absoft be liable for any incidental,*
!C	*	indirect, special, or consequential damages arising	*
!C	*	out of the use of this code.				*
!C	*								*
!C	*****************************************************************
!C
!C Routines to test real and double values against NaN and INF
!C
!C            NANCHK(X) - tests REAL*4 value X against NaN
!!C            DNANCHK(X) - tests REAL*8 value X against NaN
!!C            INFCHK(X) - tests REAL*4 value X against INF
!!C            DINFCHK(X) - test REAL*8 value X against INF
!C
!C For little endian machines (Intel x86), compile with
!C
!C      f77 -c -DBYTE_SWAPPED=1 naninfchk.f
!C	or
!C      f90 -c -DBYTE_SWAPPED=1 naninfchk.f -YBOZTYPE=INT
!C
!C For big endian machines (PowerPC), compile with
!C
!C      f77 -c naninfchk.f
!C	or
!C      f90 -c naninfchk.f -YBOZTYPE=INT
!C
!C***********************************************************************
RECURSIVE LOGICAL FUNCTION NANCHK(X)
!DEC$ ATTRIBUTES DLLEXPORT :: NANCHK

IMPLICIT NONE
REAL,INTENT(IN)      :: X
REAL                 :: Y
INTEGER              :: I
EQUIVALENCE(Y,I)

Y = X
NANCHK = isnan(Y) !((I .AND. z'7f800000') .EQ. z'7f800000') .AND.((I .AND. z'007fffff') .NE. z'00000000')

RETURN
END

end module ECCImod

max_stars_repo_name: vmonardo/EMsoft

max_stars_repo_licenses: Intel BSD-3-Clause MIT-CMU Unlicense

max_issues_repo_name: vmonardo/EMsoft

max_issues_repo_licenses: Intel BSD-3-Clause MIT-CMU Unlicense

max_forks_repo_name: vmonardo/EMsoft

max_forks_repo_licenses: Intel BSD-3-Clause MIT-CMU Unlicense


text:
subroutine sup_SigmaMatrix(ndime,tens,sig,sigmatrix)
   !This subroutine convert the vector sigma=(sig_x, sig_y, sig_xy)
   !in a symmetric matrix.
   use Mod_Element
   use typre
   implicit none
   integer(ip) :: idime, jdime, cont,ndime,tens
   real(rp) :: sigmatrix(ndime,ndime), sig(tens)

   if (ndime==2) then
      sigmatrix(1,1)=sig(1)
      sigmatrix(2,2)=sig(2)
      sigmatrix(1,2)=sig(3)
      sigmatrix(2,1)=sig(3)
   else if (ndime==3) then
      sigmatrix(1,1)=sig(1)
      sigmatrix(2,2)=sig(2)
      sigmatrix(3,3)=sig(3)
      sigmatrix(2,3)=sig(4)
      sigmatrix(3,2)=sig(4)
      sigmatrix(1,3)=sig(5)
      sigmatrix(3,1)=sig(5)
      sigmatrix(1,2)=sig(6)
      sigmatrix(2,1)=sig(6)
   end if

end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20071107 18:06:20 nbsp  Welcome to the Wiki. Please read Welcome to the Wiki/Business Owner; it will help explain how you can make the wiki a positive experience for you without clashing with established wiki social norms; for instance deleting information that you consider negative. Users/JasonAller

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! (C) Copyright 2020 UCAR
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.

module fv3jedi_vc_vertremap_mod

! fckit
use fckit_configuration_module, only: fckit_configuration

! fms
use constants_mod,         only: grav
use field_manager_mod,     only: MODEL_ATMOS
use mpp_domains_mod,       only: mpp_update_domains
use tracer_manager_mod,    only: get_number_tracers, get_tracer_names, get_tracer_index, NO_TRACER, &
                                 set_tracer_profile

! fv3
use external_ic_mod,       only: remap_scalar, remap_dwinds, source
use fv_arrays_mod,         only: fv_atmos_type, deallocate_fv_atmos_type, R_GRID
use fv_grid_utils_mod,     only: mid_pt_sphere, get_unit_vect2, get_latlon_vector, inner_prod
use test_cases_mod,        only: checker_tracers

! fv3jedi
use fv_prec_mod,           only: kind_fv3
use fv_init_mod,           only: fv_init
use fv3jedi_geom_mod,      only: fv3jedi_geom
use fv3jedi_fieldfail_mod, only: field_fail
use fv3jedi_field_mod,     only: copy_subset, field_clen, fv3jedi_field
use fv3jedi_kinds_mod,     only: kind_real
use fv3jedi_state_mod,     only: fv3jedi_state

implicit none
private
public :: fv3jedi_vc_vertremap

type :: fv3jedi_vc_vertremap
 type(fv_atmos_type), allocatable :: Atm(:)
 logical, allocatable :: grids_on_this_pe(:)
 logical :: from_cold_start, checker_tr
 integer :: nt_checker
 contains
   procedure :: create
   procedure :: delete
   procedure :: changevar
end type fv3jedi_vc_vertremap

! --------------------------------------------------------------------------------------------------

contains

! --------------------------------------------------------------------------------------------------

subroutine create(self, geom, conf)

class(fv3jedi_vc_vertremap), intent(inout) :: self
type(fv3jedi_geom),          intent(in)    :: geom
type(fckit_configuration),   intent(in)    :: conf

integer :: gtile, p_split = 1, n
logical :: checks_passed
character(len=:), allocatable :: str

! Create Atm structure
call fv_init(self%Atm, 300.0_kind_real, self%grids_on_this_pe, p_split, gtile, .false.)

! Flag to use cold starts
if( .not. conf%get('input is cold starts', self%from_cold_start) ) self%from_cold_start = .true.

! Check tracers
if( .not. conf%get('check tracers',    self%checker_tr) ) self%checker_tr = .false.
if( .not. conf%get('check tracers nt', self%nt_checker) ) self%nt_checker = 0

! Set global source variable in external_ic
if (.not. conf%get("source of inputs", str)) then
  str = 'FV3GFS GAUSSIAN NETCDF FILE'
endif
source = trim(str)

! Remapping needs nggps_ic to be true
self%Atm(1)%flagstruct%nggps_ic = .true.

! ak/bk
self%Atm(1)%ak = real(geom%ak,kind_fv3)
self%Atm(1)%bk = real(geom%bk,kind_fv3)
self%Atm%ptop = real(geom%ak(1),kind_fv3)

! Sanity checks
checks_passed = .true.
if (checks_passed) checks_passed = geom%npx == self%Atm(1)%npx
if (checks_passed) checks_passed = geom%npy == self%Atm(1)%npy
if (checks_passed) checks_passed = geom%npz == self%Atm(1)%npz
if (checks_passed) checks_passed = geom%isd == self%Atm(1)%bd%isd
if (checks_passed) checks_passed = geom%ied == self%Atm(1)%bd%ied
if (checks_passed) checks_passed = geom%jsd == self%Atm(1)%bd%jsd
if (checks_passed) checks_passed = geom%jed == self%Atm(1)%bd%jed
if (.not.checks_passed) call abor1_ftn("fv3jedi_vc_vertremap_mod.field_fail: Geometry generated"// &
                                       " here does not match fv3-jedi geometry.")

if (.not.size(self%Atm)==1) call abor1_ftn("fv3jedi_vc_vertremap_mod.field_fail: Atm strucutre"// &
                                           " with size > 1 not supported.")

end subroutine create

! --------------------------------------------------------------------------------------------------

subroutine delete(self)

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

integer :: n

do n = 1,size(self%Atm)
  call deallocate_fv_atmos_type(self%Atm(n))
enddo
deallocate(self%Atm)
deallocate(self%grids_on_this_pe)

end subroutine delete

! --------------------------------------------------------------------------------------------------

subroutine changevar(self, xin, xout)

class(fv3jedi_vc_vertremap), target, intent(inout) :: self
type(fv3jedi_state),                 intent(in)    :: xin
type(fv3jedi_state),                 intent(inout) :: xout

! Copy versus transform control
integer :: f
character(len=field_clen), allocatable :: fields_to_do(:)
type(fv3jedi_field), pointer :: field_ptr

! Local versions of fv3jedi state
real(kind=kind_real), allocatable :: orog_filt(:,:,:)
real(kind=kind_real), allocatable :: ps_cold(:,:,:)
real(kind=kind_real), allocatable :: zh_cold(:,:,:)
real(kind=kind_real), allocatable ::  w_cold(:,:,:)
real(kind=kind_real), allocatable ::  t_cold(:,:,:)
real(kind=kind_real), allocatable :: q_tmp(:,:,:)
real(kind=kind_real), allocatable :: ud_cold(:,:,:)
real(kind=kind_real), allocatable :: vd_cold(:,:,:)

logical :: have_remapped
type(fv_atmos_type), pointer :: Atm
integer:: i, j, k, nt, ntracers, ntprog, itoa, levp, isc, iec, jsc, jec, npz, nts
integer:: liq_wat, ice_wat, rainwat, snowwat, graupel, ntclamt
character(len=64) :: tracer_name
real(kind=kind_fv3), allocatable :: ak(:), bk(:)
real(kind=kind_fv3), allocatable ::  q_cold(:,:,:,:)
real(kind=kind_fv3) :: wt, qt, m_fac

! Field names
character(len=field_clen) :: ps_fname
character(len=field_clen) :: zh_fname
character(len=field_clen) :: w_fname
character(len=field_clen) :: t_fname
character(len=field_clen) :: ud_fname
character(len=field_clen) :: vd_fname

! Identity part of the change of variables
! ----------------------------------------
call copy_subset(xin%fields, xout%fields, fields_to_do)

! If variable change is the identity early exit
! ---------------------------------------------
if (.not.allocated(fields_to_do)) return

! Field names
! -----------
ps_fname = 'ps'
zh_fname = 'zh'
w_fname = 'w'
t_fname = 't'
ud_fname = 'ud'
vd_fname = 'vd'
if (self%from_cold_start) then
  ps_fname = trim(ps_fname)//'_cold'
  zh_fname = trim(zh_fname)//'_cold'
  w_fname = trim(w_fname)//'_cold'
  t_fname = trim(t_fname)//'_cold'
  ud_fname = trim(ud_fname)//'_cold'
  vd_fname = trim(vd_fname)//'_cold'
endif

! Remap to new Lagrangian vertical coordinate
! -------------------------------------------
have_remapped = .false.
if ( xin%has_field(ps_fname) .and. xin%has_field(zh_fname) .and. &
     xin%has_field(w_fname) .and. &
     xin%has_field('orog_filt') .and. &
     xin%has_field(ud_fname) .and. xin%has_field(vd_fname) ) then

  ! Shortcuts
  Atm => self%Atm(1)
  isc = Atm%bd%is
  iec = Atm%bd%ie
  jsc = Atm%bd%js
  jec = Atm%bd%je
  npz = Atm%npz

  ! Orography
  call xin%get_field('orog_filt', orog_filt)
  Atm%phis(isc:iec,jsc:jec) = real(orog_filt(isc:iec,jsc:jec,1),kind_fv3)*grav  ! Convert to phis

  ! Dynamics fields
  call xin%get_field(ps_fname, ps_cold)
  call xin%get_field(zh_fname, zh_cold)
  call xin%get_field( w_fname,  w_cold)

  ! akbk from cold starts have different levels (for now extra levels are zero)
  levp = size(w_cold,3)
  itoa = levp - npz + 1

  allocate(ak(levp+1))
  allocate(bk(levp+1))
  ak = 0.0_kind_fv3
  bk = 0.0_kind_fv3
  ak(itoa:levp+1) = Atm%ak(1:npz+1)
  bk(itoa:levp+1) = Atm%bk(1:npz+1)
  ak(1) = max(1.e-9_kind_fv3, ak(1))

  ! Tracers
  ! -------
  ! Number of tracers in Atm
  call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers, num_prog=ntprog)
  call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers, num_prog=ntprog)

  ! initialize all tracers to default values prior to being input
  do nt = 1, ntprog
    call get_tracer_names(MODEL_ATMOS, nt, tracer_name)
    ! set all tracers to an initial profile value
    call set_tracer_profile (MODEL_ATMOS, nt, Atm%q(:,:,:,nt)  )
  enddo
  do nt = ntprog+1, ntracers
    call get_tracer_names(MODEL_ATMOS, nt, tracer_name)
    ! set all tracers to an initial profile value
    call set_tracer_profile (MODEL_ATMOS, nt, Atm%qdiag(:,:,:,nt)  )
  enddo

  ! Copy from fv3-jedi to temporary cold q structure
  allocate (q_cold(isc:iec, jsc:jec, levp, ntracers))
  q_cold = 0.0_kind_fv3
  do nt = 1, ntracers
    call get_tracer_names(MODEL_ATMOS, nt, tracer_name)
    if (self%from_cold_start) tracer_name = trim(tracer_name)//"_cold"
    if (xin%has_field(trim(tracer_name))) then
      call xin%get_field(trim(tracer_name), q_tmp)
      q_cold(:,:,:,nt) = real(q_tmp, kind_fv3)
    endif
  enddo

  ! Call remapping non-wind variables
  ! ---------------------------------
  if (xin%has_field(t_fname)) then
    call xin%get_field( t_fname,  t_cold)
    call remap_scalar(Atm, levp, npz, ntracers, ak, bk, real(ps_cold(:,:,1),kind_fv3), q_cold, &
                      real(zh_cold,kind_fv3), real(w_cold,kind_fv3), real(t_cold,kind_fv3))
  else
    call remap_scalar(Atm, levp, npz, ntracers, ak, bk, real(ps_cold(:,:,1),kind_fv3), q_cold, &
                      real(zh_cold,kind_fv3), real(w_cold,kind_fv3))
  endif

  ! Call remapping wind variables
  ! -----------------------------
  call xin%get_field(ud_fname, ud_cold)
  call xin%get_field(vd_fname, vd_cold)

  call remap_dwinds(levp, npz, ak, bk, real(ps_cold(:,:,1),kind_fv3), real(ud_cold,kind_fv3), &
                    real(vd_cold,kind_fv3), Atm)

  ! Tracer weighting
  ! ----------------
  liq_wat = get_tracer_index(MODEL_ATMOS, 'liq_wat')
  ice_wat = get_tracer_index(MODEL_ATMOS, 'ice_wat')
  rainwat = get_tracer_index(MODEL_ATMOS, 'rainwat')
  snowwat = get_tracer_index(MODEL_ATMOS, 'snowwat')
  graupel = get_tracer_index(MODEL_ATMOS, 'graupel')
  ntclamt = get_tracer_index(MODEL_ATMOS, 'cld_amt')

  if (self%from_cold_start) then
    do k = 1,npz
      do j = jsc,jec
        do i = isc,iec
          wt = Atm%delp(i,j,k)
          if ( Atm%flagstruct%nwat == 6 ) then
            qt = wt*(1.0_kind_fv3 + Atm%q(i,j,k,liq_wat) + Atm%q(i,j,k,ice_wat) + &
                                    Atm%q(i,j,k,rainwat) + Atm%q(i,j,k,snowwat) + &
                                    Atm%q(i,j,k,graupel))
          else
            qt = wt*(1.0_kind_fv3 + sum(Atm%q(i,j,k,2:Atm%flagstruct%nwat)))
          endif
          Atm%delp(i,j,k) = qt
          if (ntclamt > 0) Atm%q(i,j,k,ntclamt) = 0.0
        enddo
      enddo
    enddo
  else
    ! TODO Question is do we need this if just remapping after adding the increment, say?
    do k = 1,npz
      do j = jsc,jec
        do i = isc,iec
          wt = Atm%delp(i,j,k)
          if ( Atm%flagstruct%nwat == 6 ) then
            qt = wt*(1.0_kind_fv3 + Atm%q(i,j,k,liq_wat) + Atm%q(i,j,k,ice_wat) + &
                                    Atm%q(i,j,k,rainwat) + Atm%q(i,j,k,snowwat) + &
                                    Atm%q(i,j,k,graupel))
          else
             qt = wt*(1.0_kind_fv3 + sum(Atm%q(i,j,k,2:Atm%flagstruct%nwat)))
          endif
          m_fac = wt / qt
          do nt=1,ntracers
            Atm%q(i,j,k,nt) = m_fac * Atm%q(i,j,k,nt)
          enddo
          Atm%delp(i,j,k) = qt
          if (ntclamt > 0) Atm%q(i,j,k,ntclamt) = 0.0
        enddo
      enddo
    enddo
  endif

  if (self%checker_tr) then
    nts = ntracers - self%nt_checker+1
    call checker_tracers(isc, iec, jsc, jec, Atm%bd%isd, Atm%bd%ied, Atm%bd%jsd, Atm%bd%jed, &
                         self%nt_checker, npz, Atm%q(:,:,:,nts:ntracers), &
                         Atm%gridstruct%agrid_64(isc:iec,jsc:jec,1),     &
                         Atm%gridstruct%agrid_64(isc:iec,jsc:jec,2), &
                         real(9.0_kind_real,kind_fv3), real(9.0_kind_real,kind_fv3))
  endif

  have_remapped = .true.

endif

! Loop over the fields not found in the input state and work through cases
! ------------------------------------------------------------------------
do f = 1, size(fields_to_do)

  ! Remapping needs to have happened in order to copy anything
  if (.not. have_remapped) call field_fail(fields_to_do(f))

  call xout%get_field(trim(fields_to_do(f)), field_ptr)

  select case (trim(fields_to_do(f)))

  case ("ud")

    field_ptr%array(isc:iec,jsc:jec+1,:) = Atm%u(isc:iec,jsc:jec+1,:)

  case ("vd")

    field_ptr%array(isc:iec+1,jsc:jec,:) = Atm%v(isc:iec+1,jsc:jec,:)

  case ("t")

    field_ptr%array(isc:iec,jsc:jec,:) = Atm%pt(isc:iec,jsc:jec,:)

  case ("delp")

    field_ptr%array(isc:iec,jsc:jec,:) = Atm%delp(isc:iec,jsc:jec,:)

  case ("ps")

    field_ptr%array(isc:iec,jsc:jec,1) = Atm%ps(isc:iec,jsc:jec)

  case ("delz")

    field_ptr%array(isc:iec,jsc:jec,:) = Atm%delz(isc:iec,jsc:jec,:)

  case ("w")

    field_ptr%array(isc:iec,jsc:jec,:) = Atm%w(isc:iec,jsc:jec,:)

  case ("phis")

    field_ptr%array(isc:iec,jsc:jec,1) = Atm%phis(isc:iec,jsc:jec)

  case ('sphum','liq_wat','ice_wat','o3mr','graupel','snowwat','rainwat')

    nt = get_tracer_index(MODEL_ATMOS, trim(fields_to_do(f)))
    if (nt == NO_TRACER) call field_fail(fields_to_do(f))
    field_ptr%array(isc:iec,jsc:jec,:) = Atm%q(isc:iec,jsc:jec,:,nt)

  case ('sgs_tke','cld_amt')

    field_ptr%array = 0.0_kind_real

  case default

    call abor1_ftn("fv3jedi_vc_coldstartwinds_mod.changevar unknown field: "//trim(fields_to_do(f))&
                   //". Not in input field and no transform case specified.")

  end select

enddo

! Copy calendar infomation
! ------------------------
xout%calendar_type = xin%calendar_type
xout%date_init = xin%date_init

end subroutine changevar

! --------------------------------------------------------------------------------------------------

end module fv3jedi_vc_vertremap_mod


text:
!> This unit tests the gs2 main interface,
!! and by repeatedly initializing and finalizing
!! gs2, tests that gs2 is being properly tidied up,
!! variables deallocated etc.
!
program test_gs2_main
  use gs2_main
  use unit_tests
  use mp, only: init_mp, finish_mp, mp_comm
  implicit none
  real :: eps
  type(gs2_program_state_type) :: state

  eps = 1.0e-7
  if (precision(eps).lt. 11) eps = eps * 1000.0

  call init_mp

  call announce_module_test("gs2_main")

  state%mp_comm_external = .true.
  state%mp_comm = mp_comm

  call initialize_gs2(state)
  call finalize_gs2(state)

  call initialize_gs2(state)
  call finalize_gs2(state)

  call initialize_gs2(state)
  call initialize_equations(state)
  call finalize_equations(state)
  call finalize_gs2(state)

  call initialize_gs2(state)
  call initialize_equations(state)
  call finalize_equations(state)
  call finalize_gs2(state)

  !!program gs2
    !!type(gs2_program_state_type) :: state
    call initialize_gs2(state)
    call initialize_equations(state)
    call initialize_diagnostics(state)
    !if (state%eigsolve) then
      !call solve_eigenproblem
    !else
    call evolve_equations(state, state%nstep/2)
    call evolve_equations(state, state%nstep/2)
    !! This call should do nothing and print a warning
    call evolve_equations(state, state%nstep/2)
    call evolve_equations(state, state%nstep/2)

    call calculate_outputs(state)

    call finalize_diagnostics(state)
    call finalize_equations(state)
    call finalize_gs2(state)
  !!end program gs2

    !!call init_mp

    !!test_driver_flag = .true.
    !!functional_test_flag = .true.

  call finalize_overrides(state)
  call close_module_test("gs2_main")

  !call finish_gs2

  call finish_mp

contains

end program test_gs2_main

ext: fpp


text:
C
C $Id: icloem.f,v 1.5 2008-07-27 00:17:30 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      FUNCTION ICLOEM (MESSG)
C
        CHARACTER*(*) MESSG
C
C ICLOEM(MESSG) is the index of the last non-blank character in MESSG.
C
        DO 101 I=LEN(MESSG),1,-1
         IF (MESSG(I:I).NE.' ') THEN
           ICLOEM=I
           RETURN
         END IF
  101   CONTINUE
C
        ICLOEM=1
C
        RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


      subroutine water_allocation_read

      use input_file_module
      use water_allocation_module
      use mgt_operations_module
      use maximum_data_module

      implicit none

      character (len=80) :: titldum   !           |title of file
      character (len=80) :: header    !           |header of file
      integer :: eof                  !           |end of file
      integer :: imax                 !none       |determine max number for array (imax) and total number in file
      logical :: i_exist              !none       |check to determine if file exists
      integer :: i                    !none       |counter
      integer :: k                    !none       |counter
      integer :: isrc                 !none       |counter
      integer :: iwro                 !none       |counter
      integer :: num_objs
      integer :: idmd
      integer :: idb

      eof = 0
      imax = 0

      !! read water allocation inputs
      inquire (file=in_watrts%water_rights, exist=i_exist)
      if (.not. i_exist .or. in_watrts%water_rights == "null") then
        allocate (wallo(0:0))
      else
      do
        open (107,file=in_watrts%water_rights)
        read (107,*,iostat=eof) titldum
        if (eof < 0) exit
        read (107,*,iostat=eof) imax
        if (eof < 0) exit

        allocate (wallo(imax))

        do iwro = 1, imax
          read (107,*,iostat=eof) header
          if (eof < 0) exit
          read (107,*,iostat=eof) wallo(iwro)%name, wallo(iwro)%rule_typ, wallo(iwro)%res_lim,         &
            wallo(iwro)%comp, wallo(iwro)%dmd_obs
          if (eof < 0) exit
          read (107,*,iostat=eof) header
          if (eof < 0) exit
          num_objs = wallo(iwro)%dmd_obs
          db_mx%wallo_db = num_objs
          allocate (wallo(iwro)%dmd(num_objs))
          allocate (wallod_out(num_objs))
          allocate (wallom_out(num_objs))
          allocate (walloy_out(num_objs))
          allocate (walloa_out(num_objs))

          do idmd = 1, num_objs
            read (107,*,iostat=eof) i
            if (eof < 0) exit
            backspace (107)
            read (107,*,iostat=eof) k, wallo(iwro)%dmd(i)%ob_typ, wallo(iwro)%dmd(i)%obtyp_num,                 &
              wallo(iwro)%dmd(i)%dmd_typ, wallo(iwro)%dmd(i)%amount, (wallo(iwro)%dmd(i)%src(isrc), isrc = 1, 2)

            !! xwalk with irrigation database
            if (wallo(iwro)%dmd(i)%ob_typ == "hru") then
            do idb = 1, db_mx%irrop_db
              if (wallo(iwro)%dmd(i)%dmd_typ == irrop_db(idb)%name) then
                wallo(iwro)%dmd(i)%irr_typ = idb
                exit
              end if
            end do
            end if

            !! zero output variables for summing
            wallod_out(idmd)%src(1) = walloz
            wallod_out(idmd)%src(2) = walloz
            wallom_out(idmd)%src(1) = walloz
            wallom_out(idmd)%src(2) = walloz
            walloy_out(idmd)%src(1) = walloz
            walloy_out(idmd)%src(2) = walloz
            walloa_out(idmd)%src(1) = walloz
            walloa_out(idmd)%src(2) = walloz

          end do
        end do

        exit
      enddo
      endif
      close(107)

      return
      end subroutine water_allocation_read


text:
C=======================================================================
C  SECROP, Subroutine
C
C  Determines crop sensitivity analysis options
C-----------------------------------------------------------------------
C  Revision history
C
C  01/01/1990 GH  Written
C  05/28/1993 PWW Header revision and minor change
C  05/01/2003 GH  Added citrus
C-----------------------------------------------------------------------
C  INPUT  : FILEG,IDETA,CROP,CROPD,PATHGE
C
C  LOCAL  : IDETA,BLANK,CROPC,ERRKEY,MENU,NLOOP
C
C  OUTPUT :
C-----------------------------------------------------------------------
C  Called : SENS
C
C  Calls  : ERROR CLEAR SWCROP IDCROP
C-----------------------------------------------------------------------
C                         DEFINITIONS
C
C  HDLAY  :
C=======================================================================

      SUBROUTINE SECROP (FILEC,FILEE,FILEG,RNMODE,CROP,CROPD,PATHCR)

      IMPLICIT NONE

      CHARACTER*1  RNMODE
      CHARACTER*2  CROP,CROPC
      CHARACTER*6  ERRKEY
      CHARACTER*16 CROPD
      CHARACTER*12 FILEC,FILEE,FILEG
      CHARACTER*80 PATHCR

      INTEGER      MENU,NLOOP

      PARAMETER (ERRKEY = 'SECROP')

      CROPC = CROP
      NLOOP = 0

  100 CONTINUE
      NLOOP = NLOOP + 1
      IF (NLOOP .GT. 25) CALL ERROR (ERRKEY,1,' ',0)
      IF (INDEX('IE',RNMODE).GT. 0)  THEN
         CALL CLEAR
         WRITE (*,200) CROPD,FILEC
         IF (FILEC(8:8) .NE. '0') WRITE (*,250) FILEC
         IF (CROPC .NE. CROP) WRITE (*,260) CROPD
         WRITE (*,275)
      ENDIF
C
C     Get menu choice
C
      READ (5,'(I2)',ERR = 100) MENU
C
C     Branch to menu choice
C
      IF (MENU .EQ. 0) THEN
          RETURN
      ELSE IF (MENU .EQ. 1) THEN
          CALL SWCROP (FILEC,FILEE,FILEG,RNMODE,CROP,CROPD,PATHCR)
      ELSE IF (MENU .EQ. 2) THEN
          CALL IDCROP (FILEC,RNMODE,PATHCR)
      ENDIF

      GO TO 100

C-----------------------------------------------------------------------
C     Format Strings
C-----------------------------------------------------------------------

  200 FORMAT (//,
     &       15X,'CROP SELECTION AND MODIFICATION',/,
     &       15X,'===============================',//,
     & 5X,' 0. Return to Main Menu ',//,
     & 5X,' 1. Crop Selection ....................[ ',A16,/,
     & 5X,' 2. Crop Parameter File ...............[ ',A12,/)
  250 FORMAT (9X,
     & 'A NON-STANDARD crop parameter file is selected [',A12,']',/,9X,
     & 'Simulated results can potentially be different.')
  260 FORMAT (9X,
     & 'A new crop ',A16,' is selected !!',/,9X,
     & 'Please reselect cultivar and other management options ',
     & 'in the main menu.')
  275 FORMAT (//,9X,'SELECTION ? [ Default=0 ] ===> ',$)

      END SUBROUTINE SECROP

C=======================================================================
C  SWCROP, Subroutine
C
C  Determines crop sensitivity analysis options
C-----------------------------------------------------------------------
C  Revision history
C
C  01/01/1991 GH  Written
C  05/28/1993 PWW Header revision and minor changes
C-----------------------------------------------------------------------
C  INPUT  : FILEG,IDETA,CROP,CROPD,PATHGE
C
C  LOCAL  : LINE,BLANK,CROPID,ERRKEY,CROPS,CRPTEM,FILECC,I,NLOOP,PATHL,
C           FEXIST,CID,FLAG
C
C  OUTPUT :
C-----------------------------------------------------------------------
C  Called : SECROP
C
C  Calls  : CLEAR ERROR VERIFY
C-----------------------------------------------------------------------
C                         DEFINITIONS
C
C  HDLAY  :
C=======================================================================

      SUBROUTINE SWCROP (FILEC,FILEE,FILEG,RNMODE,CROP,CROPD,PATHCR)

      IMPLICIT NONE

      CHARACTER*1  RNMODE,LINE(80),BLANK
      CHARACTER*2  CROP,CROPID(50)
      CHARACTER*6  ERRKEY
      CHARACTER*16 CROPD
      CHARACTER*12 FILEC,CRPTEM,FILEE,FILEG
      CHARACTER*80 PATHCR
      CHARACTER*92 FILECC

      INTEGER      I,NLOOP,PATHL
      LOGICAL      FEXIST
      REAL         CID,FLAG

      PARAMETER (ERRKEY = 'SWCROP')
      PARAMETER (BLANK  = ' ')

      DATA CROPID/'AR','BA','BH','BM','BN','BR','C3','C4','CB','CH',
     &            'CO','CP','CS','CT','FA','FB','G1','G2','G3','G4',
     &            'G5','G6','G7','G8','GB','LT','ML','MZ','NP','PE',
     &            'PI','PN','PP','PR','PT','RI','SB','SC','SG','SU',
     &            'TM','TN','TR','VB','WH','  ','  ','  ','  ','  '/

      CRPTEM = FILEC
      PATHL  = INDEX (PATHCR,BLANK)
      WRITE (CRPTEM(8:8),'(A1)') '5'

      NLOOP  = 0

      IF (INDEX('IE',RNMODE) .GT. 0) THEN
         CALL CLEAR
         WRITE (*,400)
         DO I = 1, 35
            WRITE (CRPTEM(1:2),600) CROPID(I)
            IF (PATHL .LE. 1) THEN
               FILECC = CRPTEM
             ELSE
               FILECC = PATHCR(1:(PATHL-1)) // CRPTEM
            ENDIF
	      INQUIRE (FILE = FILECC, EXIST = FEXIST)
            IF (FEXIST) THEN
	         CROP  = CROPID(I)
	         CALL GET_CROPD(CROP, CROPD)
               WRITE (*,700) I,CROPD
            ENDIF
         END DO
      ENDIF

 1020 CONTINUE

      NLOOP = NLOOP + 1
      IF (NLOOP .GT. 25) CALL ERROR (ERRKEY,1,FILEC,0)
      LINE(1) = ' '
      IF (INDEX('IE',RNMODE) .GT. 0) WRITE (*,900)
      READ (5,1000) LINE
      CALL VERIFY (LINE,CID,FLAG)
      IF (CID .LE. 0.0) THEN
         IF (PATHL .LE. 1) THEN
            FILECC = FILEC
          ELSE
            FILECC = PATHCR(1:(PATHL-1)) // FILEC
         ENDIF
         INQUIRE (FILE = FILECC,EXIST = FEXIST)
         IF (.NOT. FEXIST) THEN
            WRITE (*,1040)
            GO TO 1020
         ENDIF
       ELSE
         IF (FLAG .GT. 0 .OR. CID .GT. 35) THEN
            WRITE (*,1010)
            GO TO 1020
          ELSE
            WRITE (CRPTEM(1:2),'(A2)') CROPID(NINT(CID))
            IF (PATHL .LE. 1) THEN
               FILECC = CRPTEM
             ELSE
               FILECC = PATHCR(1:(PATHL-1)) // CRPTEM
            ENDIF
            INQUIRE (FILE = FILECC,EXIST = FEXIST)
            IF (.NOT. FEXIST) THEN
               WRITE (*,1040)
               GO TO 1020
            ENDIF
         ENDIF

         CROP  = CROPID (NINT(CID))
         CALL GET_CROPD(CROP, CROPD)

         IF (FILEC .NE. CRPTEM) THEN
            WRITE (FILEG(1:2),600) CROP
            WRITE (FILEE(1:2),600) CROP
         ENDIF
         FILEC = CRPTEM
      ENDIF

      RETURN

C-----------------------------------------------------------------------
C     Format Strings
C-----------------------------------------------------------------------

  400 FORMAT (/,30X,'DSSAT V4.5 CROPS',/,30X,16('='))
  600 FORMAT (A2)
  700 FORMAT (30X,I2,')',5X,A16)
  900 FORMAT (/,30X,'SELECT CROP (#) ====> ',$)
 1000 FORMAT (80A1)
 1010 FORMAT (10X,'Please enter a number between 1 and 20 ! ')
 1040 FORMAT (/,10X,'You selected a crop which does not exist.',/,
     &          10X,'Please select another crop !')

      END SUBROUTINE SWCROP

C=======================================================================
C  IDCROP, Subroutine
C
C  Determines selection of appropriate crop parameter file
C-----------------------------------------------------------------------
C  Revision history
C
C  01/01/1991 GH  Written
C  05/28/1993 PWW Header revision and minor changes
C-----------------------------------------------------------------------
C  INPUT  : FILEC,PATHCR
C
C  LOCAL  : LINE,BLANK,CROPFL,FILECC,I,PATHL,FEXIST,CPID,FLAG
C
C  OUTPUT :
C-----------------------------------------------------------------------
C  Called : SECROP
C
C  Calls  : CLEAR VERIFY
C-----------------------------------------------------------------------
C                         DEFINITIONS
C
C  HDLAY  :
C=======================================================================

      SUBROUTINE IDCROP (FILEC,RNMODE,PATHCR)

      IMPLICIT     NONE

      CHARACTER*1  LINE(80),RNMODE,BLANK
      CHARACTER*12 FILEC,CROPFL
      CHARACTER*80 PATHCR
      CHARACTER*92 FILECC

      INTEGER      I,PATHL
      LOGICAL      FEXIST
      REAL         CPID,FLAG

      PARAMETER (BLANK = ' ')

      CROPFL = FILEC
      PATHL  = INDEX (PATHCR,BLANK)

      IF (INDEX('IE',RNMODE) .GT. 0) THEN
         CALL CLEAR
         WRITE (*,100)
      ENDIF

      DO I = 1, 10
         WRITE (CROPFL(8:8),'(I1)') I-1
         IF (PATHL .LE. 1) THEN
            FILECC = CROPFL
          ELSE
            FILECC = PATHCR(1:(PATHL-1)) // CROPFL
         ENDIF
         INQUIRE (FILE = FILECC, EXIST = FEXIST)
         IF (FEXIST) THEN
            IF (INDEX('IE',RNMODE) .GT. 0) WRITE (*,30) I-1,CROPFL
         ENDIF
      END DO

 15   CONTINUE
      LINE(1) = ' '
      IF (INDEX('IE',RNMODE) .GT. 0) WRITE (*,50) FILEC(8:8)

      READ (5,'(80A1)') LINE
      CALL VERIFY (LINE,CPID,FLAG)

      IF (CPID .LE. 0) THEN
          FILECC = FILEC
       ELSE IF (FLAG .GT. 0 .OR. CPID .GT. 20.0) THEN
         WRITE (*,70)
         GO TO 15
       ELSE
         WRITE (FILEC(8:8),'(I1)') NINT(CPID)
         IF (PATHL .LE. 1) THEN
            FILECC = FILEC
          ELSE
            FILECC = PATHCR (1:(PATHL-1)) // FILEC
         ENDIF
         INQUIRE (FILE = FILECC,EXIST = FEXIST)
         IF (.NOT. FEXIST) THEN
            WRITE (*,80)
            GO TO 15
         ENDIF
      ENDIF

      RETURN

C-----------------------------------------------------------------------
C     Format Strings
C-----------------------------------------------------------------------

   30 FORMAT (2X,I2,')',5X,A12)
   50 FORMAT (/,10X,'CROP PARAMETER FILE SELECTED        ===>',1X,A1,
     &        /,10X,'ALTERNATE SELECTION ? [Default = 0] ---> ',$)
   70 FORMAT (10X,'Please enter a number between 0 and 20! ')
   80 FORMAT (10X,'You selected a file which does not exist.',/,
     &        10X,'Please select another file !')
  100 FORMAT (10X,'CROP PARAMETER FILES AVAILABLE',/,10X,30('-'))

      END SUBROUTINE IDCROP

ext: for

max_stars_repo_name: vstocca/dssat-csm-os

max_issues_repo_name: vstocca/dssat-csm-os

max_forks_repo_name: vstocca/dssat-csm-os


!=================================================
! SUBROUTINE
!	MNEWT(x1, y1, z1, x2, y2, z2, x3, y3, z3, di, da2l)
!
!=================================================
module mod_mnewt

	implicit none

contains
!==========
subroutine MNEWT(ls, x1,y1,z1, x2,y2,z2, x3,y3,z3, di,da2l)

	use prm_matrix
	use prm_inv
	use mod_wbb

	implicit none

	integer, intent(in)  :: ls
	real(8), intent(in)  :: x1, y1, z1, x2, y2, z2, x3, di, da2l
	real(8), intent(out) :: y3, z3

	real(8) :: ct, cx, cy, cz, c2, d2, cd
	real(8) :: dx, dl, dr, dz, du, dy
	real(8) :: arr, arl, ari

	integer :: ix, iy, idb(16), j
	real(8) :: px, py, bsv(16)

	ct = 1.d-6
	cx = x2 - x1
	cy = y2 - y1
	cz = z2 - z1

	dx  = x3 - x1
	dl  = di
	dr  = 0.d0
	arr = 0.d0

	do

		y3 = y1 + dl

		! == INTDEP(x3,y3,z3) ==
		CALL WHERE1(x3, ix, px, Xa(ls), Xb(ls), Ju0(ls))
		CALL WHERE1(y3, iy, py, Ya(ls), Yb(ls), Jv0(ls))

		CALL BSVALI(px, py, bsv, Ndgf(ls), Nod2(ls))
		CALL BSPARI(ix, iy, idb, Ndgf(ls), Ju1(ls), Jv1(ls))

		z3 = 0.d0
		do j = 1, Nod2(ls)
			z3 = z3 + bsv(j) * Xcof(idb(j))
		end do
		! ====

		dz = z3 - z1

		! 二つのベクトルの面積計算 (cx,cy,cz) と (dx,dy,dz)
		! AREA2(cx, cy, cz, dx, dl, dz, c2, arl)
		dy = dl
		c2  = cx*cx + cy*cy + cz*cz
		d2  = dx*dx + dy*dy + dz*dz
		cd  = cx*dx + cy*dy + cz*dz
		arl = dsqrt(d2*c2 - cd*cd)
		! ====

		if(dabs(arl-da2l) < ct) return
		if(arl < da2l) then
			dr = dl
			arr = arl
			dl = dl * da2l / arl
			cycle
		end if

		do
			du = dr + (dl - dr) * (da2l - arr) / (arl - arr)
			y3 = y1 + du

			! == INTDEP(x3,y3,z3) ==
			CALL WHERE1(x3, ix, px, Xa(ls), Xb(ls), Ju0(ls))
			CALL WHERE1(y3, iy, py, Ya(ls), Yb(ls), Jv0(ls))

			CALL BSVALI(px, py, bsv, Ndgf(ls), Nod2(ls))
			CALL BSPARI(ix, iy, idb, Ndgf(ls), Ju1(ls), Jv1(ls))

			z3 = 0.d0
			do j = 1, Nod2(ls)
				z3 = z3 + bsv(j) * Xcof(idb(j))
			end do
			! ====

			dz = z3 - z1

			! 二つのベクトルの面積計算 (cx,cy,cz) と (dx,dy,dz)
			! AREA2(cx, cy, cz, dx, du, dz, c2, ari)
			dy = du
			c2  = cx*cx + cy*cy + cz*cz
			d2  = dx*dx + dy*dy + dz*dz
			cd  = cx*dx + cy*dy + cz*dz
			ari = dsqrt(d2*c2 - cd*cd)
			! ====

			if(dabs(ari-da2l) < ct) return

			if(ari > da2l) then
				arl = ari
				dl = du
			else if(ari < da2l) then
				arr = ari
				dr = du
			end if

		end do

	end do

end subroutine MNEWT

end module mod_mnewt


text: #redirect UC Davis


text:
Certain items belonging to UC Davis have property tags, which consist of a tag number and a barcode. Some of the older tags are pieces of metal, glued onto the item in question. These may lack a barcode if theyre old enough. The modern tags are just normal stickers. Some departments may use barcode scanners to do their inventory, but Users/WilliamLewis Ive never heard of it happening. Some departments, like the fire department, barcode specific kinds of equipment, like fire extinguishers, with a separate system.

http://accounting.ucdavis.edu/CA/EquipDef.cfm Items that must have property tags are
  Nonexpendible, having a expected life of a year or greater.
  Tangible, able to be appraised.
  Freestanding
  Possessing an acquisition cost of $5000 or more.

An item not meeting these specifications can still be tagged, at the option of the department owning the item. Once an item is tagged, it remains on inventory until the university disposes of it, even if its value drops below $500. That number, by the way, was lower in the past. So, you can have items with an acquisition cost of $500 on inventory long after it has depreciated to $0.

The tag number on a property tag tells you a bit about the item.  The first two numbers of the tag tell you the year the equipment was acquired if the year was before 2000. If the equipment was acquired in the year 2000 or later, the first four digits of the tag number tell you the year the equipment was acquired. The next three digits say something, but I dont know what. The code 300 is used on lots of things. However, you will occasionally see letters here, like CTS (classroom technology services) or MW (Mediaworks). Its hypothesized that this code indicates who the equipment in question belongs to, but 300 seems to be used everywhere.

 Department specific tags

In the case of the UC Davis Fire Department Fire Department, red bar code property tags have been used since the Summer of 2009 to keep track of all the 6500 fire extinguishers, and their proper locations, on the UC Davis campus. The numbers on each property tag represent an individual extinguisher or location in the computer database.  A tags, which are on the fire extinguishers themselves, are for each individual extinguisher. L tags, which are stuck onto the wall, hook, or cabinet where the extinguisher lives, represent the registered location of the fire extinguisher in a building. A tags and L tags often have the same numbers in a given building but as fire extinguishers are exchanged for regular maintenance and recharging this will change. The new tags were painstakingly placed on nearly all 6500 campus extinguishers over the summer of 2009 by two hard working Student Resident Firefighters.

The lucky people who do fire extinguisher inventory for the Fire Department get to use barcode scanners to make their lives easier.

All other FD items including vehicles are tagged with standard UCD property tags.

20091201 12:21:06 nbsp The presence of a UC property tag does not, in itself, indicate that the item is currently UC property.  Ive purchased several items from the Bargain Barn that still had property tags attached. Users/JimStewart

20091201 13:13:38 nbsp 300 is also the identifying tag for fire department stuff (other than fire extinguishers). Perhaps the identifier 300 more generically represents office of Admin stuff. Users/ARWENNHOLD
 I may be mistaken... thats just what my boss told me. Hmm.... it would be nice to figure this out. Ive also seen 300 on printers in the college of engineering, but figured they were just taken care of IET. Thanks for bringing this to my attention. Please edit if you learn anything new. ) wl

20091201 15:43:56 nbsp yawn Users/CovertProfessor


text:
#
# CDDL HEADER START
#
# The contents of this file are subject to the terms of the
# Common Development and Distribution License (the "License").
# You may not use this file except in compliance with the License.
#
# You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
# or http://www.opensolaris.org/os/licensing.
# See the License for the specific language governing permissions
# and limitations under the License.
#
# When distributing Covered Code, include this CDDL HEADER in each
# file and include the License file at usr/src/OPENSOLARIS.LICENSE.
# If applicable, add the following below this CDDL HEADER, with the
# fields enclosed by brackets "[]" replaced with your own identifying
# information: Portions Copyright [yyyy] [name of copyright owner]
#
# CDDL HEADER END
#
#
# Copyright 2008 Sun Microsystems, Inc.  All rights reserved.
# Use is subject to license terms.
#
# ident	"%Z%%M%	%I%	%E% SMI"

# Single-threaded writes to initially empty file.
# I/O size is set to 8KB. After every 1024 writes
# (i.e., 8MB written) fsync is called.
# The run finishes after 1GB is fully written.

set $dir=/mnt/pmem_emul
set $iosize=16k
set $filesize=4g
set $writeiters=1
set $fsynccount=128
set $iters=1
set $nthreads=1

set mode quit firstdone

define file name=bigfile,path=$dir,size=$filesize,prealloc

define process name=filewriter,instances=1
{
  thread name=filewriterthread,memsize=10m,instances=1
  {
    flowop write name=write-file,filename=bigfile,iosize=$iosize,iters=$iters
    flowop fsync name=sync-file
    flowop finishoncount name=finish,value=1000000,target=sync-file
  }
}

echo  "FileMicro-WriteFsync Version 2.1 personality successfully loaded"

run

max_stars_repo_name: ShawnZhong/SplitFS

max_issues_repo_name: ShawnZhong/SplitFS

max_forks_repo_name: ShawnZhong/SplitFS


text:
interface
   subroutine foobar() bind(c, name="foobar")
   end subroutine
end interface

end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! H0 X
! H0 X   libAtoms+QUIP: atomistic simulation library
! H0 X
! H0 X   Portions of this code were written by
! H0 X     Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
! H0 X     Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
! H0 X
! H0 X   Copyright 2006-2010.
! H0 X
! H0 X   These portions of the source code are released under the GNU General
! H0 X   Public License, version 2, http://www.gnu.org/copyleft/gpl.html
! H0 X
! H0 X   If you would like to license the source code under different terms,
! H0 X   please contact Gabor Csanyi, gabor@csanyi.net
! H0 X
! H0 X   Portions of this code were written by Noam Bernstein as part of
! H0 X   his employment for the U.S. Government, and are not subject
! H0 X   to copyright in the USA.
! H0 X
! H0 X
! H0 X   When using this software, please cite the following reference:
! H0 X
! H0 X   http://www.libatoms.org
! H0 X
! H0 X  Additional contributions by
! H0 X    Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
! H0 X
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

!X
!X  ApproxFermi module
!X
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
module ApproxFermi_module

use System_module
use linearalgebra_module ! for .feq.
use Units_module
use Functions_module

implicit none

private

public :: ApproxFermi
type ApproxFermi
  integer :: n_poles = 0
  complex(dp), allocatable :: a(:), z(:)
  real(dp) :: Fermi_E, band_width
end type ApproxFermi

public :: Initialise
interface Initialise
  module procedure ApproxFermi_Initialise, ApproxFermi_Initialise_lowlevel
end interface Initialise

public :: Finalise
interface Finalise
  module procedure ApproxFermi_Finalise
end interface Finalise

public :: Wipe
interface Wipe
  module procedure ApproxFermi_Wipe
end interface Wipe

public :: Print
interface Print
  module procedure ApproxFermi_Print
end interface Print

public :: approx_f_fermi, approx_f_fermi_deriv

contains

subroutine ApproxFermi_Finalise(this)
  type(ApproxFermi), intent(inout) :: this

  call Wipe(this)

end subroutine

subroutine ApproxFermi_Wipe(this)
  type(ApproxFermi), intent(inout) :: this

  if (allocated(this%a)) deallocate(this%a)
  if (allocated(this%z)) deallocate(this%z)
  this%n_poles = 0

end subroutine

subroutine ApproxFermi_Print(this,out)
  type(ApproxFermi),    intent(in)           :: this
  type(Inoutput), intent(inout),optional:: out

  integer::i

  call Print('ApproxFermi: Fermi_E ' // this%Fermi_E // ' band_width ' // this%band_width // ' n_poles ' // this%n_poles, file=out)

  if (this%n_poles > 0) then
    do i=1, size(this%a)
      call Print ('GreensFunc : a z ' //  i // ' ' // this%a(i) // ' ' // this%z(i), PRINT_VERBOSE, out)
    end do
  endif

end subroutine ApproxFermi_Print

subroutine ApproxFermi_Initialise(this, Fermi_E, Fermi_T, band_width)
  type(ApproxFermi), intent(inout) :: this
  real(dp), intent(in) :: Fermi_E, Fermi_T, band_width

  call Finalise(this)

  call calc_poles_from_T(this, Fermi_E, Fermi_T, band_width)

end subroutine

subroutine ApproxFermi_Initialise_lowlevel(this, n_poles, Fermi_E, band_width)
  type(ApproxFermi), intent(inout) :: this
  integer, intent(in) :: n_poles
  real(dp), intent(in) :: Fermi_E, band_width

  call calc_poles(this, Fermi_E, n_poles, band_width)

end subroutine

subroutine calc_poles(this, Fermi_E, n_poles, band_width)
  type(ApproxFermi), intent(inout) :: this
  real(dp), intent(in) :: Fermi_E, band_width
  integer, intent(in) :: n_poles

  if (n_poles == 0 .or. (band_width .feq. 0.0_dp)) then
    call system_abort("Can't calc_poles with n_poles or band_width = 0")
  endif

  this%Fermi_E = Fermi_E
  this%n_poles = n_poles
  this%band_width = band_width

  if (allocated(this%a)) then
    if (size(this%a) /= this%n_poles) deallocate(this%a)
  endif
  if (allocated(this%z)) then
    if (size(this%z) /= this%n_poles) deallocate(this%z)
  endif

  if (.not.allocated(this%z)) allocate(this%z(this%n_poles))
  if (.not.allocated(this%a)) allocate(this%a(this%n_poles))

  call calc_pole_values(this)

end subroutine calc_poles

subroutine calc_poles_from_T(this, Fermi_E, Fermi_T, band_width)
  type(ApproxFermi), intent(inout) :: this
  real(dp), intent(in) :: Fermi_E, Fermi_T, band_width

  if (band_width .feq. 0.0_dp) then
    call system_abort("Can't calc_poles_from_T with band_width = 0")
  endif

  this%n_poles = guess_n_poles(Fermi_T, band_width, this%band_width)

  call calc_poles(this, Fermi_E, this%n_poles, this%band_width)

end subroutine calc_poles_from_T

function guess_n_poles(Fermi_T, band_width_in, band_width)
  real(dp), intent(in) :: Fermi_T, band_width_in
  real(dp), intent(out) :: band_width
  integer :: guess_n_poles

  type(ApproxFermi) t_AF
  real(dp) :: Fermi_slope, AF_slope
  integer np

  Fermi_slope = f_Fermi_deriv(0.0_dp, Fermi_T, 0.0_dp)

  do np=4, 1000, 4
    call Initialise(t_AF, np, 0.0_dp, band_width_in)
    AF_slope = approx_f_Fermi_deriv(t_AF, 0.0_dp)
    if (AF_slope < Fermi_slope) then
      guess_n_poles = np
      band_width = band_width_in*(AF_slope/Fermi_slope)
      return
    endif
    call Finalise(t_AF)
  end do

  call system_abort ("Couldn't find appropriate band_width and n_poles in guess_n_poles")

end function guess_n_poles

elemental function approx_f_Fermi(this, E)
  type(ApproxFermi), intent(in) :: this
  real(dp), intent(in) :: E
  real(dp) :: approx_f_Fermi

  integer i

  approx_f_Fermi = 0.0_dp
  do i=1, this%n_poles
    approx_f_Fermi = approx_f_Fermi + this%a(i)/(E-this%z(i))
  end do

end function approx_f_Fermi

elemental function approx_f_Fermi_deriv(this, E)
  type(ApproxFermi), intent(in) :: this
  real(dp), intent(in) :: E
  real(dp) :: approx_f_Fermi_deriv

  integer i

  approx_f_Fermi_deriv = 0.0_dp
  do i=1, this%n_poles
    approx_f_Fermi_deriv = approx_f_Fermi_deriv - this%a(i)/((E-this%z(i))**2)
  end do

end function approx_f_Fermi_deriv

subroutine calc_pole_values(this)
  type(ApproxFermi), intent(inout) :: this

  integer :: n_poles
  real(dp) :: n_poles_d

  complex(dp) :: a, b, c, aa, bb, cc, x
  integer i, j
  real(dp) i_d

  real(dp) :: gamma, x_bot, scale, shift
  ! real(dp) :: f_max

  n_poles = this%n_poles * 2
  n_poles_d = n_poles

  gamma = 3.0_dp-sqrt(8.0_dp)
  ! f_max = 1.0_dp/( ((4.0_dp*((1.0_dp-gamma)**-2) - 1.0_dp)**(n_poles_d/2)) + 1.0_dp )
  x_bot = -4.0_dp*(2.0_dp*n_poles_d-12.0_dp+6.0_dp*gamma)/((1.0_dp+gamma)**2)

  do i=1, this%n_poles/2
    i_d = i-1
    a = (1.0_dp + gamma) / (2.0_dp*n_poles_d)
    b = (1.0_dp - gamma) / (n_poles_d)
    c = cmplx( cos(PI/(n_poles_d/2.0_dp) + i_d*PI/(n_poles_d/4.0_dp)), &
	       sin(PI/(n_poles_d/2.0_dp) + i_d*PI/(n_poles_d/4.0_dp)), dp )

    aa = a*a
    bb = (2.0_dp*a + b*c)
    cc = (1.0_dp-c)

    this%z(i) = conjg ( (-bb - sqrt(bb*bb - 4.0_dp*aa*cc) ) / (2.0_dp * aa) )
    this%z(n_poles/2-i+1) = (-bb + sqrt(bb*bb - 4.0_dp*aa*cc) ) / (2.0_dp * aa)
  end do

  do i=1, this%n_poles/2
    x = this%z(i)
    this%z(i) = this%z(this%n_poles-i+1)
    this%z(this%n_poles-i+1) = x
  end do

  scale = -this%band_width/x_bot
  shift = this%Fermi_E

  do i=1, this%n_poles
    this%z(i) = this%z(i) * scale + shift
  end do

  do i=1, this%n_poles
    x = this%z(i)

    this%a(i) = (scale**(n_poles/2))*((1.0_dp-((x-shift)/scale)*(1.0_dp-gamma)/n_poles_d)**(n_poles/2))

    do j=1, this%n_poles
      if (j == i) then
	  this%a(i) = scale*(((n_poles_d*2.0_dp)/(1.0_dp+gamma))**2)*this%a(i) /  &
			      (x-conjg(this%z(j)))
      else
	  this%a(i) = scale*(((n_poles_d*2.0_dp)/(1.0_dp+gamma))**2)*this%a(i) / &
			       ( (x-this%z(j))*(x-conjg(this%z(j))) )
      endif
    end do

    this%a(i) = this%a(i) * 2.0_dp
  end do

end subroutine

end module ApproxFermi_module

max_stars_repo_name: albapa/QUIP

max_issues_repo_name: albapa/QUIP

max_forks_repo_name: albapa/QUIP