text:
module spherical_mod

!Balaji: parallel spectral model using transpose method
!initialize spherical operates on global domain
!all other operators attempt to detect whether argument is global (0:M,0:N)
!                                                       or local (ms:me,ns:ne)

use       fms_mod, only: mpp_pe, mpp_root_pe, error_mesg, FATAL, &
                         write_version_number

use  spec_mpp_mod, only: get_spec_domain

!-------------------------------------------------------------------------
!   provides operations on spectral spherical harmonics fields that do not
!      require transforms
!
!   spectral fields are complex with horizontal dimensions
!        (0:num_fourier,0:num_spherical)
!       where the zonal wavenumber of (m,n) is M = m*fourier_inc
!             the "meridional" wavenumber is n
!             the "total", 2D, spherical wavenumber L = M+n
!
!-----------------------------------------------------------------------

implicit none
private

character(len=128), parameter :: version = '$Id: spherical.f90,v 13.0 2006/03/28 21:18:30 fms Exp $'
character(len=128), parameter :: tagname ='$Name: latest $'

interface compute_lon_deriv_cos
   module procedure compute_lon_deriv_cos_2d,  &
                    compute_lon_deriv_cos_3d
end interface

interface compute_lat_deriv_cos
   module procedure compute_lat_deriv_cos_2d,  &
                    compute_lat_deriv_cos_3d
end interface

interface compute_gradient_cos
   module procedure compute_gradient_cos_2d,  &
                    compute_gradient_cos_3d
end interface

interface compute_laplacian
   module procedure compute_laplacian_2d,  &
                    compute_laplacian_3d
end interface

interface compute_ucos_vcos
   module procedure compute_ucos_vcos_2d,  &
                    compute_ucos_vcos_3d
end interface

interface compute_vor_div
   module procedure compute_vor_div_2d,  &
                    compute_vor_div_3d
end interface

interface compute_vor
   module procedure compute_vor_2d,  &
                    compute_vor_3d
end interface

interface compute_div
   module procedure compute_div_2d,  &
                    compute_div_3d
end interface

interface compute_alpha_operator
   module procedure compute_alpha_operator_2d,  &
                    compute_alpha_operator_3d
end interface

interface triangular_truncation
   module procedure triangular_truncation_2d,  &
                    triangular_truncation_3d
end interface

interface rhomboidal_truncation
   module procedure rhomboidal_truncation_2d,  &
                    rhomboidal_truncation_3d
end interface

public :: spherical_init, spherical_end
public :: compute_lon_deriv_cos, compute_lat_deriv_cos
public :: compute_laplacian, compute_vor, compute_div
public :: get_spherical_wave, get_fourier_wave
public :: get_eigen_laplacian, compute_gradient_cos
public :: compute_ucos_vcos, compute_vor_div
public :: triangular_truncation, rhomboidal_truncation

real,    allocatable, dimension(:,:) :: eigen_laplacian, epsilon, r_epsilon
integer, allocatable, dimension(:,:) :: fourier_wave, spherical_wave

integer :: fourier_max, fourier_inc, num_fourier, num_spherical
integer :: ms, me, ns, ne
logical :: module_is_initialized = .false.

real, allocatable, dimension(:,:) :: coef_uvm
real, allocatable, dimension(:,:) :: coef_uvc
real, allocatable, dimension(:,:) :: coef_uvp
real, allocatable, dimension(:,:) :: coef_alpm
real, allocatable, dimension(:,:) :: coef_alpp
real, allocatable, dimension(:,:) :: coef_dym
real, allocatable, dimension(:,:) :: coef_dx
real, allocatable, dimension(:,:) :: coef_dyp
real, allocatable, dimension(:,:) :: triangle_mask

contains

!--------------------------------------------------------------------------
subroutine spherical_init(radius, num_fourier_in, fourier_inc_in, num_spherical_in)

real,    intent(in) :: radius
integer, intent(in) :: num_fourier_in, fourier_inc_in, num_spherical_in

integer :: m,n

if(module_is_initialized) return

call write_version_number(version, tagname)

call get_spec_domain(ms, me, ns, ne)

num_fourier = num_fourier_in
fourier_inc = fourier_inc_in
fourier_max = fourier_inc*num_fourier
num_spherical = num_spherical_in

allocate(eigen_laplacian(0:num_fourier,0:num_spherical))
allocate(epsilon        (0:num_fourier,0:num_spherical))
allocate(r_epsilon      (0:num_fourier,0:num_spherical))
allocate(fourier_wave   (0:num_fourier,0:num_spherical))
allocate(spherical_wave (0:num_fourier,0:num_spherical))
allocate(coef_uvm       (0:num_fourier,0:num_spherical))
allocate(coef_uvc       (0:num_fourier,0:num_spherical))
allocate(coef_uvp       (0:num_fourier,0:num_spherical))
allocate(coef_alpm      (0:num_fourier,0:num_spherical))
allocate(coef_alpp      (0:num_fourier,0:num_spherical))
allocate(coef_dym       (0:num_fourier,0:num_spherical))
allocate(coef_dx        (0:num_fourier,0:num_spherical))
allocate(coef_dyp       (0:num_fourier,0:num_spherical))
allocate(triangle_mask  (0:num_fourier,0:num_spherical))

module_is_initialized = .true.

triangle_mask = 1.0
do n=0,num_spherical
  do m=0,num_fourier
    fourier_wave(m,n)   = m*fourier_inc
    spherical_wave(m,n) = fourier_wave(m,n) + n
    if(spherical_wave(m,n).gt.num_spherical-1) triangle_mask(m,n) = 0.0
  end do
end do

epsilon= (spherical_wave**2 - fourier_wave**2)/(4.0*spherical_wave**2 - 1.0)
epsilon = sqrt(epsilon)

eigen_laplacian = spherical_wave*(spherical_wave + 1.0)/(radius*radius)

where (spherical_wave > 0)
  coef_uvm = -radius*epsilon/spherical_wave
  coef_uvc = -radius*fourier_wave/(spherical_wave*(spherical_wave + 1.0))
else where
  coef_uvm = 0.0
  coef_uvc = 0.0
end where

coef_uvp(:,0:num_spherical-1) =    &
    -radius*epsilon(:,1:num_spherical)/   &
     (spherical_wave(:,0:num_spherical-1) +1.0)

coef_alpm= (spherical_wave + 1.0)*epsilon/radius
coef_alpp(:,0:num_spherical-1) =    &
     spherical_wave(:,0:num_spherical-1)*epsilon(:,1:num_spherical)/radius

coef_dym = (spherical_wave - 1.0)*epsilon/radius
coef_dx = float(fourier_wave)/radius
coef_dyp(:,0:num_spherical-1) =                                           &
  (spherical_wave(:,0:num_spherical-1) + 2.0)*epsilon(:,1:num_spherical)/radius

return
end subroutine spherical_init

!---------------------------------------------------------------------------
subroutine get_spherical_wave(spherical_wave_out)
!---------------------------------------------------------------------------

integer, intent(out), dimension (:,:) :: spherical_wave_out

if( size(spherical_wave_out,1).EQ.num_fourier+1 .AND. size(spherical_wave_out,2).EQ.num_spherical+1 )then
    spherical_wave_out = spherical_wave
else if( size(spherical_wave_out,1).EQ.me-ms+1 .AND. size(spherical_wave_out,2).EQ.ne-ns+1 )then
    spherical_wave_out = spherical_wave(ms:me,ns:ne)
else
    call error_mesg( 'get_spherical_wave', 'invalid argument size', FATAL )
endif

return
end subroutine get_spherical_wave

!---------------------------------------------------------------------------
subroutine get_fourier_wave(fourier_wave_out)
!---------------------------------------------------------------------------

integer, intent(out), dimension (:,:) :: fourier_wave_out

if( size(fourier_wave_out,1).EQ.num_fourier+1 .AND. size(fourier_wave_out,2).EQ.num_spherical+1 )then
    fourier_wave_out = fourier_wave
else if( size(fourier_wave_out,1).EQ.me-ms+1 .AND. size(fourier_wave_out,2).EQ.ne-ns+1 )then
    fourier_wave_out = fourier_wave(ms:me,ns:ne)
else
    call error_mesg( 'get_fourier_wave', 'invalid argument size', FATAL )
endif

return
end subroutine get_fourier_wave

!---------------------------------------------------------------------------
subroutine get_eigen_laplacian(eigen_laplacian_out)
!---------------------------------------------------------------------------

real, intent(out), dimension (:,:) :: eigen_laplacian_out

if( size(eigen_laplacian_out,1).EQ.num_fourier+1 .AND. size(eigen_laplacian_out,2).EQ.num_spherical+1 )then
    eigen_laplacian_out = eigen_laplacian
else if( size(eigen_laplacian_out,1).EQ.me-ms+1 .AND. size(eigen_laplacian_out,2).EQ.ne-ns+1 )then
    eigen_laplacian_out = eigen_laplacian(ms:me,ns:ne)
else
    call error_mesg( 'get_eigen_laplacian', 'invalid argument size', FATAL )
endif

return
end subroutine get_eigen_laplacian

!---------------------------------------------------------------------------
function compute_lon_deriv_cos_3d(spherical) result(deriv_lon)
!---------------------------------------------------------------------------

complex, intent(in), dimension (:,:,:) :: spherical
complex, dimension (size(spherical,1), size(spherical,2), size(spherical,3)) :: deriv_lon

integer :: k

if(.not. module_is_initialized ) then
  call error_mesg('compute_lon_deriv','module spherical not initialized', FATAL)
end if

if( size(spherical,1).EQ.num_fourier+1 .AND. size(spherical,2).EQ.num_spherical+1 )then
    do k=1,size(spherical,3)
       deriv_lon(:,:,k) = coef_dx(:,:)*cmplx(-aimag(spherical(:,:,k)),real(spherical(:,:,k)))
    end do
else if( size(spherical,1).EQ.me-ms+1 .AND. size(spherical,2).EQ.ne-ns+1 )then
    do k=1,size(spherical,3)
       deriv_lon(:,:,k) = coef_dx(ms:me,ns:ne)*cmplx(-aimag(spherical(:,:,k)),real(spherical(:,:,k)))
    end do
else
    call error_mesg( 'compute_lon_deriv_cos_3d', 'invalid argument size', FATAL )
endif

return
end function compute_lon_deriv_cos_3d

!--------------------------------------------------------------------------
function compute_lat_deriv_cos_3d(spherical) result(deriv_lat)
!--------------------------------------------------------------------------

complex, intent(in), dimension (:,0:,:) :: spherical
complex, dimension (size(spherical,1), 0:size(spherical,2)-1, size(spherical,3)) :: deriv_lat

integer :: k

if(.not. module_is_initialized ) then
  call error_mesg('compute_lat_deriv','module spherical not initialized', FATAL)
end if

if( size(spherical,2).EQ.num_spherical+1 )then
!could be global domain, or only global in N
    if( size(spherical,1).EQ.num_fourier+1 )then
        deriv_lat(:,0,:) = cmplx(0.,0.)
        do k=1,size(spherical,3)
           deriv_lat(:,1:num_spherical,k) = &
                - spherical(:,0:num_spherical-1,k)*coef_dym(:,1:num_spherical)
           deriv_lat(:,0:num_spherical-1,k) =  deriv_lat(:,0:num_spherical-1,k)   &
                + spherical(:,1:num_spherical,k)*coef_dyp(:,0:num_spherical-1)
        end do
    else if( size(spherical,1).EQ.me-ms+1 )then
        deriv_lat(:,0,:) = cmplx(0.,0.)
        do k=1,size(spherical,3)
           deriv_lat(:,1:num_spherical,k) = &
                - spherical(:,0:num_spherical-1,k)*coef_dym(ms:me,1:num_spherical)
           deriv_lat(:,0:num_spherical-1,k) =  deriv_lat(:,0:num_spherical-1,k)   &
                + spherical(:,1:num_spherical,k)*coef_dyp(ms:me,0:num_spherical-1)
        end do
    endif
else if( size(spherical,1).EQ.me-ms+1 .AND. size(spherical,2).EQ.ne-ns+1 )then
!need to write stuff to acquire data at ns-1,ne+1
    call abort()
else
    call error_mesg( 'compute_lat_deriv_cos_3d', 'invalid argument size', FATAL )
endif

return
end function compute_lat_deriv_cos_3d

!-------------------------------------------------------------------
subroutine compute_gradient_cos_3d(spherical, deriv_lon, deriv_lat)
!-------------------------------------------------------------------

complex, intent(in), dimension (:,:,:) :: spherical
complex, intent(out), dimension (:,:,:) :: deriv_lat
complex, intent(out), dimension (:,:,:) :: deriv_lon

deriv_lon = compute_lon_deriv_cos(spherical)
deriv_lat = compute_lat_deriv_cos(spherical)

return
end subroutine compute_gradient_cos_3d

!----------------------------------------------------------------------
function compute_laplacian_3d(spherical, power) result(laplacian)
!----------------------------------------------------------------------

  complex, intent(in), dimension (:,:,:) :: spherical
  integer, optional :: power

  complex, dimension (size(spherical,1), size(spherical,2), size(spherical,3)) :: laplacian

  integer :: k
  real, dimension(size(spherical,1), size(spherical,2)) :: factor

  if(.not. module_is_initialized ) then
      call error_mesg('compute_laplacian','module spherical not initialized', FATAL)
  end if

  if( size(spherical,1).EQ.num_fourier+1 .AND. size(spherical,2).EQ.num_spherical+1 )then
      if(present(power)) then
          if(power >= 0) then
              factor = (-eigen_laplacian)**power
          else
              where (eigen_laplacian .ne. 0.0)
                  factor = (-eigen_laplacian)**power
              else where
                  factor = 0.0
              end where
          end if
      else
          factor = -eigen_laplacian
      end if
  else if( size(spherical,1).EQ.me-ms+1 .AND. size(spherical,2).EQ.ne-ns+1 )then
      if(present(power)) then
          if(power >= 0) then
              factor = (-eigen_laplacian(ms:me,ns:ne))**power
          else
              where (eigen_laplacian(ms:me,ns:ne) .ne. 0.0)
                  factor = (-eigen_laplacian(ms:me,ns:ne))**power
              else where
                  factor = 0.0
              end where
          end if
      else
          factor = -eigen_laplacian(ms:me,ns:ne)
      end if
  else
      call error_mesg( 'compute_laplacian', 'invalid argument size', FATAL )
  endif

  do k= 1,size(spherical,3)
     laplacian(:,:,k) = spherical(:,:,k)*factor
  end do

  return
end function compute_laplacian_3d

!----------------------------------------------------------------------
subroutine compute_ucos_vcos_3d(vorticity , divergence, u_cos, v_cos)
!----------------------------------------------------------------------

complex, intent(in), dimension (:,0:,:) :: vorticity
complex, intent(in), dimension (:,0:,:) :: divergence
complex, intent(out), dimension (:,0:,:) :: u_cos
complex, intent(out), dimension (:,0:,:) :: v_cos

integer :: k

if(.not. module_is_initialized ) then
  call error_mesg('compute_ucos_vcos','module spherical not initialized', FATAL)
end if

if( size(vorticity,2).EQ.num_spherical+1 )then
!could be global domain, or only global in N
    if( size(vorticity,1).EQ.num_fourier+1 )then
        do k=1,size(vorticity,3)
           u_cos(:,:,k) = coef_uvc(:,:)*                                     &
                cmplx(-aimag(divergence(:,:,k)),real(divergence(:,:,k)))
           v_cos(:,:,k) = coef_uvc(:,:)*                                     &
                cmplx(-aimag(vorticity(:,:,k)),real(vorticity(:,:,k)))

           u_cos(:,1:num_spherical,k) = u_cos(:,1:num_spherical,k) +         &
                coef_uvm(:,1:num_spherical)*vorticity(:,0:num_spherical-1,k)
           v_cos(:,1:num_spherical,k) = v_cos(:,1:num_spherical,k) -         &
                coef_uvm(:,1:num_spherical)*divergence(:,0:num_spherical-1,k)

           u_cos(:,0:num_spherical-1,k) = u_cos(:,0:num_spherical-1,k) -     &
                coef_uvp(:,0:num_spherical-1)*vorticity(:,1:num_spherical,k)
           v_cos(:,0:num_spherical-1,k) = v_cos(:,0:num_spherical-1,k) +     &
                coef_uvp(:,0:num_spherical-1)*divergence(:,1:num_spherical,k)
        end do
    else if( size(vorticity,1).EQ.me-ms+1 )then
        do k=1,size(vorticity,3)
           u_cos(:,:,k) = coef_uvc(ms:me,:)*                                     &
                cmplx(-aimag(divergence(:,:,k)),real(divergence(:,:,k)))
           v_cos(:,:,k) = coef_uvc(ms:me,:)*                                     &
                cmplx(-aimag(vorticity(:,:,k)),real(vorticity(:,:,k)))

           u_cos(:,1:num_spherical,k) = u_cos(:,1:num_spherical,k) +         &
                coef_uvm(ms:me,1:num_spherical)*vorticity(:,0:num_spherical-1,k)
           v_cos(:,1:num_spherical,k) = v_cos(:,1:num_spherical,k) -         &
                coef_uvm(ms:me,1:num_spherical)*divergence(:,0:num_spherical-1,k)

           u_cos(:,0:num_spherical-1,k) = u_cos(:,0:num_spherical-1,k) -     &
                coef_uvp(ms:me,0:num_spherical-1)*vorticity(:,1:num_spherical,k)
           v_cos(:,0:num_spherical-1,k) = v_cos(:,0:num_spherical-1,k) +     &
                coef_uvp(ms:me,0:num_spherical-1)*divergence(:,1:num_spherical,k)
        end do
    endif
else if( size(vorticity,1).EQ.me-ms+1 .AND. size(vorticity,2).EQ.ne-ns+1 )then
!need to write stuff to acquire data at ns-1,ne+1
    call abort()
else
    call error_mesg( 'compute_ucos_vcos', 'invalid argument size', FATAL )
endif

return
end subroutine compute_ucos_vcos_3d

!-------------------------------------------------------------------------
subroutine compute_vor_div_3d(u_cos, v_cos, vorticity, divergence)
!-------------------------------------------------------------------------

complex, intent(in), dimension (:,:,:) :: u_cos
complex, intent(in), dimension (:,:,:) :: v_cos
complex, intent(out), dimension (:,:,:) :: vorticity
complex, intent(out), dimension (:,:,:) :: divergence

vorticity  = compute_alpha_operator(v_cos, u_cos, -1)
divergence = compute_alpha_operator(u_cos, v_cos, +1)

return
end subroutine compute_vor_div_3d

!-------------------------------------------------------------------------
function compute_vor_3d(u_cos, v_cos) result(vorticity)
!-------------------------------------------------------------------------

complex, intent(in), dimension (:,:,:) :: u_cos
complex, intent(in), dimension (:,:,:) :: v_cos
complex, dimension (size(u_cos,1), size(u_cos,2), size(u_cos,3)) :: vorticity

vorticity = compute_alpha_operator(v_cos, u_cos, -1)

return
end function compute_vor_3d

!-------------------------------------------------------------------------
function compute_div_3d(u_cos, v_cos) result(divergence)
!-------------------------------------------------------------------------

complex, intent(in), dimension (:,:,:) :: u_cos
complex, intent(in), dimension (:,:,:) :: v_cos
complex, dimension (size(u_cos,1), size(u_cos,2), size(u_cos,3)) :: divergence

divergence = compute_alpha_operator(u_cos, v_cos, +1)

return
end function compute_div_3d

!--------------------------------------------------------------------------------
function compute_alpha_operator_3d(spherical_a, spherical_b, isign) result(alpha)
!--------------------------------------------------------------------------------

complex, intent(in), dimension (:,0:,:) :: spherical_a
complex, intent(in), dimension (:,0:,:) :: spherical_b
integer,intent(in) :: isign

complex, dimension (size(spherical_a,1), 0:size(spherical_a,2)-1, size(spherical_a,3)) :: alpha

integer :: k

if(.not. module_is_initialized ) then
  call error_mesg('compute_vor or div','module spherical not initialized', FATAL)
end if

alpha = cmplx(0.,0.)

if( size(spherical_a,2).EQ.num_spherical+1 )then
!could be global domain, or only global in N
    if( size(spherical_a,1).EQ.num_fourier+1 )then
        do k=1,size(spherical_a,3)
           alpha(:,:,k) = coef_dx(:,:)*    &
                cmplx(-aimag(spherical_a(:,:,k)),real(spherical_a(:,:,k)))
           alpha(:,1:num_spherical,k) = alpha(:,1:num_spherical,k) -  &
                isign*coef_alpm(:,1:num_spherical)  &
                *spherical_b(:,0:num_spherical-1,k)
           alpha(:,0:num_spherical-1,k) = alpha(:,0:num_spherical-1,k) +  &
                isign*coef_alpp(:,0:num_spherical-1)*spherical_b(:,1:num_spherical,k)
        end do
    else if( size(spherical_a,1).EQ.me-ms+1 )then
        do k=1,size(spherical_a,3)
           alpha(:,:,k) = coef_dx(ms:me,:)*    &
                cmplx(-aimag(spherical_a(:,:,k)),real(spherical_a(:,:,k)))
           alpha(:,1:num_spherical,k) = alpha(:,1:num_spherical,k) -  &
                isign*coef_alpm(ms:me,1:num_spherical)  &
                *spherical_b(:,0:num_spherical-1,k)
           alpha(:,0:num_spherical-1,k) = alpha(:,0:num_spherical-1,k) +  &
                isign*coef_alpp(ms:me,0:num_spherical-1)*spherical_b(:,1:num_spherical,k)
        end do
    endif
else if( size(spherical_a,1).EQ.me-ms+1 .AND. size(spherical_a,2).EQ.ne-ns+1 )then
!need to write stuff to acquire data at ns-1,ne+1
    call abort()
else
    call error_mesg( 'compute_alpha_operator_3d', 'invalid argument size', FATAL )
endif

return
end function compute_alpha_operator_3d

!-----------------------------------------------------------------------
subroutine triangular_truncation_3d(spherical, trunc)
!-----------------------------------------------------------------------

complex, intent(inout), dimension (:,:,:) :: spherical
integer, intent(in), optional :: trunc
integer :: k

if( size(spherical,1).EQ.num_fourier+1 .AND. size(spherical,2).EQ.num_spherical+1 )then
    if(present(trunc)) then
        do k=1, size(spherical,3)
           where (spherical_wave > trunc)
               spherical(:,:,k) = cmplx(0.,0.)
           end where
        end do
    else
        do k=1, size(spherical,3)
           spherical(:,:,k) = spherical(:,:,k)*triangle_mask
        end do
    end if
else if( size(spherical,1).EQ.me-ms+1 .AND. size(spherical,2).EQ.ne-ns+1 )then
    if(present(trunc)) then
        do k=1, size(spherical,3)
           where (spherical_wave(ms:me,ns:ne) > trunc)
               spherical(:,:,k) = cmplx(0.,0.)
           end where
        end do
    else
        do k=1, size(spherical,3)
           spherical(:,:,k) = spherical(:,:,k)*triangle_mask(ms:me,ns:ne)
        end do
    end if
else
    call error_mesg( 'triang_trunc', 'invalid argument size', FATAL )
endif

return
end subroutine triangular_truncation_3d

!-----------------------------------------------------------------------------
subroutine rhomboidal_truncation_3d(spherical, trunc_fourier, trunc_spherical)
!-----------------------------------------------------------------------------

complex, intent(inout), dimension (0:,0:,:) :: spherical
integer, intent(in), optional:: trunc_fourier, trunc_spherical

integer :: k,n,m

if( size(spherical,1).EQ.num_fourier+1 .AND. size(spherical,2).EQ.num_spherical+1 )then
    if(present(trunc_fourier).and.present(trunc_spherical)) then
        do k=1,size(spherical,3)
           do n=0,size(spherical,2)-1
              do m=0,size(spherical,1)-1
                 if(fourier_wave(m,n) .GT. trunc_fourier.OR.n .GT. trunc_spherical)  &
                      spherical(m,n,k) = cmplx(0.,0.)
              end do
           end do
        end do
    else
        spherical(:,num_spherical,:) = cmplx(0.,0.)
    end if
else if( size(spherical,1).EQ.me-ms+1 .AND. size(spherical,2).EQ.ne-ns+1 )then
!wavenumbers are numbered from (0,0) which is actually (ms,ns)
    if(present(trunc_fourier).and.present(trunc_spherical)) then
        do k=1,size(spherical,3)
           do n=0,size(spherical,2)-1
              do m=0,size(spherical,1)-1
                 if(fourier_wave(m+ms,n+ns).GT.trunc_fourier .OR. n+ns.GT.trunc_spherical)  &
                      spherical(m,n,k) = cmplx(0.,0.)
              end do
           end do
        end do
    else
!truncate on n=num_spherical, only on domain containing n=num_spherical
        if( ne.EQ.num_spherical )spherical(:,ne-ns,:) = cmplx(0.,0.)
    end if
else
    call error_mesg( 'rhomboid_trunc', 'invalid argument size', FATAL )
endif

return
end subroutine rhomboidal_truncation_3d

!rest are 2d versions of 3d routines

!---------------------------------------------------------------------------
function compute_lon_deriv_cos_2d(spherical) result(deriv_lon)
!---------------------------------------------------------------------------

complex, intent(in),  dimension (:,:) :: spherical
complex, dimension (size(spherical,1), size(spherical,2)) :: deriv_lon

complex, dimension (size(spherical,1), size(spherical,2), 1) :: spherical_3d
complex, dimension (size(spherical,1), size(spherical,2), 1) :: deriv_lon_3d

spherical_3d(:,:,1) = spherical(:,:)
deriv_lon_3d = compute_lon_deriv_cos_3d(spherical_3d)
deriv_lon(:,:) = deriv_lon_3d(:,:,1)

return
end function compute_lon_deriv_cos_2d

!--------------------------------------------------------------------------
function compute_lat_deriv_cos_2d(spherical) result(deriv_lat)
!--------------------------------------------------------------------------

complex, intent(in),  dimension (:,:) :: spherical
complex, dimension (size(spherical,1), size(spherical,2)) :: deriv_lat

complex, dimension (size(spherical,1), size(spherical,2), 1) :: spherical_3d
complex, dimension (size(spherical,1), size(spherical,2), 1) :: deriv_lat_3d

spherical_3d(:,:,1) = spherical(:,:)
deriv_lat_3d = compute_lat_deriv_cos_3d(spherical_3d)
deriv_lat(:,:) = deriv_lat_3d(:,:,1)

return
end function compute_lat_deriv_cos_2d

!-------------------------------------------------------------------
subroutine compute_gradient_cos_2d(spherical, deriv_lon, deriv_lat)
!-------------------------------------------------------------------

complex, intent(in),  dimension(:,:) :: spherical
complex, intent(out), dimension(size(spherical,1), size(spherical,2)) :: deriv_lat
complex, intent(out), dimension(size(spherical,1), size(spherical,2)) :: deriv_lon

complex, dimension(size(spherical,1), size(spherical,2), 1) :: spherical_3d
complex, dimension(size(spherical,1), size(spherical,2), 1) :: deriv_lat_3d
complex, dimension(size(spherical,1), size(spherical,2), 1) :: deriv_lon_3d

spherical_3d(:,:,1) = spherical(:,:)
call compute_gradient_cos_3d(spherical_3d, deriv_lon_3d, deriv_lat_3d)
deriv_lon(:,:) = deriv_lon_3d(:,:,1)
deriv_lat(:,:) = deriv_lat_3d(:,:,1)

return
end subroutine compute_gradient_cos_2d

!----------------------------------------------------------------------
function compute_laplacian_2d(spherical, power) result(laplacian)
!----------------------------------------------------------------------

complex, intent(in), dimension (:,:) :: spherical
integer, optional :: power

complex, dimension (size(spherical,1), size(spherical,2)) :: laplacian

complex, dimension(size(spherical,1), size(spherical,2), 1) :: spherical_3d
complex, dimension(size(spherical,1), size(spherical,2), 1) :: laplacian_3d

spherical_3d(:,:,1) = spherical(:,:)
laplacian_3d = compute_laplacian_3d(spherical_3d, power)
laplacian(:,:) = laplacian_3d(:,:,1)

return
end function compute_laplacian_2d

!----------------------------------------------------------------------
subroutine compute_ucos_vcos_2d(vorticity , divergence, u_cos, v_cos)
!----------------------------------------------------------------------

complex, intent(in),  dimension (:,:) :: vorticity
complex, intent(in),  dimension (size(vorticity,1), size(vorticity,2)) :: divergence
complex, intent(out), dimension (size(vorticity,1), size(vorticity,2)) :: u_cos
complex, intent(out), dimension (size(vorticity,1), size(vorticity,2)) :: v_cos

complex, dimension (size(vorticity,1), size(vorticity,2), 1) :: vorticity_3d
complex, dimension (size(vorticity,1), size(vorticity,2), 1) :: divergence_3d
complex, dimension (size(vorticity,1), size(vorticity,2), 1) :: u_cos_3d
complex, dimension (size(vorticity,1), size(vorticity,2), 1) :: v_cos_3d

vorticity_3d(:,:,1)  = vorticity(:,:)
divergence_3d(:,:,1) = divergence(:,:)
call compute_ucos_vcos_3d(vorticity_3d, divergence_3d, u_cos_3d, v_cos_3d)
u_cos(:,:) = u_cos_3d(:,:,1)
v_cos(:,:) = v_cos_3d(:,:,1)

return
end subroutine compute_ucos_vcos_2d

!-------------------------------------------------------------------------
subroutine compute_vor_div_2d(u_div_cos, v_div_cos, vorticity, divergence)
!-------------------------------------------------------------------------

complex, intent(in),  dimension (:,:) :: u_div_cos
complex, intent(in),  dimension (size(u_div_cos,1), size(u_div_cos,2)) :: v_div_cos
complex, intent(out), dimension (size(u_div_cos,1), size(u_div_cos,2)) :: vorticity
complex, intent(out), dimension (size(u_div_cos,1), size(u_div_cos,2)) :: divergence

complex, dimension (size(u_div_cos,1), size(u_div_cos,2), 1) :: u_div_cos_3d
complex, dimension (size(u_div_cos,1), size(u_div_cos,2), 1) :: v_div_cos_3d
complex, dimension (size(u_div_cos,1), size(u_div_cos,2), 1) :: vorticity_3d
complex, dimension (size(u_div_cos,1), size(u_div_cos,2), 1) :: divergence_3d

u_div_cos_3d(:,:,1) = u_div_cos(:,:)
v_div_cos_3d(:,:,1) = v_div_cos(:,:)
call compute_vor_div_3d(u_div_cos_3d, v_div_cos_3d, vorticity_3d, divergence_3d)
vorticity(:,:)  = vorticity_3d(:,:,1)
divergence(:,:) = divergence_3d(:,:,1)

return
end subroutine compute_vor_div_2d

!-------------------------------------------------------------------------
function compute_vor_2d(u_div_cos, v_div_cos) result(vorticity)
!-------------------------------------------------------------------------

complex, intent(in),  dimension (:,:) :: u_div_cos
complex, intent(in),  dimension (:,:) :: v_div_cos
complex, dimension (size(u_div_cos,1), size(u_div_cos,2)) :: vorticity

complex, dimension (size(u_div_cos,1), size(u_div_cos,2), 1) :: u_div_cos_3d
complex, dimension (size(u_div_cos,1), size(u_div_cos,2), 1) :: v_div_cos_3d
complex, dimension (size(u_div_cos,1), size(u_div_cos,2), 1) :: vorticity_3d

u_div_cos_3d(:,:,1) = u_div_cos(:,:)
v_div_cos_3d(:,:,1) = v_div_cos(:,:)
vorticity_3d = compute_vor_3d(u_div_cos_3d, v_div_cos_3d)
vorticity(:,:) = vorticity_3d(:,:,1)

return
end function compute_vor_2d

!-------------------------------------------------------------------------
function compute_div_2d(u_div_cos, v_div_cos) result(divergence)
!-------------------------------------------------------------------------

complex, intent(in),  dimension (:,:) :: u_div_cos
complex, intent(in),  dimension (:,:) :: v_div_cos
complex, dimension (size(u_div_cos,1), size(u_div_cos,2)) :: divergence

complex, dimension (size(u_div_cos,1), size(u_div_cos,2), 1) :: u_div_cos_3d
complex, dimension (size(u_div_cos,1), size(u_div_cos,2), 1) :: v_div_cos_3d
complex, dimension (size(u_div_cos,1), size(u_div_cos,2), 1) :: divergence_3d

u_div_cos_3d(:,:,1) = u_div_cos(:,:)
v_div_cos_3d(:,:,1) = v_div_cos(:,:)
divergence_3d = compute_div_3d(u_div_cos_3d, v_div_cos_3d)
divergence(:,:) = divergence_3d(:,:,1)

return
end function compute_div_2d

!--------------------------------------------------------------------------------
function compute_alpha_operator_2d(spherical_a, spherical_b, isign) result(alpha)
!--------------------------------------------------------------------------------

complex, intent(in),  dimension (:,:) :: spherical_a
complex, intent(in),  dimension (:,:) :: spherical_b
integer, intent(in) :: isign

complex, dimension (size(spherical_a,1), size(spherical_a,2)) :: alpha

complex, dimension (size(spherical_a,1), size(spherical_a,2), 1) :: spherical_a_3d
complex, dimension (size(spherical_a,1), size(spherical_a,2), 1) :: spherical_b_3d
complex, dimension (size(spherical_a,1), size(spherical_a,2), 1) :: alpha_3d

spherical_a_3d(:,:,1) = spherical_a(:,:)
spherical_b_3d(:,:,1) = spherical_b(:,:)
alpha_3d = compute_alpha_operator_3d(spherical_a_3d, spherical_b_3d, isign)
alpha(:,:) = alpha_3d(:,:,1)

return
end function compute_alpha_operator_2d

!-----------------------------------------------------------------------
subroutine triangular_truncation_2d(spherical, trunc)
!-----------------------------------------------------------------------

complex, intent(inout), dimension (:,:) :: spherical
integer, intent(in), optional :: trunc
complex, dimension (size(spherical,1), size(spherical,2), 1) :: spherical_3d

spherical_3d(:,:,1) = spherical(:,:)
call triangular_truncation_3d(spherical_3d, trunc)
spherical(:,:) = spherical_3d(:,:,1)

return
end subroutine triangular_truncation_2d

!--------------------------------------------------------------------------------
subroutine rhomboidal_truncation_2d(spherical, trunc_fourier, trunc_spherical)
!--------------------------------------------------------------------------------

complex, intent(inout), dimension (:,:) :: spherical
integer, intent(in), optional :: trunc_fourier, trunc_spherical

complex, dimension (size(spherical,1), size(spherical,2), 1) :: spherical_3d

spherical_3d(:,:,1) = spherical(:,:)
call rhomboidal_truncation_3d(spherical_3d, trunc_fourier, trunc_spherical)
spherical(:,:) = spherical_3d(:,:,1)

return
end subroutine rhomboidal_truncation_2d

!-----------------------------------------------------------------------
subroutine spherical_end

if(.not.module_is_initialized) return

deallocate(eigen_laplacian, epsilon, r_epsilon, fourier_wave, spherical_wave)
deallocate(coef_uvm, coef_uvc, coef_uvp, coef_alpm, coef_alpp, coef_dym, coef_dx, coef_dyp, triangle_mask)
module_is_initialized = .false.

return
end subroutine spherical_end
!-----------------------------------------------------------------------

end module spherical_mod

max_stars_repo_name: JoyMonteiro/FMSDynamicalCore


text:
      subroutine mboost(m,vec,beta,vin,vout)
c     boosts the m vectors vin(0:3,m) into the vectors vout(0:3,m) (that can
c     be the same) in the direction of vec(3) (|vec|=1) with velocity
c     beta.  Lorents convention: (t,x,y,z).
      implicit none
      integer m
      real * 8 vec(3),beta,vin(0:3,m),vout(0:3,m)
      real * 8 betav,gamma
      real * 8 vdotb
      integer ipart,idim
      gamma=1/sqrt(1-beta**2)
      do ipart=1,m
         vdotb=vin(1,ipart)*vec(1)
     #         +vin(2,ipart)*vec(2)+vin(3,ipart)*vec(3)
         do idim=1,3
            vout(idim,ipart)=vin(idim,ipart)
     #           +vec(idim)*((gamma-1)*vdotb
     #           +gamma*beta*vin(0,ipart))
         enddo
         vout(0,ipart)=gamma*(vin(0,ipart)+vdotb*beta)
      enddo
      end

      subroutine mrotate(dir,sinphi,cosphi,vec)
c Rotates vector vec counterclockwise around the direction
c dir (|dir|=1) with angle phi, given sin phi and cos phi.
      implicit none
      real * 8 sinphi,cosphi,dir(3),vec(3)
      real * 8 dircrossvec(3),dirdotvec
      integer i
      dircrossvec(1)=dir(2)*vec(3)-dir(3)*vec(2)
      dircrossvec(2)=dir(3)*vec(1)-dir(1)*vec(3)
      dircrossvec(3)=dir(1)*vec(2)-dir(2)*vec(1)
      dirdotvec=dir(1)*vec(1)+dir(2)*vec(2)+dir(3)*vec(3)
      do i=1,3
         vec(i)=vec(i)+sinphi*dircrossvec(i)
     #        -(1-cosphi)*(vec(i)-dir(i)*dirdotvec)
      enddo
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text: #redirect wiki:woodland:Dragon Fire Toy Company

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
The Davis Blue Devils Blue Devil Football Backers is a NonProfit Organizations notforprofit Parent/Guardian run organization dedicated to making each Davis High School DHS Football football season successful, instructive, and most importantly, fun for all of the players, spirit leaders, parents, guardians and coaches.

20070223 12:09:23 nbsp If your a backer you get a spot on the 50 yard line along with all the parents of players! wooo prime seats! Users/StevenDaubert


text:
c>>   1996-06-24  DREXSORT  Krogh   Added code for conversion to C.
c>>   1995-05-28  DREXSORT  Krogh   Converted SFTRAN to Fortran
c>>   1990-02-09  DREXSORT  Snyder  Initial code.
c
c     Test driver for EXSORT.
c
c     Sort LDATA=10000 random numbers using EXSORT.
c     Check whether they are in order.
c
      integer LENBUF
      parameter (LENBUF=1000)
      integer L(1:LENBUF),OPTION,OUTFIL
      parameter (OPTION=1)
      external DATAOP
      real R(1:LENBUF)
      common /RCOM/ R
c
      call exsort (dataop,LENBUF,l,option,outfil)
      stop
      end

      subroutine DATAOP (IOP, I, J, IFLAG)
      integer LDATA, LENBUF
      parameter (LDATA=10000, LENBUF=1000)
      integer IOP, I, J, IFLAG
      logical ISopen(4)
      save ISOPEN
      integer NCOMP, NDATA
      logical OK
      save NCOMP, NDATA, OK
      real PREV,R(1:LENBUF),SRANU
      save PREV
      external SRANU
      common /RCOM/ R
      data ISOPEN, NCOMP, NDATA, OK, PREV /4*.FALSE., 0, 0, .TRUE., -1./
c++ CODE for .C. is inactive
c%% static float end_of_seq[1] = {-1.0e0};
c%% static char *fname[4]={"scratch1","scratch2","scratch3","scratch4"};
c%% static FILE *fp[4];
c++ END
c
      go to (10, 20, 30, 40, 50, 60, 70, 80), IOP
      return
c      case 1 @ initial input into record J
   10  continue
         ndata = ndata + 1
         if (ndata.le.LDATA) then
            r(j)=sranu()
            iflag=0
         else
            iflag=1
         end if
         return
c      case 2 @ write scratch from J onto file I
   20  continue
c%%      fwrite( &rcom.r[j-1], sizeof(rcom.r[0]), 1L, fp[i-1]);
         write (i+10) r(j)
         return
c      case 3 @ write end-of-sequence onto file I
   30  continue
c%%      fwrite( end_of_seq, sizeof(rcom.r[0]), 1L, fp[i-1]);
         write (i+10) -1.0E0
         return
c      case 4 @ read scratch into J from file I
   40  continue
c%%      fread( &rcom.r[j-1], sizeof(rcom.r[0]), 1L, fp[i-1]);
         read (i+10) r(j)
         iflag=0
         if (r(j).lt.0.0) iflag=1
         return
c      case 5 @ rewind file I
   50  continue
         if (.not. isopen(i)) then
c%%         fp[i-1] = fopen(fname[i-1], "wb+");
            open (i+10, status='scratch',form='unformatted')
            isopen(i)=.true.
         end if
c%%      rewind( fp[i-1] );
         rewind (i+10)
         return
c      case 6 @ output from record J
   60  continue
         if (i.ne.0) then
            if (ok) then
               print '('' EXSORT succeeded using'',i7,'' compares'')',
     1            ncomp
            else
               print *,'EXSORT failed'
            end if
         else
            if (r(j).lt.prev) ok=.false.
            prev=r(j)
         end if
         return
c      case 7 @ move record I to record J
   70  continue
         r(j)=r(i)
         return
c      case 8 @ iflage I and J
   80  continue
         ncomp=ncomp+1
         if (r(i)-r(j)) 110,120,130
  110    iflag=-1
         return
  120    iflag=0
         return
  130    iflag=+1
         return
c
      end


text:
C
C	$Id: gflcnv.f,v 1.6 2008-07-27 00:21:06 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
      SUBROUTINE GFLCNV (FNUM,INUM)
C
C     This subroutine converts the floating-point number in
C     FNUM into two 16-bit integers in the low order bits of
C     INUM.  These two 16-bit quantities are suitable for
C     insertion into the CGM.  The first 16-bit integer is
C     given by the low-order 15 bits resulting from INUM(1)=FUNM.
C     This number is then expressed in 2's complement format.  The
C     second 16-bit integer is the decimal part of FNUM
C     multiplied by 2**16 (65536.)
C
      REAL     FNUM
      INTEGER  INUM(2)
C
      INTEGER  JNUM,JABS
C
C  Find the smallest integer less than or equal to FNUM and store
C  it in JNUM.
C
      JNUM = INT(FNUM)
      IF (FNUM.LT.0. .AND. REAL(JNUM).NE.FNUM) JNUM = JNUM-1
C
C  Find the fractional part.
C
      INUM(2) = 65536.*(FNUM-REAL(JNUM))
C
C  Define 2's complement of whole part.
C
      JABS = ABS(JNUM)
      INUM(1) = IAND(32767,JABS)
      IF (JNUM .LT. 0) THEN
        INUM(1) = (32767-INUM(1))+1
        INUM(1) = IOR(INUM(1),ISHIFT(1,15))
      ENDIF
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      program t
      implicit none
! inquire stmt all 'access' and 'form' specifiers (sequential formatted)
! and character-variables in open specifiers
      character*11::a,b,c,x,y,z
      character*10::access='sequential'
      open (56, status='new', file='tmpfile', access=access, for
     "m='formatted')
      inquire (56, access=a, sequential=b, direct=c, form=x, form
     *atted=y, unformatted=z)

      print *,a
      print *,b
      print *,c
      print *,x
      print *,y
      print *,z

      close (56,status='delete')

      endprogram t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
C
C      **********************************************************************
C      ********************** SUBROUTINE CENS  ******************************
C      **********************************************************************
C
        SUBROUTINE CENS(X,Y,IND,NTOT)
C
C      *         THIS SUBROUTINE ADDS OR SUBTRACTS 0.00001 TIMES THE        *
C      *         VALUE FROM CENSORED POINTS SO THAT IF THERE ARE TIES WITH  *
C      *         DETECTED POINTS, THE CENSORED POINTS CAN BE                *
C      *         DISTINGUISHED FROM THE DETECTIONS.                         *
C      *         IF THE SMALLEST DIGIT IS LESS THAN OR EQUAL TO             *
C      *         0.0001, THEN 'CONST' SHOULD BE CHANGED TO A                *
C      *         SMALLER VALUE.                                             *
C      *                                                                    *
C      *         INPUT AND OUTPUT:                                          *
C      *                   X(I)    : INDEPENDENT VARIABLE                   *
C      *                   Y(I)    : DEPENDENT VARIABLE                     *
C      *                  IND(I)   : CENSORED STATUS INDICATOR              *
C
        IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N)
        DIMENSION X(NTOT),Y(NTOT),IND(NTOT)
C
C
        CONST=1.000001
        DO 500 I=1,NTOT
C
C      *                 UPPER LIMITS CASES                                 *
C
        IF(IND(I).EQ.-1) THEN
              Y(I)=Y(I)/CONST

          ELSEIF(IND(I).EQ.-2) THEN
              X(I)=X(I)/CONST

          ELSEIF(IND(I).EQ.-3) THEN
              X(I)=X(I)/CONST
              Y(I)=Y(I)/CONST

          ELSEIF(IND(I).EQ.-4) THEN
              X(I)=X(I)/CONST
              Y(I)=Y(I)*CONST

C
C      *                  LOWER LIMIT CASES                                 *
C
          ELSEIF(IND(I).EQ.1)  THEN
              Y(I)=Y(I)*CONST

          ELSEIF(IND(I).EQ.2)  THEN
              X(I)=X(I)*CONST

          ELSEIF(IND(I).EQ.3)  THEN
              X(I)=X(I)*CONST
              Y(I)=Y(I)*CONST

          ELSEIF(IND(I).EQ.4)  THEN
              X(I)=X(I)*CONST
              Y(I)=Y(I)/CONST

        ENDIF
  500   CONTINUE
        RETURN
        END

max_stars_repo_name: iraf-community/stsdas


text:
C
C
C  Copyright (C) 2000, 2001 Silicon Graphics, Inc.  All Rights Reserved.
C
C  This program is free software; you can redistribute it and/or modify it
C  under the terms of version 2.1 of the GNU Lesser General Public License
C  as published by the Free Software Foundation.
C
C  This program is distributed in the hope that it would be useful, but
C  WITHOUT ANY WARRANTY; without even the implied warranty of
C  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
C
C  Further, this software is distributed without any warranty that it is
C  free of the rightful claim of any third person regarding infringement
C  or the like.  Any license provided herein, whether implied or
C  otherwise, applies only to this software file.  Patent licenses, if
C  any, provided herein do not apply to combinations of this program with
C  other software, or any other product whatsoever.
C
C  You should have received a copy of the GNU Lesser General Public
C  License along with this program; if not, write the Free Software
C  Foundation, Inc., 59 Temple Place - Suite 330, Boston MA 02111-1307,
C  USA.
C
C  Contact information:  Silicon Graphics, Inc., 1600 Amphitheatre Pky,
C  Mountain View, CA 94043, or:
C
C  http://www.sgi.com
C
C  For further information regarding this notice, see:
C
C  http://oss.sgi.com/projects/GenInfo/NoticeExplan
C
C

	subroutine eoshift_spec_wd1_p1@(result, source, dim, sext, shftval,
     &		bndval)
*
CDIR$ ID "@(#) libfi/mpp_shift/eoshift_p_spec_1@.f	92.0	"
CDIR$ ID "10/08/98 14:37:14"
*
**************************************************************************
*
* Purpose:	Perform a end-off shift of an array in argument SOURCE
*		by SHIFT elements.  This routine handles 1-word source
*		source arrays of rank 1.  This routine is called from the
*		C routine _EOSHIFT_[JLS]P1.
*
* Arguments:	result	- integer array of rank 1.
*		source	- integer array of rank 1.
*		dim	- integer scalar: dimension to shift over
*		sext	- integer array of rank 1.
*		shftval	- integer scalar: scalar shift value
*		bndval	- integer scalar: boundary value
*
* Variables:	mype	- integer scalar: logical PE number.
*		blksize	- integer scalar: block size on PE
*		pe_low	- integer scalar: PE receiving data on +1 shift
*		pe_hi	- integer scalar: PE receiving data on -1 shift
*		low	- integer scalar: low index on PE
*		hi	- integer scalar: high index on PE
*
* Intrinsic functions called:
*		my_pe()
*		blkct()
*		lowidx()
*		hiidx()
*		home()
*
* Documentation:
*		T3D Fortran 90 Array Intrinsics Interface Description
*		T3D Fortran 90 Array Intrinsics Design Description
*************************************************************************
	implicit	none
	intrinsic	my_pe, blkct, lowidx, hiidx, home
	integer		sext(1)
	integer		result(sext(1))
	integer		source(sext(1))
	integer		dim
	integer		shftval
	integer		bndval
	integer		blksize, mype
	integer		pe_low, pe_hi
	integer		low, hi
*
cdir$	unknown_shared	source, result
*
	mype = my_pe()
	if (blkct (source, 1, mype) .gt. 0) then
	    low = lowidx(source,1,mype,1)
	    hi = hiidx (source,1,mype,1)
	    blksize = hi - low + 1
	    if (low .gt. 1) then
		pe_low = home (source(low-1))
	    else
		pe_low = -1
	    endif
	    if (hi .lt. sext(1)) then
		pe_hi = home (source(hi+1))
	    else
		pe_hi = -1
	    endif
*
	    call eoshift_spec_i_wd1_p1@ (result, source, sext, shftval, bndval,
     &		blksize, pe_low, pe_hi)
	endif
*
cdir$	barrier
	call shmem_udcflush()
*
	return
	end

	subroutine eoshift_spec_i_wd1_p1@ (result, source, sext, shftval,
     &		bndval,	blksize, pe_low, pe_hi)
*
**************************************************************************
*
* Purpose:	Coerce a shared array to private, and call a work routine
*		which will perform a eoshift on the private arrays.
*
* Arguments:	result	- integer array of rank 1.
*		source	- integer array of rank 1.
*		dim	- integer scalar: dimension to shift over
*		sext	- integer scalar: extent
*		shftval	- integer scalar: scalar shift value
*		bndval	- integer scalar: boundary value
*		blksize	- integer scalar: block sizes of array
*		pe_low	- integer scalar: PE receiving data on +1 shift
*		pe_hi	- integer scalar: PE receiving data on -1 shift
*
* Variables:	None
*
* Intrinsic functions called:
*		my_pe()
*
* Documentation:
*		T3D Fortran 90 Array Intrinsics Interface Description
*		T3D Fortran 90 Array Intrinsics Design Description
*************************************************************************
*
	implicit	none
	integer		sext(1)
	integer		result(sext(1))
	integer		source(sext(1))
	integer		shftval
	integer		bndval
	integer		blksize
	integer		pe_low, pe_hi
*
cdir$	unknown_shared	source, result
cdir$	pe_resident	source, result
*
	call eoshift_spec_w_wd1_p1@ (result(1), source(1), blksize, shftval,
     &		bndval, pe_low, pe_hi)
*
	return
	end

	subroutine eoshift_spec_w_wd1_p1@ (result, source, sext, shftval,
     &		bndval, pe_low, pe_hi)
*
**************************************************************************
*
* Purpose:	Perform a eoshift on a coerced shared to private section
*		of a rank 1 array.  This routine will only be called for
*		scalar shift counts of -1, 0, or 1.  The PE number of the
*		PE which will receive the data passed out of this block
*		is found in pe_low for shifts of 1, and pe_hi for shifts
*		of -1.  Data will be transferred to the other PE's via
*		shmem_iputs.
*
*		This routine is divided into three main sections: shift
*		counts of -1, 1, and 0.  Within those sections are two
*		major divisions: Entire extent of dimension being eoshifted
*		over on 1 PE, and extent being spread over 2 or more PEs.
*		Lastly, these divisions have two parts: dim 1 and dim 2.
*
* Arguments:	result	- integer array of rank 1.
*		source	- integer array of rank 1.
*		dim	- integer scalar: dimension to shift over
*		sext	- integer scalar: extent
*		shftval	- integer scalar: scalar shift value
*		bndval	- integer scalar: boundary value
*		blksize	- integer scalar: block size of array
*		pe_low	- integer scalar: PE receiving data on +1 shift
*		pe_hi	- integer scalar: PE receiving data on -1 shift
*
* Variables:	mype	- integer scalar: PE on which we are working
*		str1	- integer scalar: number of words in 1 element
*
* Intrinsic functions called:
*		my_pe()
*
* Documentation:
*		T3D Fortran 90 Array Intrinsics Interface Description
*		T3D Fortran 90 Array Intrinsics Design Description
*************************************************************************
*
	implicit	none
	intrinsic	my_pe
	integer		sext
	integer		result(sext)
	integer		source(sext)
	integer		shftval
	integer		bndval
	integer		pe_low, pe_hi
	integer		i, j
	integer		str1
	integer		mype
*
cdir$	pe_private	source, result
*
	mype = my_pe()
*
	str1 = 1
	include "eoshift_spec_1@.fh"
*
	return
	end

	subroutine eoshift_spec_wd2_p1@(result, source, dim, sext, shftval,
     &		bndval)
*
**************************************************************************
*
* Purpose:	Perform a end-off shift of an array in argument SOURCE
*		by SHIFT elements.  This routine handles
*		complex source arrays of rank 1.  This routine is called
*		from the C routine _EOSHIFT_JP1.
*
* Arguments:	result	- complex array of rank 1.
*		source	- complex array of rank 1.
*		dim	- integer scalar: dimension to shift over
*		sext	- integer array of rank 1.
*		shftval	- integer scalar: scalar shift value
*		bndval	- complex scalar: boundary value
*
* Variables:	mype	- integer scalar: logical PE number.
*		blksize	- integer scalar: block size on PE
*		pe_low	- integer scalar: PE receiving data on +1 shift
*		pe_hi	- integer scalar: PE receiving data on -1 shift
*		low	- integer scalar: low index on PE
*		hi	- integer scalar: high index on PE
*
* Intrinsic functions called:
*		my_pe()
*		blkct()
*		lowidx()
*		hiidx()
*		home()
*
* Documentation:
*		T3D Fortran 90 Array Intrinsics Interface Description
*		T3D Fortran 90 Array Intrinsics Design Description
*************************************************************************
	implicit	none
	intrinsic	my_pe, blkct, lowidx, hiidx, home
	integer		sext(1)
	complex		result(sext(1))
	complex		source(sext(1))
	integer		dim
	integer		shftval
	complex		bndval
	integer		blksize, mype
	integer		pe_low, pe_hi
	integer		low, hi
*
cdir$	unknown_shared	source, result
*
	mype = my_pe()
	if (blkct (source, 1, mype) .gt. 0) then
	    low = lowidx(source,1,mype,1)
	    hi = hiidx (source,1,mype,1)
	    blksize = hi - low + 1
	    if (low .gt. 1) then
		pe_low = home (source(low-1))
	    else
		pe_low = -1
	    endif
	    if (hi .lt. sext(1)) then
		pe_hi = home (source(hi+1))
	    else
		pe_hi = -1
	    endif
*
	    call eoshift_spec_i_wd2_p1@ (result, source, sext, shftval, bndval,
     &		blksize, pe_low, pe_hi)
	endif
*
cdir$	barrier
	call shmem_udcflush()
*
	return
	end

	subroutine eoshift_spec_i_wd2_p1@ (result, source, sext, shftval,
     &		bndval,	blksize, pe_low, pe_hi)
*
**************************************************************************
*
* Purpose:	Coerce a shared array to private, and call a work routine
*		which will perform a eoshift on the private arrays.
*
* Arguments:	result	- complex array of rank 1.
*		source	- complex array of rank 1.
*		dim	- integer scalar: dimension to shift over
*		sext	- integer scalar: extent
*		shftval	- integer scalar: scalar shift value
*		bndval	- complex scalar: boundary value
*		blksize	- integer scalar: block sizes of array
*		pe_low	- integer scalar: PE receiving data on +1 shift
*		pe_hi	- integer scalar: PE receiving data on -1 shift
*
* Variables:	None
*
* Intrinsic functions called:
*		my_pe()
*
* Documentation:
*		T3D Fortran 90 Array Intrinsics Interface Description
*		T3D Fortran 90 Array Intrinsics Design Description
*************************************************************************
*
	implicit	none
	integer		sext(1)
	complex		result(sext(1))
	complex		source(sext(1))
	integer		shftval
	complex		bndval
	integer		blksize
	integer		pe_low, pe_hi
*
cdir$	unknown_shared	source, result
cdir$	pe_resident	source, result
*
	call eoshift_spec_w_wd2_p1@ (result(1), source(1), blksize, shftval,
     &		bndval, pe_low, pe_hi)
*
	return
	end

	subroutine eoshift_spec_w_wd2_p1@ (result, source, sext, shftval,
     &		bndval, pe_low, pe_hi)
*
**************************************************************************
*
* Purpose:	Perform a eoshift on a coerced shared to private section
*		of a rank 1 array.  This routine will only be called for
*		scalar shift counts of -1, 0, or 1.  The PE number of the
*		PE which will receive the data passed out of this block
*		is found in pe_low for shifts of 1, and pe_hi for shifts
*		of -1.  Data will be transferred to the other PE's via
*		shmem_iputs.
*
*		This routine is divided into three main sections: shift
*		counts of -1, 1, and 0.  Within those sections are two
*		major divisions: Entire extent of dimension being eoshifted
*		over on 1 PE, and extent being spread over 2 or more PEs.
*		Lastly, these divisions have two parts: dim 1 and dim 2.
*
* Arguments:	result	- complex array of rank 1.
*		source	- complex array of rank 1.
*		dim	- integer scalar: dimension to shift over
*		sext	- integer scalar: extent
*		shftval	- integer scalar: scalar shift value
*		bndval	- complex scalar: boundary value
*		blksize	- integer scalar: block size of array
*		pe_low	- integer scalar: PE receiving data on +1 shift
*		pe_hi	- integer scalar: PE receiving data on -1 shift
*
* Variables:	mype	- integer scalar: PE on which we are working
*		str1	- integer scalar: number of words in 1 element
*
* Intrinsic functions called:
*		my_pe()
*
* Documentation:
*		T3D Fortran 90 Array Intrinsics Interface Description
*		T3D Fortran 90 Array Intrinsics Design Description
*************************************************************************
*
	implicit	none
	intrinsic	my_pe
	integer		sext
	complex		result(sext)
	complex		source(sext)
	integer		shftval
	integer		bndval
	integer		pe_low, pe_hi
	integer		i, j
	integer		str1
	integer		mype
*
cdir$	pe_private	source, result
*
	mype = my_pe()
*
	str1 = 2
	include "eoshift_spec_1@.fh"
*
	return
	end

max_stars_repo_name: sharugupta/OpenUH

max_issues_repo_name: sharugupta/OpenUH

max_forks_repo_name: sharugupta/OpenUH


text:
Davis Commons is the small Shopping Centers shopping center on the corner of 1st St. and Richards Blvd. at E St.

Davis Commons was completed in 1997, amid protest that the arrival of Davis first corporate chain bookstore, Borders, would drive out smaller independent booksellers. (Several did actually Departed Businesses close shortly thereafter, including The Next Chapter, Coffee & Classics, and Gayles Books). And then in 2011 Borders (the whole chain) itself closed, with Whole Foods Market Whole Foods as the new tenant, amid worry that its arrival would threaten the Davis Food Coop, Nugget, and other preexisting grocery stores, especially while the town was still absorbing the recent arrival of Trader Joes (often see as a competitor to Whole Foods).

Davis Commons was featured on the ABC Nightly News as an example of appropriate development for the 21st Century. (In all fairness, though, Davis Commons is preferable to the standard shoeboxstyle strip mall, like Fifth Street Plaza.)

The tract of land that is now Davis Commons was once a Native American village, which was another reason for controversy around the building of the center. Interestingly, the burial ground for the village was also built over; its called the Mondavi Center.

Mark Friedman, president of Fulcrum Properties, built and owns Davis Commons; it was developed by two local firms, Pyramid Construction and Cunningham Engineering Corporation.

It consists of the following stores and restaurants:

  Bath and Body Works
  Habit Burger
  The Gap
  The Melt
  Mikuni, a sushi/Japanese restaurant chain
  Pinkberry
  Plutos
  Verizon Wireless EVO2 Store
  Whole Foods Market

Former Stores:
  AT&T Wireless  Became Cingular and then became AT&T again.
  Cingular Wireless  Formerly AT&T Wireless.
  Ben & Jerrys
  Borders
  DeVons Jewelers
  Fuzio
  Jamba Juice
  Papyrus
  Spirit Halloween Occupied the space left by Borders after Borders closed until just after Halloween.

The Arboretum Terrace on the east side of the complex has some great outside seating. At the back of the parking lot you can access the Arboretum arboretum Bike Paths bikeways.

Does anyone have a picture of the sculpture of a man in a box that used to be in at the end of the pictured walkway?

  I remember it being a woman?  Of course, that was a long time ago.
  Im pretty sure it was a nude man  Users/StevenDaubert


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 BLDSRF(XA,YA,ZA,XYZSRF)

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

C BLDSRF CREATES ARRAYS XYZSRF

C Called by MAPVAR

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

      include 'amesh.blk'

      DIMENSION XA(*),YA(*),ZA(*)
      DIMENSION XYZSRF(NODESA,3)

      DO 10 I = 1, NODESA
        XYZSRF(I,1) = XA(I)
        XYZSRF(I,2) = YA(I)
        IF (NDIMA .EQ. 2)THEN
          XYZSRF(I,3) = 0.
        ELSE
          XYZSRF(I,3) = ZA(I)
        END IF
   10 CONTINUE

      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
      subroutine asgrfm(da,fit)
c
c     assign next family member for random i/o
c
      implicit integer(a-z)
      dimension fit(8)
c
      common/frfcm1/mxfrf,ifrf,buflen,fcsize,dskloc ,curlen,kop,ier
      character*8 frfn,frn,kfn
      common/frfcm2/frfn(2,16),frn,kfn
      common/double/iprec,ncpw,unit
c
      logical fxist
      character nfn*8,msg*48
c.... parameter giving number of record units per integer word
c     for most systems a single character is used as a record unit
      data msg/'read requested from nonexistent family member - '/
c

c.... compute family member index & bias disk address for correct access
      i=da/(fcsize)
      da=da-i*fcsize

c.... get the name of the requested family member
      call nrfnam(frn,i,nfn)
c.... return if current family member is the desired one
      if(kfn.eq.nfn) return
c.... flush the buffer if data is present which is not on disk
      if (fit(3).lt.0) then
      fit(3)=-fit(3)
      call wdiska (fit(1),fit(8),buflen,dskloc )
      fit(7)=max(fit(7),dskloc +buflen)
      endif
c.... determine if requested family member exists
      if (kop.ne.0) go to 20
      inquire (file=nfn,exist=fxist)
      if (fxist) go to 20
      if (ier.ne.0) return
c      call abort (msg)
c.... close the current family member
   20 continue
      if (kfn.ne.'        ') then
      close (fit(1),status='keep')
      endif
c.... open/create the requested family member
c     lrecl=buflen
      lrecl=ncpw*buflen
      open (fit(1),file=nfn,access='direct',form='unformatted'
     1,recl=lrecl,status='unknown')
      kfn=nfn
      frfn(2,ifrf)=kfn
      fit(5)=fit(4)
      dskloc =fit(5)
      fit(6)=0
      curlen=fit(6)
      fit(7)=0
      if (fxist) fit(7)=fit(4)
      ier=0
      return
      end

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode


text:
      SUBROUTINE STEP2LON (XSTA,T,XCORSTA)
*+
*  - - - - - - - - - - -
*   S T E P 2 L O N
*  - - - - - - - - - - -
*
*  This routine is part of the International Earth Rotation and
*  Reference Systems Service (IERS) Conventions software collection.
*
*  This subroutine gives the in-phase and out-of-phase corrections
*  induced by mantle anelasticity in the long period band.
*
*  In general, Class 1, 2, and 3 models represent physical effects that
*  act on geodetic parameters while canonical models provide lower-level
*  representations or basic computations that are used by Class 1, 2, or
*  3 models.
*
*  Status: Class 1
*
*     Class 1 models are those recommended to be used a priori in the
*     reduction of raw space geodetic data in order to determine
*     geodetic parameter estimates.
*     Class 2 models are those that eliminate an observational
*     singularity and are purely conventional in nature.
*     Class 3 models are those that are not required as either Class
*     1 or 2.
*     Canonical models are accepted as is and cannot be classified as a
*     Class 1, 2, or 3 model.
*
*  Given:
*     XSTA          d(3)   Geocentric position of the IGS station (Note 1)
*     T             d      Centuries since J2000
*
*  Returned:
*     XCORSTA       d(3)   In phase and out of phase station corrections
*                          for diurnal band (Note 2)
*
*  Notes:
*
*  1) The IGS station is in ITRF co-rotating frame.  All coordinates are
*     expressed in meters.
*
*  2) All coordinates are expressed in meters.
*
*  Test case:
*     given input: XSTA(1) = 4075578.385D0 meters
*                  XSTA(2) =  931852.890D0 meters
*                  XSTA(3) = 4801570.154D0 meters
*                  T       = 0.1059411362080767D0 Julian centuries
*
*     expected output:  XCORSTA(1) = -0.9780962849562107762D-04 meters
*                       XCORSTA(2) = -0.2236349699932734273D-04 meters
*                       XCORSTA(3) =  0.3561945821351565926D-03 meters
*
*  References:
*
*     Mathews, P. M., Dehant, V., and Gipson, J. M., 1997, ''Tidal station
*     displacements," J. Geophys. Res., 102(B9), pp. 20,469-20,477
*
*     Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
*     IERS Technical Note No. 36, BKG (2010)
*
*  Revisions:
*  1996 March    23 V. Dehant      Original code
*  2009 August   07 B.E. Stetzler  Initial standardization of code
*                                  and found unnecessary variables tau
*                                  and fhr
*  2009 August   07 B.E. Stetzler  Provided a test case
*  2009 August   07 B.E. Stetzler  Capitalized all variables for
*                                  Fortran 77 compatibility
*  2010 October  20 B.E. Stetzler  Input T corrected to be number of
*                                  centuries since J2000
*-----------------------------------------------------------------------

      IMPLICIT NONE
      INTEGER I, J
      DOUBLE PRECISION XSTA(3), XCORSTA(3), DATDI(9,5), DEG2RAD,
     .                 T, S, PR, H, P, ZNS, PS, RSTA, SINPHI,
     .                 COSPHI, COSLA, SINLA, THETAF, DR, DN, DE,
     .                 DR_TOT, DN_TOT
      DOUBLE PRECISION D2PI
      PARAMETER ( D2PI = 6.283185307179586476925287D0 )

      DATA ((DATDI(I,J),I=1,9),J=1,5)/
     .   0, 0, 0, 1, 0,   0.47D0, 0.23D0, 0.16D0, 0.07D0,
     .   0, 2, 0, 0, 0,  -0.20D0,-0.12D0,-0.11D0,-0.05D0,
     .   1, 0,-1, 0, 0,  -0.11D0,-0.08D0,-0.09D0,-0.04D0,
     .   2, 0, 0, 0, 0,  -0.13D0,-0.11D0,-0.15D0,-0.07D0,
     .   2, 0, 0, 1, 0,  -0.05D0,-0.05D0,-0.06D0,-0.03D0/

      DEG2RAD = D2PI/360D0

*  Compute the phase angles in degrees.
      S = 218.31664563D0
     .  + (481267.88194D0
     .  + (-0.0014663889D0
     .  + (0.00000185139D0)*T)*T)*T

      PR = (1.396971278D0
     .   + (0.000308889D0
     .   + (0.000000021D0
     .   + (0.000000007D0)*T)*T)*T)*T

      S = S + PR

      H = 280.46645D0
     .  + (36000.7697489D0
     .  + (0.00030322222D0
     .  + (0.000000020D0
     .  + (-0.00000000654D0)*T)*T)*T)*T

      P = 83.35324312D0
     .  + (4069.01363525D0
     .  + (-0.01032172222D0
     .  + (-0.0000124991D0
     .  + (0.00000005263D0)*T)*T)*T)*T

      ZNS = 234.95544499D0
     .    + (1934.13626197D0
     .    + (-0.00207561111D0
     .    + (-0.00000213944D0
     .    + (0.00000001650D0)*T)*T)*T)*T

      PS = 282.93734098D0
     .   + (1.71945766667D0
     .   + (0.00045688889D0
     .   + (-0.00000001778D0
     .   + (-0.00000000334D0)*T)*T)*T)*T

      RSTA=DSQRT(XSTA(1)**2+XSTA(2)**2+XSTA(3)**2)
      SINPHI=XSTA(3)/RSTA
      COSPHI=DSQRT(XSTA(1)**2+XSTA(2)**2)/RSTA

      COSLA=XSTA(1)/COSPHI/RSTA
      SINLA=XSTA(2)/COSPHI/RSTA

* Reduce angles to between the range 0 and 360.
      S =  DMOD(S,360D0)
*      TAU = DMOD(TAU,360D0)
      H =  DMOD(H,360D0)
      P =  DMOD(P,360D0)
      ZNS = DMOD(ZNS,360D0)
      PS = DMOD(PS,360D0)

      DR_TOT = 0D0
      DN_TOT = 0D0

      DO 99 I=1,3
      XCORSTA(I)=0D0
99    CONTINUE
      DO 98 J=1,5
      THETAF=(DATDI(1,J)*S+DATDI(2,J)*H+DATDI(3,J)*P+
     . DATDI(4,J)*ZNS+DATDI(5,J)*PS)*DEG2RAD

      DR=DATDI(6,J)*(3D0*SINPHI**2-1D0)/2D0*COS(THETAF)+
     .     DATDI(8,J)*(3D0*SINPHI**2-1D0)/2D0*SIN(THETAF)

      DN=DATDI(7,J)*(COSPHI*SINPHI*2D0)*COS(THETAF)+
     .     DATDI(9,J)*(COSPHI*SINPHI*2D0)*SIN(THETAF)

      DE = 0D0
      DR_TOT = DR_TOT+DR
      DN_TOT = DN_TOT+DN

      XCORSTA(1)=XCORSTA(1)+DR*COSLA*COSPHI-DE*SINLA
     . -DN*SINPHI*COSLA
      XCORSTA(2)=XCORSTA(2)+DR*SINLA*COSPHI+DE*COSLA
     . -DN*SINPHI*SINLA
      XCORSTA(3)=XCORSTA(3)+DR*SINPHI+DN*COSPHI
98    CONTINUE

      DO 97 I=1,3
      XCORSTA(I)=XCORSTA(I)/1000D0
97    CONTINUE
      RETURN

*  Finished.

*+----------------------------------------------------------------------
*
*  Copyright (C) 2008
*  IERS Conventions Center
*
*  ==================================
*  IERS Conventions Software License
*  ==================================
*
*  NOTICE TO USER:
*
*  BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING TERMS AND CONDITIONS
*  WHICH APPLY TO ITS USE.
*
*  1. The Software is provided by the IERS Conventions Center ("the
*     Center").
*
*  2. Permission is granted to anyone to use the Software for any
*     purpose, including commercial applications, free of charge,
*     subject to the conditions and restrictions listed below.
*
*  3. You (the user) may adapt the Software and its algorithms for your
*     own purposes and you may distribute the resulting "derived work"
*     to others, provided that the derived work complies with the
*     following requirements:
*
*     a) Your work shall be clearly identified so that it cannot be
*        mistaken for IERS Conventions software and that it has been
*        neither distributed by nor endorsed by the Center.
*
*     b) Your work (including source code) must contain descriptions of
*        how the derived work is based upon and/or differs from the
*        original Software.
*
*     c) The name(s) of all modified routine(s) that you distribute
*        shall be changed.
*
*     d) The origin of the IERS Conventions components of your derived
*        work must not be misrepresented; you must not claim that you
*        wrote the original Software.
*
*     e) The source code must be included for all routine(s) that you
*        distribute.  This notice must be reproduced intact in any
*        source distribution.
*
*  4. In any published work produced by the user and which includes
*     results achieved by using the Software, you shall acknowledge
*     that the Software was used in obtaining those results.
*
*  5. The Software is provided to the user "as is" and the Center makes
*     no warranty as to its use or performance.   The Center does not
*     and cannot warrant the performance or results which the user may
*     obtain by using the Software.  The Center makes no warranties,
*     express or implied, as to non-infringement of third party rights,
*     merchantability, or fitness for any particular purpose.  In no
*     event will the Center be liable to the user for any consequential,
*     incidental, or special damages, including any lost profits or lost
*     savings, even if a Center representative has been advised of such
*     damages, or for any claim by any third party.
*
*  Correspondence concerning IERS Conventions software should be
*  addressed as follows:
*
*                     Gerard Petit
*     Internet email: gpetit[at]bipm.org
*     Postal address: IERS Conventions Center
*                     Time, frequency and gravimetry section, BIPM
*                     Pavillon de Breteuil
*                     92312 Sevres  FRANCE
*
*     or
*
*                     Brian Luzum
*     Internet email: brian.luzum[at]usno.navy.mil
*     Postal address: IERS Conventions Center
*                     Earth Orientation Department
*                     3450 Massachusetts Ave, NW
*                     Washington, DC 20392
*
*
*-----------------------------------------------------------------------
      END

max_stars_repo_name: fsoe/RTKLIB

max_stars_repo_licenses: Intel BSD-2-Clause Unlicense

max_issues_repo_name: fsoe/RTKLIB

max_issues_repo_licenses: Intel BSD-2-Clause Unlicense

max_forks_repo_name: fsoe/RTKLIB

max_forks_repo_licenses: Intel BSD-2-Clause Unlicense


text:
!This module contains all the must have physics variables. It is imported by both
!the other physics modules (so they can modify the physics) and the chemistry
!module (so it can use the physics for reaction rates).

MODULE physicscore
    USE constants
    IMPLICIT NONE
    !Use main loop counters in calculations so they're kept here
    INTEGER :: dstep,points
    !Switches for processes are also here, 1 is on/0 is off.
    LOGICAL :: endAtFinalDensity,freefall,instantSublimation

    !evap changes evaporation mode (see chem_evaporate), ion sets c/cx ratio (see initializeChemistry)
    !Flags let physics module control when evap takes place.flag=0/1/2 corresponding to not yet/evaporate/done
    INTEGER :: ion

    !Optional CR attentuation with column density and better H2 dissociation rates
    LOGICAL ::cosmicRayAttenuation,improvedH2CRPDissociation
    REAL(dp) :: h2CRPRate,zetaScale
    CHARACTER(len=1) :: ionModel="L"

    !variables either controlled by physics or that user may wish to change
    REAL(dp) :: initialDens,timeInYears,targetTime,currentTime,currentTimeold,finalDens,finalTime,initialTemp
    REAL(dp) ::  freefallFactor,cloudSize,rout,rin,baseAv,zeta
    REAL(dp), allocatable :: av(:),coldens(:),gasTemp(:),dustTemp(:),density(:)

    !Arrays fopr calculating rates
    !ckLIon and ckHIon are the coefficients for the cosmic ray ionization rate
    !ckLDiss and ckHDiss are the coefficients for the H2 dissociation rate
    !if ionModel = L use the L model coefficients, if = H use the H model
    REAL(dp),PARAMETER :: ckLIon(10)=(/1.545456645800d7, -6.307708626617d6, 1.142680666041d6, -1.205932302621d5,&
    8.170913352693d3, -3.686121296079d2,1.107203722057d1, -2.135293914267d-1,&
    2.399219033781d-3, -1.196664901916d-5/)
    REAL(dp), PARAMETER :: ckLDiss(10)=(/1.582911005330d7,-6.465722684896d6, 1.172189025424d6, -1.237950798073d5, &
        8.393404654312d3, -3.788811358130d2, 1.138688455029d1, -2.197136304567d-1, &
        2.469841278950d-3, -1.232393620924d-5/)
    REAL(dp),PARAMETER :: ckHIon(10)=(/1.223529865309d7, -5.013766644305d6, 9.120125566763d5, -9.665446168847d4,&
        6.576930812109d3, -2.979875686226d2,8.989721355058d0, -1.741300519598d-1,&
        1.965098116126d-3, -9.844203439473d-6/)
    REAL(dp), PARAMETER :: ckHDiss(10)=(/1.217227462831d7,-4.989649250304d6, 9.079152156645d5, -9.624890825395d4, &
        6.551161486120d3, -2.968976216187d2, 8.959037875226d0, -1.735757324445d-1, &
        1.959267277734d-3, -9.816996707980d-6/)

CONTAINS
    !basic initialization of physics. All physics modules should call this and then
    !do their own initialization.
    SUBROUTINE coreInitializePhysics(successFlag)
        INTEGER, INTENT(OUT) :: successFlag
        timeInYears=currentTime/SECONDS_PER_YEAR

        ! Modules not restarted in python wraps so best to reset everything manually.
        IF (ALLOCATED(av)) DEALLOCATE(av,coldens,gasTemp,dustTemp,density)
        ALLOCATE(av(points),coldens(points),gasTemp(points),dustTemp(points),density(points))

        cloudSize = (rout-rin)*pc
        gasTemp=initialTemp
        dustTemp=gasTemp
        density=initialDens
        currentTimeOld=0.0
        IF (.not. ((ionModel .eq. 'L') .or. (ionModel .eq. 'H'))) THEN
            successFlag=-1
            write(*,*) "Error: ionModel must be either L or H"
            RETURN
        END IF
        IF ((improvedH2CRPDissociation) .and. (.not. cosmicRayAttenuation)) THEN
            successFlag=-1
            write(*,*) "Error: improvedH2CRPDissociation requires cosmicRayAttentuation to also be True"
            RETURN
        END IF

        !calculate initial column density as distance from core edge to current point * density
        DO dstep=1,points
            coldens(dstep)=real(points-dstep+1)*cloudSize/real(points)*initialDens
        END DO
          !calculate the Av using an assumed extinction outside of core (baseAv), depth of point and density
        av= baseAv +coldens/1.6d21
        zetaScale=zeta
    END SUBROUTINE coreInitializePhysics

    SUBROUTINE coreUpdatePhysics
        !calculate column density. Remember dstep counts from core centre to edge
        !and coldens should be amount of gas from edge to parcel.
        coldens(dstep)=cloudSize/real(points)*density(dstep)

        ! add previous column densities to current as we move into cloud to get total
        IF (dstep .lt. points) coldens(dstep)=coldens(dstep)+coldens(dstep-1)

        !calculate the Av using an assumed extinction outside of core (baseAv), depth of point and density
        av(dstep)= baseAv +coldens(dstep)/1.6d21
        dustTemp=gasTemp

        IF (cosmicRayAttenuation) CALL ionizationDependency
    END SUBROUTINE coreUpdatePhysics

    pure FUNCTION densdot(density)
    ! Returns the time derivative of the density.
    ! Analytical function taken from Rawlings et al. 1992
    ! Called from chemistry.f90, density integrated with abundances so this gives ydot
    REAL(dp), INTENT(IN) :: density
    REAL(dp) :: densdot
    !Rawlings et al. 1992 freefall collapse. With freefallFactor for B-field etc
    IF ((density .lt. finalDens) .and. (freefall)) THEN
        densdot=freefallFactor*(density**4./initialDens)**0.33*&
        &(8.4d-30*initialDens*((density/initialDens)**0.33-1.))**0.5
    ELSE
        densdot=0.0
    ENDIF
    END FUNCTION densdot

    pure FUNCTION dByDnDensdot(density)
    !Defunct function which provides the necessary derivative d(dn/dt)/dn
    !in the case one uses a Jacobian.
    REAL(dp), INTENT(IN) :: density
    REAL(dp) :: dByDnDensdot
    !Rawlings et al. 1992 freefall collapse. With freefallFactor for B-field etc
    IF (density .lt. finalDens) THEN
        dByDndensdot=freefallFactor*8.4d-30*(density**3)*((9.0d0*((density/initialDens)**0.33))-8.0d0)
        dByDnDensdot=dByDnDensdot/(6.0d0*(((density**4.0)/initialDens)**0.66))
        dByDnDensdot=dByDnDensdot/dsqrt(initialDens*8.4d-30*(((density/initialDens)**0.33))-1.0d0)
    ELSE
        dByDnDensdot=0.0
    ENDIF
    END FUNCTION dByDnDensdot

    SUBROUTINE ionizationDependency
        REAL(dp) :: dissSum,dRate,zSum,ionRate
        INTEGER :: k
        !Attenuate CR by column density
        zeta= 1.0
        zSum = 0
        DO k=0,9,1
            IF (ionModel .eq. 'L') THEN
                ionRate=ckLIon(k+1)*log10(coldens(dstep))**k
            ELSEIF (ionModel .eq. 'H') THEN
                ionRate=ckHIon(k+1)*log10(coldens(dstep))**k
            ELSE
                write(*,*) "WARNING: ionModel switch must be 0 or 1"
            ENDIF
            zSum=zSum+ionRate
        END DO

        zeta=((10**zSum)/1.3d-17)* zetaScale

        !rate calculation for H2 dissociation
        IF (improvedH2CRPDissociation) THEN
            dissSum = 0
            DO k=0,9,1
                IF (ionModel .eq. 'L') THEN
                    dRate=ckLDiss(k+1)*log10(coldens(dstep))**k
                ELSEIF (ionModel .eq. 'H') THEN
                    dRate=ckHDiss(k+1)*log10(coldens(dstep))**k
                ELSE
                    write(*,*) "WARNING: ionModel switch must be L or H"
                ENDIF
                dissSum=dissSum+dRate
            END DO
            h2CRPRate=(10**dissSum)*zetaScale
        END IF
    END SUBROUTINE ionizationDependency

END MODULE physicscore

max_stars_repo_name: jonholdship/UCLCHEM

max_issues_repo_name: jonholdship/UCLCHEM

max_forks_repo_name: jonholdship/UCLCHEM


text:
! { dg-do compile }
! { dg-options "-fcoarray=single -fdump-tree-original" }
!
! Coarray support -- allocatable array coarrays
! PR fortran/18918
!
integer,allocatable :: a(:)[:,:]
nn = 5
mm = 7
allocate(a(nn)[mm,*])
end

subroutine testAlloc3
  implicit none
  integer, allocatable :: ab(:,:,:)[:,:]
  integer, allocatable, dimension(:),codimension[:] :: b(:,:,:)[:,:]
  integer, allocatable, dimension(:,:),codimension[:,:,:] :: c
  integer, allocatable, dimension(:,:),codimension[:,:,:] :: d[:,:]
  integer, allocatable, dimension(:,:,:),codimension[:,:,:] :: e(:,:)
  integer, allocatable, dimension(:,:,:),codimension[:,:,:] :: f(:,:)[:,:]

  allocate(ab(1,2,3)[4,*])
  allocate(b(1,2,3)[4,*])
  allocate(c(1,2)[3,4,*])
  allocate(d(1,2)[3,*])
  allocate(e(1,2)[3,4,*])
  allocate(f(1,2)[3,*])
end subroutine testAlloc3

subroutine testAlloc4()
  implicit none
  integer, allocatable :: xxx(:)[:,:,:,:]
  integer :: mmm
  mmm=88
  allocate(xxx(1)[7,-5:8,mmm:2,*])
end subroutine testAlloc4

subroutine testAlloc5()
  implicit none
  integer, allocatable :: yyy(:)[:,:,:,:]
  integer :: ooo, ppp
  ooo=88
  ppp=42
  allocate(yyy(1)[7,-5:ppp,1,ooo:*])
end subroutine testAlloc5

! { dg-final { scan-tree-dump-times "a.dim.0..lbound = 1;"     1 "original" } }
! { dg-final { scan-tree-dump-times "a.dim.0..ubound = .*nn;" 1 "original" } }
! { dg-final { scan-tree-dump-times "a.dim.1..lbound = 1;"     1 "original" } }
! { dg-final { scan-tree-dump-times "a.dim.1..ubound = .*mm;" 1 "original" } }
! { dg-final { scan-tree-dump-times "a.dim.2..lbound = 1;"     1 "original" } }
! { dg-final { scan-tree-dump-times "a.dim.2..ubound"          0 "original" } }

! { dg-final { scan-tree-dump-times "xxx.dim.0..lbound = 1;"     1 "original" } }
! { dg-final { scan-tree-dump-times "xxx.dim.0..ubound = 1;"     1 "original" } }
! { dg-final { scan-tree-dump-times "xxx.dim.1..lbound = 1;"     1 "original" } }
! { dg-final { scan-tree-dump-times "xxx.dim.1..ubound = 7;"     1 "original" } }
! { dg-final { scan-tree-dump-times "xxx.dim.2..lbound = -5;"    1 "original" } }
! { dg-final { scan-tree-dump-times "xxx.dim.2..ubound = 8;"     1 "original" } }
! { dg-final { scan-tree-dump-times "xxx.dim.3..lbound = .*mmm;" 1 "original" } }
! { dg-final { scan-tree-dump-times "xxx.dim.3..ubound = 2;"     1 "original" } }
! { dg-final { scan-tree-dump-times "xxx.dim.4..lbound = 1;"     1 "original" } }
! { dg-final { scan-tree-dump-times "xxx.dim.4..ubound"          0 "original" } }

! { dg-final { scan-tree-dump-times "yyy.dim.0..lbound = 1;"     1 "original" } }
! { dg-final { scan-tree-dump-times "yyy.dim.0..ubound = 1;"     1 "original" } }
! { dg-final { scan-tree-dump-times "yyy.dim.1..lbound = 1;"     1 "original" } }
! { dg-final { scan-tree-dump-times "yyy.dim.1..ubound = 7;"     1 "original" } }
! { dg-final { scan-tree-dump-times "yyy.dim.2..lbound = -5;"    1 "original" } }
! { dg-final { scan-tree-dump-times "yyy.dim.2..ubound = .*ppp;" 1 "original" } }
! { dg-final { scan-tree-dump-times "yyy.dim.3..lbound = 1;"     1 "original" } }
! { dg-final { scan-tree-dump-times "yyy.dim.3..ubound = 1;"     1 "original" } }
! { dg-final { scan-tree-dump-times "yyy.dim.4..lbound = .*ooo;" 1 "original" } }
! { dg-final { scan-tree-dump-times "yyy.dim.4..ubound"          0 "original" } }

! { dg-final { cleanup-tree-dump "original" } }

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*	PROGRAM D_CGRADIENT
*
* Central reconstruction with conjugate gradients
*
* INPUT:
*	*.PHIT  First approx. of the object
*	*.PCR   First approx. of the object: real part of FT
*	*.PCI   First approx. of the object: imag. part of FT
*	*.V  	Support of the object (direct space)
*	*.G	Regularizing function
*
* OUTPUT:
*	*.PHIN	Final estimation of the object
*
* Log file: *.SUM
*
*
* Pointers :
*  ---------
* p_Phit	PhiT		partie reelle			R4 -
* p_PhitCr	PhiT chapeau	partie reelle			R4 -
* p_PhitCi	  "     "	partie imaginaire		R4 -
* p_V		support de l'Objet				R4 -
* p_G		Fonction de regularisation			R4 -
* p_phin	Objet reconstruit (module)			R4 -
* p_Zn		|	Espace de travail pour les fonctions	R4   vm
* p_Rn		|	zn et rn (Cf. These page 42)		R4   vm
* p_ddnr	|	et pour la fonction dn		R8	     vm
* p_ddnI	|					R8	     vm
* p_ddntemp	|	tableau temporaire pour dn-1	R8	     vm
* p_dphitr	PhiT		partie reelle		R8	     vm
* p_dphiti	PhiT		partie imaginaire	R8	     vm
*
* From cgradient.for, S. ROQUES - J. VIGNEAU, Version 01-12-89
*
* JLP
* Version 02-06-92
*-------------------------------------------------------------------------------
	PROGRAM D_CGRADIENT
	implicit none
	integer*4 npix(2),n,n4,n8,p_phit,p_phitcr, p_phitci
	integer*4 p_v,p_g,p_zn,p_rn,p_phin,p_ddntemp,istatus
	integer*4 p_ddnr,p_ddni,p_dphitr,p_dphiti,long,lu,i
	real*4 eps(3)
	character name*20,comments*80,filesum*40,buffer*90
 	character*40 fileout,filephit,filepcr,filepci,filev,fileg
	integer*4 madrid(1)
	common /vmr/madrid
*
* Initialisation
	call jlp_begin
        write(6,112)
112     format(/,'******* Program d_cgradient  Version 12-06-92')
	call jlp_inquifmt

*
* ouverture du fichier de la recapitulation
55	write(6,*)' Generic name (should end with a dot) ?'
	read(5,10) name
10	format(a)
	long=index(name,'.')-1
	if(long.le.0)goto 55
*
	filesum=name(:long)//'.SUM'
	lu=1
C "APPEND" is not available for IBM...
C	open(lu,file=filesum,status='OLD',access='APPEND')
	open(lu,file=filesum,status='OLD',access='SEQUENTIAL')
        do i=1,10000
        read(lu,10,end=989) buffer
        enddo
989	write (lu,980)
980	format(/,' Program CGRADIENT',/,
     1	' ----------------',/)
*
* Ouverture du fichier contenant la premiere approximation de l'Objet
*
	FilePhit = name(:long)//'.PHIT'
	write(6,*)' Reading the first approx. of the object:',filephit
	call jlp_vm_readimag(p_phit,npix(1),npix(2),filephit,comments)
* Recenter image to get gravity center at (1,1):
        call recentre(madrid(p_phit),madrid(p_phit),
     1     npix(1),npix(2),npix(1))
*
	N = npix(1)*npix(2)
*
* reservons de la memoire et passons en D.P. pour les parties R et I
*
        N8 = N*8
	call jlp_getvm (p_dphitr,N8)
	call jlp_getvm (p_dphiti,N8)
*
	call S_to_D (madrid(p_phit), madrid(p_dphitr), N)
*
* Ouverture des fichiers contenant la TF de la premiere approximation de
* l'Objet (partie reelle et partie imaginaire)
*
	FilePcr = name(:long)//'.PCR'
	write(6,*)
     1  ' Reading Fourier transf. (real) of the first approx. of the object',
     1	' (.PCR) : ',filepcr
	call jlp_vm_readimag(p_phitcr,npix(1),npix(2),filepcr,comments)
* Recenter image to get zero frequency at (1,1):
        call recentre(madrid(p_phitcr),madrid(p_phitcr),
     1   npix(1),npix(2),npix(1))
*
	FilePci = name(:long)//'.PCI'
	write(6,*)
     1  ' Reading Fourier transf. (real) of the first approx. of the object',
     1	' (.PCI) : ',filepci
	call jlp_vm_readimag(p_phitci,npix(1),npix(2),filepci,comments)
* Recenter image to get zero frequency at (1,1):
        call recentre(madrid(p_phitci),madrid(p_phitci),
     1     npix(1),npix(2),npix(1))
*
* Ouverture du fichier contenant le support de l'objet
*
	FileV = name(:long)//'.V'
	write(6,*)' Reading the support (.V) of the object',filev
	call jlp_vm_readimag(p_v,npix(1),npix(2),filev,comments)
* Recenter image to get gravity center at (1,1):
        call recentre(madrid(p_v),madrid(p_v),npix(1),npix(2),npix(1))
*
* Ouverture du fichier contenant la fonction de regularisation
*
	FileG = name(:long)//'.G'
	write(6,*)' Reading regul. function (.G) ',fileg
	call jlp_vm_readimag(p_g,npix(1),npix(2),fileg,comments)
* Recenter image to get zero frequency at (1,1):
        call recentre(madrid(p_g),madrid(p_g),npix(1),npix(2),npix(1))
*
* et lecture du descripteur Epsilon
* (doesn't work in current version 02-05-91)
C	call jlp_get_symbol('D_EPSILON',buffer,istatus)
C	read(buffer,*,err=87) (eps(i),i=1,3)
C	if(istatus.eq.0) goto 88
87	write(6,*) ' Epsilon1, eps2, eps3 (from "erreur")?'
	read(5,*) eps(1),eps(2),eps(3)
*
88	write(6,*) ' Epsilon1, eps2, eps3:'
	write(6,*) eps(1),eps(2),eps(3)
*
* reservons de la place pour travailler...
*
        N4 = N*4
        N8 = N*8
	call jlp_getvm (p_rn,N4)
	call jlp_getvm (p_zn,N8)
 	call jlp_getvm (p_ddntemp,N8)
	call jlp_getvm (p_ddnr,N8)
	call jlp_getvm (p_ddni,N8)
*
* Ouverture du fichier contenant le resultat final
	call jlp_getvm(p_phin,N4)
*
	call restit(n,npix,madrid(p_phitcr),madrid(p_phitci),
     1	madrid(p_phit),madrid(p_dphitr),madrid(p_dphiti),
     1	madrid(p_v),madrid(p_g),madrid(p_rn),madrid(p_zn),
     1	madrid(p_ddnr),madrid(p_ddni),madrid(p_ddntemp),
     1	madrid(p_phin),p_dphitr,p_dphiti,
     1	p_ddnr,p_ddni,eps,lu)
*
* et on ferme les fichiers
*
	fileout = name(:long)//'.PHIN'
	comments = 'Resultat final'
* Recenter image to get gravity center in the middle:
        call recentre(madrid(p_phin),madrid(p_phin),npix(1),npix(2),npix(1))
	call jlp_writeimag(madrid(p_phin),npix(1),npix(2),npix(1),
     1	fileout,comments)
*
* Sortie normale
	close(lu)
	call jlp_end
	stop
	end
*------------------------------------------------------------------------------
*
	subroutine restit (n,npix,phitcr,phitci,phit,
     1	dphitr, dphiti,v,g,rn,zn,ddnr,ddni,ddntemp,
     1	phin,p_dphitr, p_dphiti,p_ddnr, p_ddnI, eps, lu)
*
*             *********************************
*
*                  RECONSTRUCTION CENTRALE
*
*             *********************************
*
*                  VALEURS PROPRES DE A*A
*
*   Application de la methode des gradients conjugues
*   Calcul du polynome dont les racines convergent vers les valeurs propres
*
*             *********************************
*
	implicit none
*
       	integer*4  N,npix(2),lu,p_dphitr, p_dphiti,p_ddnr, p_ddnI
	integer*4 k,i,m,mult(35),icn(101),itest,isig,num,nomb,nrac
	real*4	V(N),G(N),eps (3),Rn(N),phin(N),PhitcR(N), PhitcI(N),
     1	Phit(N)
 	real*4	rac(35),xabs(101),xrn(101), xdn(101),
     1	xnorfin,xx, x,xnun,pas,xnfinc,deltan, ron,
     1	tete,tetap,xmum,major
	real*8	ddnr(N), ddnI(N),dphitr(N), dphiti(N),Zn(N),ddntemp(N)
	real*8	DG2,xnorrn,psdz,omegan
        integer*4 iter_max
*
* pour eviter la limitation a .29D-38
C JLP91	real*16	test,y,xnorrnp
	real*8	test,y,xnorrnp
	integer*4 madrid(1)
	common /vmr/madrid
*
* majorant de 1/mu
	major=10.
* Max number of iterations:
        iter_max=20
*
***************************************************************************
*
*
* N = npix ** 2       nombre total de pixels
* (PhitcR,PhitcI)     tf de phit - parties reelle et imaginaire
* (Phit)	     module de phit
* p_dphitr et p_dphiti  pointeurs associes en dp
* V		   support de l'objet
* G		   fonction de regularisation
* eps		 premiere estimation d'erreur (epsi)
*
***************************************************************************
*
* initialisation polynomes
* xabs = abscisses  -  xrn = polynome rn  -  xdn = polynome dn
* icn = polynome cn des changements de signe
	do k = 1, 101
	 xabs (k) = (k-1) * 0.01
	 xrn (k) = 1.
	 xdn (k) = 1.
	 icn (k) = 0
	end do
*
* initialisation gradients conjugues
* fonction de depart phi zero = v phit
	do k = 1, N
*SP
	 phin(k)=V(k)*Phit(k)
*dp
	 dphitr(k)=dble(phin(k))
	 dphiti(k)=0.
	end do
*
* calcul de r zero = v U* g2 (phitc - phi zero chap)
* calculs de d zero et de ||rn||**2 pour n=0
*
* p_dphitr : pointeur associe a dphitr
* p_dphiti : pointeur associe a dphiti
*
	call fourier1(madrid(p_dphitr),madrid(p_dphiti),npix(1),npix(2),1)
C JLP92
C	print *,' phitr(4),phiti(4)',dphitr(4),dphiti(4)
*
	do k = 1, N
	 DG2 = dble(G(k)) * dble(G(k))
	 dphitr(k) = (dble(phitcr(k)) - dphitr(k)) * dg2
	 dphiti(k) = (dble(phitci(k)) - dphiti(k)) * dg2
	end do
*
	call fourier1(madrid(p_dphitr),madrid(p_dphiti),npix(1),npix(2),-1)
C JLP92
C	print *,' phitr(4),phiti(4)',dphitr(4),dphiti(4)
*
	call mod_tf(n,madrid(p_dphitr),madrid(p_dphiti))
*
	xnorrn = 0.
	do k = 1, N
	 dphiti(k) = 0.
	 rn(k) = sngl(dphitr(k)) * V(k)
	 ddnr(k) = dble(rn(k))
	 xnorrn = xnorrn + ddnr(k) * ddnr(k)
	end do
*
	write(6,105)
	write(lu,105)
105	format('  N         Test',/)
*
*******************************************************************************
*
*			I T E R A T I O N
*
*******************************************************************************
*
	I = 0
600	I = I + 1
*
* on garde l'ancien ddnr dans ddntemp
	do k = 1, N
 	 ddntemp(k) = ddnr(k)
	 ddnI(k) = 0.
	end do
*
* calcul de dn chapeau
	call fourier1(madrid(p_ddnr),madrid(p_ddni),npix(1),npix(2),1)
*
* calcul de zn et de omega n
*
	do k = 1, N
	 DG2 = dble(G(k)) * dble(G(k))
	 ddnr(k) = ddnr(k) * DG2
	 ddnI(k) = ddnI(k) * DG2
	end do
*
	call fourier1(madrid(p_ddnr),madrid(p_ddni),npix(1),npix(2),-1)
* Scalar product: (d_n | A*A d_n)
	psdz = 0.
*
	do k = 1, N
	 Zn(k) = ddnr(k) * V(k)
	 psdz = psdz + (ddntemp(k)* 1.e10) * (Zn(k) * 1.e10)
	end do
*
	if (psdz .ne. 0.) then
	 omegan = xnorrn * 1.e20/psdz
	else
	 write(6,109)
109	 format(' d_cgradient/Fatal error when computing Omega_n'/,
     1   ' Emergency exit *** (d_n | A*A d_n) = 0')
	 write(lu,103)
	 stop
	end if
*
* calculs de phi n+1 et de ||phi n+1||**2
* calculs de r n+1 et de ||r n+1||**2
*
	xnorfin = 0.
	xnorrnp = 0.
*
	do k = 1, N
	 phin(k) = phin(k) + omegan * ddntemp(k)
	 x = phin(k) * phin(k)
	 xnorfin = xnorfin + x
	 Rn(k) = Rn(k) - sngl(omegan * Zn(k))
C JLP91	 y = qext(Rn(k)) * qext(Rn(k))
	 y = dble(Rn(k)) * dble(Rn(k))
	 xnorrnp = xnorrnp + y
	end do
*
* formation des polynomes Rn et Cn
* nombre de degres de liberte : icn(95)
* isig correspond au polynome Sn(X)
*
	itest = 0
	do k = 1, 101
	 isig = sign (1.,xrn(k))
	 xx = sngl(omegan) * xabs(k) * xdn(k)
	 xrn(k) = xrn(k) - xx
*
* ... on evite le point d'accumulation en 1 :
*
	 if (k .le. 95) then
	  if(isig .ne. sign(1.,xrn(k))) then
	   itest = 1
	   icn(k) = icn(k) + 1
	  end if
	 end if
	end do
*
* test d'arret
*
	test = xnorrnp/xnorfin
	test = major * sqrt(test)
*
	write (6,102) i, test
	write (lu,102) i, test
102	format(1x,i2,5x,g12.6)
*
	if (test.le.1.e-03.or.i.ge.iter_max) then
*
* ... itest=0 signifie que Rn n'a change de signe en aucun point xk
*
	 if (itest.eq.0.or.i.ge.iter_max) goto 500
	end if
*
* calcul de nu n
	xnun = sngl(xnorrnp/xnorrn)
*
* calcul de d n+1
	do k = 1, N
	 ddnr(k) = dble(Rn(k)) + dble(xnun) * ddntemp(k)
	end do
*
* calcul du polynome D n+1
	do k = 1, 101
	 xdn(k) = xrn(k) + xnun * xdn(k)
	end do
*
* reinitialisations
	xnorrn = dble(xnorrnp)
*
* retour a I = I + 1
*
	if(i.lt.iter_max) goto 600
*
* ********************  FIN D'ITERATION   ********************
*
500	continue
*
	write(6,1000) icn(95)
	write(lu,1000) icn(95)
1000	format (/,' Number of degrees of freedom : ', i3,/)
*
* calcul des valeurs propres de A*A
*
* rac est la racine de rn
* num est son numero d'ordre
* mult est sa multiplicite
* nomb est le nombre de racines (y compris la multiplicite)
*
	pas = 1./(2.*100.)
	m = 94
	num = 0
	nomb = 0
*
	write(6,106)
	write(lu,106)
106	format('  N   Eigenvalues       Multiplicity',/)
	do k = 1, m
	 nrac = icn(k+1) - icn(k)
	 if (nrac .ne. 0) then
	  num = num + 1
	  rac (num) = xabs(k+1) - pas
	  mult (num) = nrac
	  nomb = nomb + nrac
	  write (lu,103) num,rac(num),mult(num)
	  write (6,103) num,rac(num),mult(num)
103	  format(1x, i2, 10x, f6.4, 10x, i2)
	 end if
	end do
*
* calcul des majorants de l'erreur
*
	do k = 1, N
	 dphitr(k) = dble(G(k)) * (dphitr(k) - phin(k))
	end do
*
	call dnormel2 (N, dphitr, dphiti, deltan)
	xnfinc = 0.
	do k = 1, N
	 xnfinc = xnfinc + phin(k) * phin(k)
	end do
	xnfinc = sqrt (xnfinc)
	ron = sqrt( abs(eps(1) - deltan*deltan))
*
* borne superieure de l'erreur
*
	tete = sqrt(rac(1) * xnfinc)
	tete = ron / tete
*
* calcul du 1/mu moyen
* calcul des ereurs correspondantes
*
	xmum = 0.
	do k = 1, num
	 xmum = xmum + (1./rac(k))*mult(k)
	end do
*
	xmum = xmum/float(nomb)
	xmum = (xmum + 1.)/2.
	xmum = sqrt(xmum)
	tetap = (xmum*ron)/xnfinc
*
	write (6, 1001) xnfinc,xmum,sqrt(eps(1))
	write (lu, 1001) xnfinc,xmum,sqrt(eps(1))
1001	format(/,' Norm of the solution Phi n  ( ||Phi n|| ) :',
     1	'. . . . ',g12.5,/,
     1	' Weighted mean of the inverses of the eigenvalues (1/mu): ',
     1	g12.5,/,' Error of image reconstruction (sqrt(eps1)) :',
     1	' . . . . ',g12.5)
*
	write (6, 1002) sqrt(eps(3)),sqrt(eps(2)),deltan,tetap
	write (lu, 1002) sqrt(eps(3)),sqrt(eps(2)),deltan,tetap
1002	format(' Error due to the choice of Hr: ',
     1	'(epsilon_o) sqrt(eps3) : ',g12.5,/,
	1 	' Erreur related to the SNR ratio: ',
     1	'(epsilon_i) sqrt(eps2): ',g12.5,/,
     1	' Term to correct the error (Delta_n) :',
     1	' . . . .  ',g12.5,/,
	1 	' Relative error for object reconstruction:',
     1	'  . . . .  ',g12.5)
*
	return
	end
C-----------------------------------------------------------------------
	subroutine dnormel2(n,reel,imagi,norme)
*
* Calcul de la norme L2 du complexe (REEL, IMAG) dans l'espace de Fourier
*
* NE DOIT PAS ETRE APPELE EN UTILISANT LE MECANISME	%VAL
* ----------------------------------------------------------
*
	implicit none
	integer*4 n,i
	real*8	reel(*), imagi(*)
	real*4	norme
*
	norme = 0.
	do i = 1, n
	 norme = norme + (reel(i)*reel(i)+imagi(i)*imagi(i))
	end do
*
	norme = sqrt(norme/float(n))
*
	return
	end
C-----------------------------------------------------------------------------
C	include 'd_utilities.for'

ext: for

max_stars_repo_name: jlprieur/deconv

max_issues_repo_name: jlprieur/deconv

max_forks_repo_name: jlprieur/deconv


  MODULE maxmin_module

  CONTAINS

      subroutine maxmin(izz,jzz,kzz,f,nstat,rstat,kmin,kmax,amax,amin)
      use input
      use mpi
      implicit none

      integer :: izz,jzz,kzz,nstat,kmin,kmax
      real, dimension(stat_out) :: rstat
      real, dimension(1-ngxy:izz+ngxy,1-ngxy:jzz+ngxy,1-ngz:kzz+ngz) :: f
      character(len=6) :: amax,amin

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

      integer :: i,j,k
      integer :: imax,jmax,imin,jmin
      integer, dimension(nk+1) :: imaxt,jmaxt,kmaxt,imint,jmint,kmint
      real, dimension(nk+1) :: tmax,tmin
      real :: fmax,fmin,rmax,rmin
      integer :: loc
      real, dimension(2) :: mmax,nmax,mmin,nmin

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

      imin = 1
      jmin = 1
      kmin = 1
      imax = 1
      jmax = 1
      kmax = 1

!$omp parallel do default(shared)    &
!$omp private(i,j,k)
      do k=1,kzz
        tmax(k)= -1.e30
        tmin(k)=  1.e30
        do j=1,jzz
        do i=1,izz
          if(f(i,j,k).gt.tmax(k))then
            tmax(k)=f(i,j,k)
            imaxt(k)=i
            jmaxt(k)=j
            kmaxt(k)=k
          endif
          if(f(i,j,k).lt.tmin(k))then
            tmin(k)=f(i,j,k)
            imint(k)=i
            jmint(k)=j
            kmint(k)=k
          endif
        enddo
        enddo
      enddo

      fmax= -1.e30
      fmin=  1.e30
      do k=1,kzz
        if(tmax(k).gt.fmax)then
          fmax=tmax(k)
          imax=imaxt(k)
          jmax=jmaxt(k)
          kmax=kmaxt(k)
        endif
        if(tmin(k).lt.fmin)then
          fmin=tmin(k)
          imin=imint(k)
          jmin=jmint(k)
          kmin=kmint(k)
        endif
      enddo

      mmax(1)=fmax
      mmax(2)=myid
      call MPI_ALLREDUCE(mmax,nmax,1,MPI_2REAL,MPI_MAXLOC,   &
                         MPI_COMM_WORLD,ierr)
      loc=nint(nmax(2))
      imax=imax+(myi-1)*ni
      jmax=jmax+(myj-1)*nj
      call MPI_BCAST(imax,1,MPI_INTEGER,loc,MPI_COMM_WORLD,ierr)
      call MPI_BCAST(jmax,1,MPI_INTEGER,loc,MPI_COMM_WORLD,ierr)
      call MPI_BCAST(kmax,1,MPI_INTEGER,loc,MPI_COMM_WORLD,ierr)

      mmin(1)=fmin
      mmin(2)=myid
      call MPI_ALLREDUCE(mmin,nmin,1,MPI_2REAL,MPI_MINLOC,   &
                         MPI_COMM_WORLD,ierr)
      loc=nint(nmin(2))
      imin=imin+(myi-1)*ni
      jmin=jmin+(myj-1)*nj
      call MPI_BCAST(imin,1,MPI_INTEGER,loc,MPI_COMM_WORLD,ierr)
      call MPI_BCAST(jmin,1,MPI_INTEGER,loc,MPI_COMM_WORLD,ierr)
      call MPI_BCAST(kmin,1,MPI_INTEGER,loc,MPI_COMM_WORLD,ierr)

      fmax=nmax(1)
      fmin=nmin(1)

    if(myid.eq.0)then
      write(6,100) amax,fmax,imax,jmax,kmax,    &
                   amin,fmin,imin,jmin,kmin
100   format(2x,a6,':',1x,g13.6,i5,i5,i5,    &
             4x,a6,':',1x,g13.6,i5,i5,i5)

      nstat = nstat + 1
      rstat(nstat) = fmax
      nstat = nstat + 1
      rstat(nstat) = fmin
    endif

      if(timestats.ge.1) time_stat=time_stat+mytime()

      end subroutine maxmin

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine maxmin2d(izz,jzz,f,nstat,rstat,amax,amin)
      use input
      use mpi
      implicit none

      integer :: izz,jzz,nstat
      real, dimension(stat_out) :: rstat
      real, dimension(1-ngxy:izz+ngxy,1-ngxy:jzz+ngxy) :: f
      character(len=6) :: amax,amin

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

      integer :: i,j
      integer :: imax,jmax,imin,jmin
      integer, dimension(jzz) :: imaxt,jmaxt,imint,jmint
      real, dimension(jzz) :: tmax,tmin
      real :: fmax,fmin,rmax,rmin
      integer :: loc
      real, dimension(2) :: mmax,nmax,mmin,nmin

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

!$omp parallel do default(shared)    &
!$omp private(i,j)
      do j=1,jzz
        tmax(j)= -1.e30
        tmin(j)=  1.e30
        do i=1,izz
          if(f(i,j).gt.tmax(j))then
            tmax(j)=f(i,j)
            imaxt(j)=i
            jmaxt(j)=j
          endif
          if(f(i,j).lt.tmin(j))then
            tmin(j)=f(i,j)
            imint(j)=i
            jmint(j)=j
          endif
        enddo
      enddo

      fmax= -1.e30
      fmin=  1.e30
      do j=1,jzz
        if(tmax(j).gt.fmax)then
          fmax=tmax(j)
          imax=imaxt(j)
          jmax=jmaxt(j)
        endif
        if(tmin(j).lt.fmin)then
          fmin=tmin(j)
          imin=imint(j)
          jmin=jmint(j)
        endif
      enddo

      mmax(1)=fmax
      mmax(2)=myid
      call MPI_ALLREDUCE(mmax,nmax,1,MPI_2REAL,MPI_MAXLOC,   &
                         MPI_COMM_WORLD,ierr)
      loc=nint(nmax(2))
      imax=imax+(myi-1)*ni
      jmax=jmax+(myj-1)*nj
      call MPI_BCAST(imax,1,MPI_INTEGER,loc,MPI_COMM_WORLD,ierr)
      call MPI_BCAST(jmax,1,MPI_INTEGER,loc,MPI_COMM_WORLD,ierr)

      mmin(1)=fmin
      mmin(2)=myid
      call MPI_ALLREDUCE(mmin,nmin,1,MPI_2REAL,MPI_MINLOC,   &
                         MPI_COMM_WORLD,ierr)
      loc=nint(nmin(2))
      imin=imin+(myi-1)*ni
      jmin=jmin+(myj-1)*nj
      call MPI_BCAST(imin,1,MPI_INTEGER,loc,MPI_COMM_WORLD,ierr)
      call MPI_BCAST(jmin,1,MPI_INTEGER,loc,MPI_COMM_WORLD,ierr)

      fmax=nmax(1)
      fmin=nmin(1)

    if(myid.eq.0)then
      write(6,100) amax,fmax,imax,jmax,1,    &
                   amin,fmin,imin,jmin,1
100   format(2x,a6,':',1x,g13.6,i5,i5,i5,    &
             4x,a6,':',1x,g13.6,i5,i5,i5)

      nstat = nstat + 1
      rstat(nstat) = fmax
      nstat = nstat + 1
      rstat(nstat) = fmin
    endif

      if(timestats.ge.1) time_stat=time_stat+mytime()

      end subroutine maxmin2d

  END MODULE maxmin_module