text:
      subroutine qqb_wbb_soft(P,msq)
c---Matrix element SUBTRACTION squared averaged over initial colors and spins
c     q(-p1)+qbar(-p2) -->  b(p5)+bb(p6) + W +g(p7)
c                                          |
c                                          --> nu(p3)+e^+(p4)
c
c   positively charged W only
c--all momenta incoming
      implicit none
      include 'constants.f'
      include 'qcdcouple.f'
      include 'masses.f'
      integer j,k

      double precision P(mxpart,4),msq(-nf:nf,-nf:nf),
     . msq0(-nf:nf,-nf:nf),s(mxpart,mxpart),facqqb,facqbq
      double precision eik17_2,eik17_5,eik17_6,eik27_5,eik27_6,eik57_6
      double precision eik27_1,eik57_1,eik67_1,eik57_2,eik67_2,eik67_5

      do j=-nf,nf
      do k=-nf,nf
      msq(j,k)=0d0
      enddo
      enddo
      call dotem(7,p,s)
c + p1Dp7^-1*p2Dp7^-1
c  * ( - n^-1*xlo^2*p1Dp2 )
c + p7Dp5^-1*p7Dp6^-1
c  * ( - n^-1*xlo^2*p5Dp6 )
c + p1Dp7^-1*p7Dp6^-1
c  * ( 2*n^-1*xlo^2*p1Dp6 )
c + p2Dp7^-1*p7Dp5^-1
c  * ( 2*n^-1*xlo^2*p2Dp5 )
c + p1Dp7^-1*p7Dp5^-1
c  * ( n*xlo^2*p1Dp5 - 2*n^-1*xlo^2*p1Dp5 )
c + p2Dp7^-1*p7Dp6^-1
c  * ( n*xlo^2*p2Dp6 - 2*n^-1*xlo^2*p2Dp6 )

c      eik12=s(1,2)/s(1,7)/s(2,7)
c      eik56=s(5,6)/s(5,7)/s(6,7)

c      eik15=s(1,5)/s(1,7)/s(5,7)
c      eik16=s(1,6)/s(1,7)/s(6,7)

c      eik25=s(2,5)/s(2,7)/s(5,7)
c      eik26=s(2,6)/s(2,7)/s(6,7)

c      facqqb=((+two*(eik17_6+eik67_1+eik27_5+eik57_2
c     . -eik27_6-eik67_2-eik17_5-eik57_1)
c     . -eik17_2-eik27_1-eik57_6-eik67_5
c     . )/xn
c     & +xn*(eik17_5+eik57_1+eik27_6+eik67_2))

      eik17_2=+s(1,2)/(s(1,7)+s(2,7))/s(1,7)*two*gsq
      eik27_1=+s(1,2)/(s(1,7)+s(2,7))/s(2,7)*two*gsq

      eik57_6=+s(5,6)/(s(5,7)+s(6,7))/s(5,7)*two*gsq
      eik67_5=+s(5,6)/(s(5,7)+s(6,7))/s(6,7)*two*gsq

      eik17_6=+s(1,6)/(-s(1,7)+s(6,7))/s(1,7)*two*gsq
      eik67_1=-s(1,6)/(-s(1,7)+s(6,7))/s(6,7)*two*gsq

      eik27_5=+s(2,5)/(-s(2,7)+s(5,7))/s(2,7)*two*gsq
      eik57_2=-s(2,5)/(-s(2,7)+s(5,7))/s(5,7)*two*gsq

      eik27_6=+s(2,6)/(-s(2,7)+s(6,7))/s(2,7)*two*gsq
      eik67_2=-s(2,6)/(-s(2,7)+s(6,7))/s(6,7)*two*gsq

      eik17_5=+s(1,5)/(-s(1,7)+s(5,7))/s(1,7)*two*gsq
      eik57_1=-s(1,5)/(-s(1,7)+s(5,7))/s(5,7)*two*gsq

c      facqbq=((+two*(eik27_6+eik67_2+eik17_5+eik57_1
c     . -eik17_6-eik67_1-eik27_5-eik57_2)
c     . -eik27_1-eik17_2-eik57_6-eik67_5
c     . )/xn
c     & +xn*(eik27_5+eik57_2+eik17_6+eik67_1))

      facqqb=(eik27_6+eik67_2+eik17_5+eik57_1)*(xn-two/xn)
     .      +(eik17_6+eik67_1+eik27_5+eik57_2)*(two/xn)
     .      +(eik27_1+eik17_2+eik57_6+eik67_5)*(-one/xn)
      facqbq=(eik17_6+eik67_1+eik27_5+eik57_2)*(xn-two/xn)
     .      +(eik27_6+eik67_2+eik17_5+eik57_1)*(two/xn)
     .      +(eik17_2+eik27_1+eik57_6+eik67_5)*(-one/xn)

c??????? exchange ???/

      call qqb_wbb(p,msq0)

c      write(6,*) 'p(1,4)',p(1,4)
c      write(6,*) 'p(2,4)',p(2,4)
c      write(6,*) 'msq0(2,-1)',msq0(2,-1),'msq0(-1,2)',msq0(-1,2)
c      write(6,*)

      do j=2,2
      do k=-1,-1
            if ((j .gt. 0) .and. (k.lt.0)) then
               msq(j,k)=facqqb*msq0(j,k)

       write(6,*) '-eik17_2*msq0(j,k)/xn',-eik17_2*msq0(j,k)/xn
       write(6,*) '-eik27_1*msq0(j,k)/xn',-eik27_1*msq0(j,k)/xn

       write(6,*) '-eik57_6*msq0(j,k)/xn',-eik57_6*msq0(j,k)/xn
       write(6,*) '-eik67_5*msq0(j,k)/xn',-eik67_5*msq0(j,k)/xn

       write(6,*) 'two/xn*(eik17_5+eik57_1)*msq0(j,k)',
     . +(xn-two/xn)*(eik17_5+eik57_1)*msq0(j,k)
       write(6,*) 'two/xn*(eik17_6+eik67_1)*msq0(j,k)',
     . two/xn*(eik17_6+eik67_1)*msq0(j,k)
       write(6,*) '(eik27_6+eik67_2)*msq0(j,k)',
     . +(xn-two/xn)*(eik27_6+eik67_2)*msq0(j,k)
       write(6,*) 'two/xn*(eik27_5+eik57_2)*msq0(j,k)',
     . +two/xn*(eik27_5+eik57_2)*msq0(j,k)

            elseif ((j .lt. 0) .and. (k.gt.0)) then
               msq(j,k)=facqbq*msq0(j,k)
            endif
      enddo
      enddo
      return
      end


module Globals
  implicit none

  ! Integer parameters
  integer, parameter :: gp_a = 48, gp_z = 20, gp_q = 7, gp_gamma = 3, agents = 10000, periods = 20001

  integer, dimension(:,:), allocatable :: shock_z, shock_q, shock_g, shock_mu

  ! Real parameters
  real(8), parameter :: cover_z = 2.d0 , cover_q = 2.d0, min_a = 0, max_a = 1440.0, &
                        tiny = 1.0d-10

  ! Real variables
  real(8) :: alpha, beta, gamma_bar, epsilon_gamma, rho_z, sigma_epsilon, sigma_q, lambda_e, &
  lambda_u, lambda_n, sigma, mu, b_0, b_bar, theta, delta, tau, int_rate, wage, &
  KLratio, average_z, T, new_KLratio, new_average_z, new_T, Erate, Urate, Nrate

  ! Real vectors
  real(8), dimension(gp_a) :: a_values
  real(8), dimension(gp_z) :: z_values, z_ssdist
  real(8), dimension(gp_q) :: q_values, q_trans
  real(8), dimension(gp_gamma) :: gamma_values, gamma_trans

  ! Real matrices
  real(8), dimension(:,:), allocatable :: shock_lm
  real(8), dimension(3,3) :: transitions
  real(8), dimension(gp_z,gp_z) :: z_trans
  real(8), dimension(0:1, 0:1) :: IB_trans
  real(8), dimension(gp_a,gp_z) :: N_vf, N_pf
  real(8), dimension(gp_a,gp_z,gp_q) :: W_vf, W_pf
  real(8), dimension(gp_a,gp_z,gp_gamma,0:1) :: U_vf, J_vf, U_pf
  real(8), dimension(gp_a,gp_z,gp_q,gp_gamma,0:1) :: V_vf

CONTAINS
  real(8) function u(consumption)
    implicit none
    real(8), intent(in) :: consumption

    u = log(max(1.0d-320,consumption))
  end function u

  ! Creates the value functions J and V given N, U, and W
  subroutine ValueFunctions(N,U,W,J,V)
    implicit none
    integer :: ind_a, ind_z, ind_g, ind_b, ind_q
    real(8), dimension(gp_a,gp_z), intent(in) :: N
    real(8), dimension(gp_a,gp_z,gp_q), intent(in) :: W
    real(8), dimension(gp_a,gp_z,gp_gamma,0:1), intent(in) :: U
    real(8), dimension(gp_a,gp_z,gp_gamma,0:1), intent(out) :: J
    real(8), dimension(gp_a,gp_z,gp_q,gp_gamma,0:1), intent(out) :: V

    do ind_a = 1, gp_a
    do ind_z = 1, gp_z
    do ind_g = 1, gp_gamma
    do ind_b = 0, 1
      J(ind_a,ind_z,ind_g,ind_b) = max(U(ind_a,ind_z,ind_g,ind_b), N(ind_a,ind_z))
      do ind_q = 1, gp_q
        V(ind_a,ind_z,ind_q,ind_g,ind_b) = max(W(ind_a,ind_z,ind_q),J(ind_a,ind_z,ind_g,ind_b))
      end do
    end do
    end do
    end do
    end do
  end subroutine ValueFunctions

  ! Computes equilibrium prices
  subroutine Prices(my_KL,my_r,my_w)
    implicit none
    real(8), intent(in) :: my_KL
    real(8), intent(out) ::  my_r, my_w

    my_w = (1.d0-theta)*(my_KL**theta)
    my_r = (theta*(my_KL**(theta-1.d0))) - delta
  end subroutine Prices

  ! Computes unemployment benefits
  real(8) function benefits(productivity)
    implicit none
    real(8), intent(in) :: productivity

    if ((productivity*b_0).lt.(average_z*b_bar)) then
      benefits = productivity*b_0*wage
    else
      benefits = average_z*b_bar*wage
    end if
  end function benefits
end module Globals


text:
integer function f()

end function

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
MODULE th

USE sdata, ONLY: DP

IMPLICIT NONE

SAVE

CONTAINS

SUBROUTINE th_iter(maxIter, oupt)
!
! Purpose:
!    To do thermal-hydrailics iteration
!

USE sdata, ONLY: nnod, ftem, mtem, cden, bcon, bpos, npow, pow, ppow,  &
                 zdel, node_nf, ix, iy, iz, th_err, node_nf, ix, iy, iz, &
                 th_niter, fer, ferc, ser, serc, get_time, th_time, Ke
USE cmfd, ONLY: outer_th
USE io, ONLY: ounit, scr
USE xsec, ONLY: XS_updt

IMPLICIT NONE

INTEGER, INTENT(IN) :: maxIter    ! maximum number of th iteration
INTEGER, INTENT(IN) :: oupt       ! print ouput option
REAL(DP), DIMENSION(nnod) :: pline
REAL(DP), DIMENSION(nnod) :: otem
INTEGER :: n, i
REAL(DP) :: st, fn

th_err = 1.
DO i = 1, maxIter
    ! Save old fuel temp
    otem = ftem

    CALL XS_updt(bcon, ftem, mtem, cden, bpos)

    ! Perform outer inner iteration
    CALL outer_th()

    !Get start th_time
    st = get_time()

    ! Calculate power density
    CALL PowDis(npow)

    ! Calculate linear power density for each nodes (W/cm)
    DO n = 1, nnod
        pline(n) = npow(n) * pow * ppow * 0.01_DP &
                 / (node_nf(ix(n),iy(n)) * zdel(iz(n)))     ! Linear power density (W/cm)
    END DO

    ! Update fuel, moderator temp. and coolant density
    CALL th_upd(pline)

    ! Get fuel absolute difference from current and previous th iteration
    CALL AbsE(ftem, otem, th_err)

    !Get th_time
    fn = get_time()
    th_time = th_time + (fn-st)

    if (oupt > 0) then
      write(ounit,'(I5,F13.6,3ES15.5)') i, Ke, ser, fer, th_err        ! Write outer iteration evolution
      if (scr) write(*,'(I5,F13.6,3ES15.5)') i, Ke, ser, fer, th_err   ! Write outer iteration evolution
    end if

    ! If error is small enough
    IF (th_err < 0.01 .and. fer < ferc .and. ser < serc) EXIT

END DO

IF (i-1 == th_niter) THEN
  WRITE(ounit,*) '  MAXIMUM TH ITERATION REACHED.'
  WRITE(ounit,*) '  CALCULATION MIGHT BE NOT CONVERGED OR CHANGE ITERATION CONTROL'
  WRITE(*,*) '  MAXIMUM TH ITERATION REACHED.'
  WRITE(*,*) '  CALCULATION MIGHT BE NOT CONVERGED OR CHANGE ITERATION CONTROL'
  STOP
END IF

END SUBROUTINE th_iter

!****************************************************************************!

subroutine PowDis (p)

!
! Purpose:
!    To calculate power for each nodes
!

USE sdata, ONLY: ng, nnod, sigf, f0, vdel, powtot, mode
USE io,    ONLY: ounit

implicit none

real(dp), dimension(:), intent(out) :: p
integer :: g, n
real(dp) :: pow

p = 0._dp
do g= 1, ng
    do n= 1, nnod
      pow = f0(n,g) * sigf(n,g) * vdel(n)
      if (pow < 0.) pow = 0.
      p(n) = p(n) + pow
    end do
end do

! Normalize to 1._DP
powtot = 0._DP
do n = 1, nnod
    powtot = powtot + p(n)
end do

if (powtot <= 0 .AND. mode /= 'FIXEDSRC') THEN
   write(ounit, *) '   ERROR: TOTAL NODES POWER IS ZERO OR LESS'
   write(ounit, *) '   STOP IN subroutine POWDIS'
   STOP
end if

if (powtot > 0.0) then
   do n = 1, nnod
       p(n) = p(n) / powtot
   end do
end if

end subroutine PowDis

!****************************************************************************!

SUBROUTINE AbsE(newF, oldF, rel)

  !
  ! Purpose:
  !    To calculate Max Relative error

USE sdata, ONLY: nnod

IMPLICIT NONE

REAL(DP), DIMENSION(:), INTENT(IN) :: newF, oldF
REAL(DP), INTENT(OUT) :: rel

REAL(DP) :: error
INTEGER :: n

rel = 0.

DO n= 1, nnod
    IF (ABS(newF(n)) > 1.e-10_DP) THEN
        error = ABS(newF(n) - oldF(n))
        IF (error > rel) rel = error
    END IF
END DO

END SUBROUTINE AbsE

!****************************************************************************!

SUBROUTINE par_ave(par, ave)
!
! Purpose:
!    To calculate average fuel temp (only for active core)
!

USE sdata, ONLY: vdel, nnod, ng, nuf

IMPLICIT NONE

REAL(DP), DIMENSION(:), INTENT(IN) :: par
REAL(DP), INTENT(OUT) :: ave
REAL(DP) :: dum, dum2
INTEGER :: n

dum = 0.; dum2 = 0.
DO n = 1, nnod
   IF (nuf(n,ng) > 0.) THEN
      dum = dum + par(n) * vdel(n)
      dum2 = dum2 + vdel(n)
   END IF
END DO

ave = dum / dum2

END SUBROUTINE par_ave

!****************************************************************************!

SUBROUTINE par_ave_out(par, ave)
!
! Purpose:
!    To calculate average fuel temp (only for active core)
!

USE sdata, ONLY: vdel, nnod, iz, nzz, nuf, ng, ix, iy, xyz, nxx, nyy, nzz, ystag

IMPLICIT NONE

REAL(DP), DIMENSION(:), INTENT(IN) :: par
REAL(DP), INTENT(OUT) :: ave
REAL(DP) :: dum, dum2
INTEGER, DIMENSION(nxx,nyy) :: zmax
INTEGER :: n, i, j, k

! get number of nodex in axial direction from bottom -> fuel
DO j = 1, nyy
  DO i = ystag(j)%smin, ystag(j)%smax
    zmax(i,j) = 0
    DO k = 1, nzz/2
      if (nuf(xyz(i,j,k),ng) < 1.e-5_dp) zmax(i,j) = zmax(i,j) + 1
    END DO
  END DO
END DO
! get number of nodex in axial direction from fuel -> top reflectors
DO j = 1, nyy
  DO i = ystag(j)%smin, ystag(j)%smax
    DO k = 1, nzz
      if (nuf(xyz(i,j,k),ng) > 1.e-5) zmax(i,j) = zmax(i,j) + 1
    END DO
  END DO
END DO

dum = 0.; dum2 = 0.
DO n = 1, nnod
   IF (iz(n) == zmax(ix(n),iy(n)) .AND. nuf(n,ng) > 1.e-5) THEN
      dum = dum + par(n) * vdel(n)
      dum2 = dum2 + vdel(n)
   END IF
END DO

ave = dum / dum2

END SUBROUTINE par_ave_out

!****************************************************************************!

SUBROUTINE par_max(par, pmax)
!
! Purpose:
!    To calculate maximum fuel tem, coolant tem, and density
!

USE sdata, ONLY: nnod

IMPLICIT NONE

REAL(DP), DIMENSION(:), INTENT(IN) :: par
REAL(DP), INTENT(OUT) :: pmax
INTEGER :: n

pmax = 0.
DO n = 1, nnod
   IF (par(n) > pmax) pmax = par(n)
END DO

END SUBROUTINE par_max

!****************************************************************************!

SUBROUTINE getent(t,ent)
!
! Purpose:
!    To get enthalpy for given coolant temp. from steam table
!

USE sdata, ONLY: stab, ntem
USE io, ONLY : ounit

IMPLICIT NONE

REAL(DP), INTENT(IN) :: t
REAL(DP), INTENT(OUT) :: ent
REAL(DP) :: t1, ent1
REAL(DP) :: t2, ent2
INTEGER :: i

IF ((t < stab(1,1)) .OR. (t > stab(ntem,1))) THEN
    WRITE(ounit,*) '  Coolant temp. : ', t, 'K'
    WRITE(ounit,*) '  ERROR : MODERATOR TEMP. IS OUT OF THE RANGE OF DATA IN THE STEAM TABLE'
    WRITE(ounit,*) '  CHECK INPUT COOLANT MASS FLOW RATE OR CORE POWER'
    WRITE(*,*) '  Coolant temp. : ', t, 'K'
    WRITE(*,*) '  ERROR : MODERATOR TEMP. IS OUT OF THE RANGE OF DATA IN THE STEAM TABLE'
    WRITE(*,*) '  CHECK INPUT COOLANT MASS FLOW RATE OR CORE POWER'
    STOP
END IF

t2 = stab(1,1); ent2 = stab(1,3)
DO i = 2, ntem
    t1 = t2
    ent1 = ent2
    t2 = stab(i,1); ent2 = stab(i,3)
    IF ((t >= t1) .AND. (t <= t2)) THEN
        ent = ent1 + (t - t1) / (t2 - t1) * (ent2 - ent1)
        EXIT
    END IF
END DO

END SUBROUTINE getent

!****************************************************************************!

SUBROUTINE gettd(ent,t,rho,prx,kvx,tcx,Rx)
!
! Purpose:
!    To get enthalpy for given coolant temp. from steam table
!

USE sdata, ONLY: stab, ntem
USE io, ONLY : ounit

IMPLICIT NONE

REAL(DP), INTENT(IN) :: ent
REAL(DP), INTENT(OUT) :: t, rho, prx, kvx, tcx
REAL(DP), INTENT(OUT), OPTIONAL :: Rx
REAL(DP) :: ratx

INTEGER :: i, i1, i2

! Get two closest interpolation points
IF (ent >= stab(1,3) .AND. ent <= stab(ntem,3)) THEN  !If enthalpy inside data range
  DO i = 2, ntem
    IF (ent >= stab(i-1,3) .AND. ent <= stab(i,3)) THEN
      i1 = i-1
      i2 = i
      EXIT
    END IF
  END DO
ELSE IF (ent < stab(1,3)  .AND. (stab(1,3) - ent) / stab(1,3) < 0.1) THEN !If 10% lower than min. steam table data
  i1 = 1
  i2 = 2
ELSE IF (ent > stab(ntem,3) .AND. (ent - stab(ntem,3)) / stab(ntem,3) < 0.1) THEN !if 10% higher than max. steam table data
  i1 = ntem-1
  i2 = ntem
ELSE
  WRITE(ounit,1557) ent/1000.
  WRITE(*,1557) ent/1000.
  WRITE(ounit,*) '   CHECK INPUT COOLANT MASS FLOW RATE OR CORE POWER'
  WRITE(*,*) '   CHECK INPUT COOLANT MASS FLOW RATE OR CORE POWER'
  STOP
  1557 FORMAT(2X, '  ERROR: ENTHALPY', F8.1 ,' KJ/Kg &
  & IS OUT OF THE RANGE IN THE STEAM TABLE')
END IF

! Interpolate
ratx = (ent - stab(i1,3)) / (stab(i2,3) - stab(i1,3))
t   = stab(i1,1) + ratx * (stab(i2,1) - stab(i1,1))
rho = stab(i1,2) + ratx * (stab(i2,2) - stab(i1,2))
prx = stab(i1,4) + ratx * (stab(i2,4) - stab(i1,4))
kvx = stab(i1,5) + ratx * (stab(i2,5) - stab(i1,5))
tcx = stab(i1,6) + ratx * (stab(i2,6) - stab(i1,6))
IF (PRESENT(Rx)) THEN
  Rx = 1000._DP * (stab(i2,2) - stab(i1,2)) / (stab(i2,3) - stab(i1,3))
END IF

END SUBROUTINE gettd

!****************************************************************************!

REAL(DP) FUNCTION getkc(t)
!
! Purpose:
!    To calculate thermal conductivity of cladding
!

IMPLICIT NONE

REAL(DP), INTENT(IN) :: t

getkc = 7.51_DP + 2.09e-2_DP*t - 1.45e-5_DP*t**2 + 7.67e-9_DP*t**3

END FUNCTION getkc

!****************************************************************************!

REAL(DP) FUNCTION getkf(t)
!
! Purpose:
!    To calculate thermal conductivity of fuel
!

IMPLICIT NONE

REAL(DP), INTENT(IN) :: t

getkf = 1.05_DP + 2150.0_DP / (t - 73.15_DP)

END FUNCTION getkf

!****************************************************************************!

REAL(DP) FUNCTION getcpc(t)
!
! Purpose:
!    To calculate specific heat capacity of cladding
!

IMPLICIT NONE

REAL(DP), INTENT(IN) :: t

getcpc = 252.54_DP + 0.11474_DP*t

END FUNCTION getcpc

!****************************************************************************!

REAL(DP) FUNCTION getcpf(t)
!
! Purpose:
!    To calculate specific heat capacity of fuel
!

IMPLICIT NONE

REAL(DP), INTENT(IN) :: t

getcpf = 162.3_DP + 0.3038_DP*t - 2.391e-4_DP*t**2 + 6.404e-8_DP*t**3

END FUNCTION getcpf

!****************************************************************************!

SUBROUTINE TridiaSolve(a,b,c,d,x)
!
! Purpose:
!    To solve tridiagonal matrix
!

IMPLICIT NONE

REAL(DP), DIMENSION(:), INTENT(INOUT) :: a, b, c, d
REAL(DP), DIMENSION(:), INTENT(OUT) :: x

INTEGER :: i, n

n = SIZE(d)

! Gauss Elimination
c(1) = c(1)/b(1)
d(1) = d(1)/b(1)
DO i = 2, n
    c(i) = c(i) / (b(i) - a(i) * c(i-1))
    d(i) = (d(i) - a(i) * d(i-1)) / (b(i) - a(i) * c(i-1))
END DO

! Back Substitution
x(n) = d(n)
DO i = n-1, 1, -1
    x(i) = d(i) - c(i) * x(i+1)
END DO

END SUBROUTINE TridiaSolve

!****************************************************************************!

REAL(DP) FUNCTION geths(xden, tc, kv, Pr)
!
! Purpose:
!    To calculate heat transfer coef.
!

USE sdata, ONLY: dh, farea, cflow

IMPLICIT NONE

REAL(DP), INTENT(IN) :: xden  ! coolant densisty
REAL(DP), INTENT(IN) :: tc  ! coolant thermal conductivity
REAL(DP), INTENT(IN) :: kv  ! kinematic viscosity
REAL(DP), INTENT(IN) :: Pr  ! Prandtl Number

REAL(DP) :: cvelo, Nu, Re

cvelo = cflow / (farea * xden * 1000._DP)        ! Calculate flow velocity (m/s)
Re = cvelo * dh / (kv * 1.e-6_DP)                 ! Calculate Reynolds Number
Nu = 0.023_DP*(Pr**0.4_DP)*(Re**0.8_DP)                ! Calculate Nusselt Number
geths = (tc / dh) * Nu                        ! Calculate heat transfer coefficient

END FUNCTION geths

!****************************************************************************!

SUBROUTINE th_trans(xpline, h)

!
! Purpose:
!    To perform fuel pin thermal transient
!

USE sdata, ONLY: nnod, mtem, cden, ftem, tin, cflow, nxx, nyy, cf, ent, heatf, &
                 tfm, nt, rpos, rdel, rf, rg, rc, farea, dia, pi, zdel, &
                 ix, iy, iz, frate, th_time, get_time

IMPLICIT NONE

REAL(DP), DIMENSION(:), INTENT(IN) :: xpline    ! Linear Power Density (W/cm)
REAL(DP), INTENT(IN) :: h                       ! Time step

INTEGER :: i, j, k, n
REAL(DP), DIMENSION(nt+1) :: a, b, c, d
REAL(DP) :: hs, hg = 1.d4, kt , kt1, kt2          ! coolant heat transfer coef., gap heat transfer coef, and thermal conductivity
REAL(DP) :: alpha = 0.7_DP
REAL(DP) :: xa, xc
REAL(DP) :: fdens = 10.412e3            ! UO2 density (kg/m3)
REAL(DP) :: cdens = 6.6e3               ! Cladding density (kg/m3)
REAL(DP) :: cp                          ! Specific heat capacity
REAL(DP) :: eps, eta
REAL(DP) :: mdens, vol                  ! Coolant density and channel volume
REAL(DP), DIMENSION(nnod) :: entp        ! previous enthalpy

REAL(DP) :: pdens      ! power densisty  (W/m3)
REAL(DP) :: enti       ! Coolant inlet enthalpy
REAL(DP), DIMENSION(nxx, nyy) :: entm, bfrate
REAL(DP) :: cpline     ! Coolant Linear power densisty (W/m)
REAL(DP) :: Pr, kv, tcon ! Coolant Prandtl Number, Kinematic viscosity, and thermal conductivity
real(dp) :: R
logical  :: first = .true.
real(dp) :: st, fn

!Get start th_time
st = get_time()

!set initial tridiagonal matrix element a, b
a = 0._dp; b = 0._dp; c = 0._dp;

if (first) then
  allocate(frate(nnod))
  frate = cflow
  first = .false.
end if
! frate = cflow
CALL getent(tin, enti)
entp = ent

DO n = 1, nnod
   mdens = cden(n) * 1000._DP                                    ! Coolant density (kg/m3)
   cpline = heatf(n) * pi * dia + cf * xpline(n) * 100._DP       ! Coolant Linear power densisty (W/m)
   vol   = farea * zdel(iz(n)) * 0.01_DP                             ! channel node volume
   i = ix(n); j = iy(n); k = iz(n)

   IF (k == 1) THEN                                              ! Calculate coolant enthalpy
       eps = mdens * vol / h
       ent(n) = (cpline * zdel(iz(n)) * 0.01_DP + 2._DP * frate(n) * enti &
              + eps * entp(n)) / (eps + 2._DP * frate(n))                             ! Calculate enthalpy
       CALL gettd(ent(n), mtem(n), cden(n), Pr, kv, tcon, R)        ! Get corresponding temp and density
       entm(i,j)   = 2._DP * ent(n) - enti
       frate(n)    = cflow - 0.5_dp * vol / h * R * (ent(n) - entp(n))
       bfrate(i,j) = 2._DP * frate(n) - cflow
   ELSE
       eps = mdens * vol / h
       ent(n) = (cpline * zdel(iz(n)) * 0.01_DP + 2._DP * frate(n) &
              * entm(ix(n),iy(n)) + eps * entp(n)) / (eps + 2._DP * frate(n))
       CALL gettd(ent(n), mtem(n), cden(n), Pr, kv, tcon, R)
       entm(i,j)   = 2._DP * ent(n) - entm(i,j)
       frate(n)    = bfrate(i,j) - 0.5_dp * vol / h * R * (ent(n) - entp(n))
       bfrate(i,j) = 2._DP * frate(n) - bfrate(i,j)
   END IF

   hs = geths(cden(n), Pr, kv, tcon)                                               ! Calculate heat transfer coef
   pdens = 100._DP * xpline(n) / (pi * rf**2)                ! Fuel pin Power Density (W/m3)

   ! Calculate tridiagonal matrix: a, b, c and source: d
   ! For nt=1 [FUEL CENTERLINE]
   kt1 = getkf(tfm(n,1))                                                     ! Get thermal conductivity
   kt2 = getkf(tfm(n,2))
   kt  = 2._DP * kt1 * kt2 / (kt1 + kt2)
   cp = getcpf(tfm(n,1))                                                           ! Get specific heat capacity
   eta = fdens * cp * rpos(1)**2 / (2._DP * h)
   xc  = kt * rpos(1) / rdel(1)
   b(1) =  xc + eta
   c(1) = -xc
   d(1) = pdens * 0.5_DP * rpos(1)**2 + eta * tfm(n,1)

   DO i = 2, nt-2
       kt1 = kt2
       kt2 = getkf(tfm(n,i+1))
       kt  = 2._DP * kt1 * kt2 / (kt1 + kt2)
       cp = getcpf(tfm(n,i))
       eta = fdens * cp * (rpos(i)**2 - rpos(i-1)**2) / (2. * h)
       xa = xc
       xc = kt * rpos(i) / rdel(i)
       a(i) = -xa
       b(i) =  xa + xc + eta
       c(i) = -xc
       d(i) = pdens * 0.5_DP * (rpos(i)**2 - rpos(i-1)**2) + eta * tfm(n,i)
   END DO

   ! For nt-1 [FUEL-GAP INTERFACE]
   cp = getcpf(tfm(n,nt-1))
   eta = fdens * cp * (rf**2 - rpos(nt-2)**2) / (2. * h)
   xa = xc
   xc = rg * hg
   a(nt-1) = -xa
   b(nt-1) =  xa + xc + eta
   c(nt-1) = -xc
   d(nt-1) = pdens * 0.5_DP * (rf**2 - rpos(nt-2)**2) + eta * tfm(n,nt-1)

   ! For nt [GAP-CLADDING INTERFACE]
   kt1 = getkc(tfm(n,nt))
   kt2 = getkc(tfm(n,nt+1))
   kt  = 2._DP * kt1 * kt2 / (kt1 + kt2)     ! For cladding
   cp = getcpc(tfm(n,nt))
   eta = cdens * cp * (rpos(nt)**2 - rg**2) / (2. * h)
   xa = xc
   xc = kt * rpos(nt) / rdel(nt)
   a(nt) = -xa
   b(nt) =  xa + xc + eta
   c(nt) = -xc
   d(nt) = eta * tfm(n,nt)

   ! For nt+1  [CLADDING-COOLANT INTERFACE]
   cp = getcpc(tfm(n,nt+1))
   eta = cdens * cp * (rc**2 - rpos(nt)**2) / (2. * h)
   xa = xc
   xc = rc * hs
   a(nt+1) = -xa
   b(nt+1) =  xa + xc + eta
   d(nt+1) = rc * hs * mtem(n) + eta * tfm(n,nt+1)

   ! Solve tridiagonal matrix
   CALL TridiaSolve(a, b, c, d, tfm(n, :))

   ! Get lumped fuel temp
   ftem(n) = (1.-alpha) * tfm(n, 1) + alpha * tfm(n, nt-1)

   ! Calculate heat flux
   heatf(n) = hs * (tfm(n, nt+1) - mtem(n))
END DO

!Get th_time
fn = get_time()
th_time = th_time + (fn-st)

END SUBROUTINE th_trans

!****************************************************************************!

SUBROUTINE th_upd(xpline)

!
! Purpose:
!    To update thermal parameters
!

USE sdata, ONLY: nnod, mtem, cden, ftem, tin, ix, iy, iz, nxx, nyy, cflow, cf,  &
                 ent, heatf, tfm, nt, rpos, rdel, rf, rg, rc, pi, zdel, dia

IMPLICIT NONE

REAL(DP), DIMENSION(:), INTENT(IN) :: xpline    ! Linear Power Density (W/cm)

INTEGER :: i, n
REAL(DP), DIMENSION(nt+1) :: a, b, c, d
REAL(DP) :: hs, Hg = 1.D4, kt, kt1, kt2
REAL(DP) :: alp = 0.7_DP
REAL(DP) :: xa, xc
REAL(DP) :: pdens      ! power densisty  (W/m3)
REAL(DP) :: enti       ! Coolant inlet enthalpy
REAL(DP), DIMENSION(nxx, nyy) :: entm   ! enthalpy at node boundary
REAL(DP) :: cpline     ! Coolant Linear power densisty (W/m)
REAL(DP) :: Pr, kv, tcon ! Coolant Prandtl Number, Kinematic viscosity, and thermal conductivity
REAL(DP) :: zd  ! zdel in meter

!set initial tridiagonal matrix element a, b
a = 0._dp; b = 0._dp; c = 0._dp;

CALL getent(tin, enti)

DO n = 1, nnod
  cpline = heatf(n) * pi * dia  + cf * xpline(n) * 100._DP          ! Coolant Linear power densisty (W/m)
  zd = zdel(iz(n)) * 0.01_DP
  IF (iz(n) == 1) THEN                                              ! For most bootom channel
    ent(n) = enti + 0.5_DP * cpline * zd / cflow    ! Calculate coolant enthalpy
    CALL gettd(ent(n), mtem(n), cden(n), Pr, kv, tcon)             ! Get corresponding temp and density
    entm(ix(n),iy(n)) = 2._DP * ent(n) - enti                      ! Extrapolate enthalpy at node boundary
  ELSE
    ent(n) = entm(ix(n),iy(n)) + 0.5_DP * cpline * zd / cflow
    CALL gettd(ent(n), mtem(n), cden(n), Pr, kv, tcon)
    entm(ix(n),iy(n)) = 2._DP * ent(n) - entm(ix(n),iy(n))
  END IF

  hs = geths(cden(n), Pr, kv, tcon)
  pdens = (1._DP - cf) * 100._DP * xpline(n) / (pi * rf**2)        ! Fuel pin Power Density (W/m3)

  ! Calculate tridiagonal matrix: a, b, c and source: d
  ! For nt=1 [FUEL CENTERLINE]
  kt1 = getkf(tfm(n,1))                                         ! Get thermal conductivity
  kt2 = getkf(tfm(n,2))
  kt  = 2._DP * kt1 * kt2 / (kt1 + kt2)
  xc  = kt * rpos(1) / rdel(1)
  b(1) =  xc
  c(1) = -xc
  d(1) = pdens * 0.5_DP * rpos(1)**2

  DO i = 2, nt-2
    kt1 = kt2
    kt2 = getkf(tfm(n,i+1))
    kt  = 2._DP * kt1 * kt2 / (kt1 + kt2)
    xa = xc
    xc = kt * rpos(i) / rdel(i)
    a(i) = -xa
    b(i) =  xa + xc
    c(i) = -xc
    d(i) = pdens * 0.5_DP * (rpos(i)**2 - rpos(i-1)**2)
  END DO

  ! For nt-1 [FUEL-GAP INTERFACE]
  xa = xc
  xc = rg * Hg
  a(nt-1) = -xa
  b(nt-1) =  xa + xc
  c(nt-1) = -xc
  d(nt-1) = pdens * 0.5_DP * (rf**2 - rpos(nt-2)**2)

  ! For nt [GAP-CLADDING INTERFACE]
  kt1 = getkc(tfm(n,nt))
  kt2 = getkc(tfm(n,nt+1))
  kt  = 2._DP * kt1 * kt2 / (kt1 + kt2)     ! For cladding
  xa = xc
  xc = kt * rpos(nt) / rdel(nt)
  a(nt) = -xa
  b(nt) =  xa + xc
  c(nt) = -xc
  d(nt) = 0.

  ! For nt+1  [CLADDING-COOLANT INTERFACE]
  xa = xc
  a(nt+1) = -xa
  b(nt+1) =  xa + hs * rc
  d(nt+1) = rc * hs * mtem(n)

  ! Solve tridiagonal matrix
  CALL TridiaSolve(a, b, c, d, tfm(n, :))

  ! Get lumped fuel temp
  ftem(n) = (1.-alp) * tfm(n, 1) + alp * tfm(n, nt-1)

  ! Calculate heat flux
  heatf(n) = hs * (tfm(n, nt+1) - mtem(n))
END DO

END SUBROUTINE th_upd

END MODULE th

max_stars_repo_name: khurrumsaleem/KOMODO

max_issues_repo_name: khurrumsaleem/KOMODO

max_forks_repo_name: khurrumsaleem/KOMODO


text:
C    Copyright(C) 1988-2017 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    Redistribution and use in source and binary forms, with or without
C    modification, are permitted provided that the following conditions are
C    met:
C
C    * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C    * Redistributions in binary form must reproduce the above
C      copyright notice, this list of conditions and the following
C      disclaimer in the documentation and/or other materials provided
C      with the distribution.
C
C    * Neither the name of NTESS nor the names of its
C      contributors may be used to endorse or promote products derived
C      from this software without specific prior written permission.
C
C    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C $Id: range.f,v 1.2 1999/02/16 21:38:01 gdsjaar Exp $
C $Log: range.f,v $
C Revision 1.2  1999/02/16 21:38:01  gdsjaar
C Converted to read exodusII database format.  Somewhat tested, not
C ready for production yet.
C
C Revision 1.1.1.1  1991/02/21 15:45:22  gdsjaar
C NUMBERS: Greg Sjaardema, initial Unix release
C
c Revision 1.1  1991/02/21  15:45:21  gdsjaar
c Initial revision
c
      SUBROUTINE RANGE (LEN, LIST, IOMIN, IOMAX)
      LOGICAL LIST(*), INRNG
      INTEGER IRANGE(3)
      CHARACTER*80 LINE

      INRNG = .FALSE.
      IRINC = 0

      DO 10 I=1, LEN
         IF (LIST(I)) THEN
            IBEG = I+1
            LASTSL = I
            GO TO 20
         END IF
   10 CONTINUE

   20 CONTINUE
      DO 40 I = IBEG, LEN
         IF (LIST(I)) THEN
            IRINC1 = I - LASTSL
            IF (IRINC .EQ. 0) IRINC = IRINC1
            IF (.NOT. INRNG) THEN
               IRBEG = LASTSL
               IRINC = IRINC1
               INRNG = .TRUE.
            ELSE IF (INRNG .AND. IRINC1 .NE. IRINC) THEN
               IREND = I
               INRNG = .FALSE.
               IRANGE(1) = IRBEG
               IRANGE(2) = IREND - IRINC1
               IRANGE(3) = IRINC
               LINE = ' '
               CALL FFADDV (IRANGE, LINE)
               DO 30 IO = IOMIN, IOMAX
                  CALL LOGERR ('CMDSPEC', LINE(:LENSTR(LINE)), IO)
   30          CONTINUE
               IRINC = IRINC1
            END IF
            LASTSL = I
         END IF
   40 CONTINUE
      IRANGE(1) = IRBEG
      IRANGE(2) = LASTSL
      IRANGE(3) = IRINC
      LINE = ' '
      CALL FFADDV (IRANGE, LINE)
      DO 50 IO = IOMIN, IOMAX
         CALL LOGERR ('CMDSPEC', LINE(:LENSTR(LINE)), IO)
   50 CONTINUE
      RETURN
      END

ext: f

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


text:
! RUN: %python %S/test_errors.py %s %flang_fc1

! Test that the interface of specific intrinsics passed as dummy arguments
! are correctly validated against actual arguments explicit interface.

  intrinsic :: abs, dabs
  interface
    subroutine foo(f)
      interface
        function f(x)
          real :: f
          real, intent(in) :: x
        end function
      end interface
    end subroutine

    subroutine foo2(f)
      interface
        function f(x)
          double precision :: f
          double precision, intent(in) :: x
        end function
      end interface
    end subroutine
  end interface

  ! OK
  call foo(abs)

  ! OK
  call foo2(dabs)

  !ERROR: Actual procedure argument has interface incompatible with dummy argument 'f='
  call foo(dabs)

  !ERROR: Actual procedure argument has interface incompatible with dummy argument 'f='
  call foo2(abs)
end

max_stars_repo_name: mkinsner/llvm

max_issues_repo_name: mkinsner/llvm

max_forks_repo_name: mkinsner/llvm


text:
c///////////////////////////////////////////////////////////////////////////
c---------------------------------------------------------------------------
      SUBROUTINE SGECO(A,LDA,N,IPVT,RCOND,Z)
c---------------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER LDA,N,IPVT(1)
      REAL A(LDA,1),Z(1)
      REAL RCOND
C
C     SGECO FACTORS A REAL MATRIX BY GAUSSIAN ELIMINATION
C     AND ESTIMATES THE CONDITION OF THE MATRIX.
C
C     IF  RCOND  IS NOT NEEDED, SGEFA IS SLIGHTLY FASTER.
C     TO SOLVE  A*X = B , FOLLOW SGECO BY SGESL.
C     TO COMPUTE  INVERSE(A)*C , FOLLOW SGECO BY SGESL.
C     TO COMPUTE  DETERMINANT(A) , FOLLOW SGECO BY SGEDI.
C     TO COMPUTE  INVERSE(A) , FOLLOW SGECO BY SGEDI.
C
C     ON ENTRY
C
C        A       REAL(LDA, N)
C                THE MATRIX TO BE FACTORED.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY  A .
C
C        N       INTEGER
C                THE ORDER OF THE MATRIX  A .
C
C     ON RETURN
C
C        A       AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
C                WHICH WERE USED TO OBTAIN IT.
C                THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE
C                L  IS A PRODUCT OF PERMUTATION AND UNIT LOWER
C                TRIANGULAR MATRICES AND  U  IS UPPER TRIANGULAR.
C
C        IPVT    INTEGER(N)
C                AN INTEGER VECTOR OF PIVOT INDICES.
C
C        RCOND   REAL
C                AN ESTIMATE OF THE RECIPROCAL CONDITION OF  A .
C                FOR THE SYSTEM  A*X = B , RELATIVE PERTURBATIONS
C                IN  A  AND  B  OF SIZE  EPSILON  MAY CAUSE
C                RELATIVE PERTURBATIONS IN  X  OF SIZE  EPSILON/RCOND .
C                IF  RCOND  IS SO SMALL THAT THE LOGICAL EXPRESSION
C                           1.0 + RCOND .EQ. 1.0
C                IS TRUE, THEN  A  MAY BE SINGULAR TO WORKING
C                PRECISION.  IN PARTICULAR,  RCOND  IS ZERO  IF
C                EXACT SINGULARITY IS DETECTED OR THE ESTIMATE
C                UNDERFLOWS.
C
C        Z       REAL(N)
C                A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT.
C                IF  A  IS CLOSE TO A SINGULAR MATRIX, THEN  Z  IS
C                AN APPROXIMATE NULL VECTOR IN THE SENSE THAT
C                NORM(A*Z) = RCOND*NORM(A)*NORM(Z) .
C
C     LINPACK. THIS VERSION DATED 08/14/78 .
C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C     SUBROUTINES AND FUNCTIONS
C
C     LINPACK SGEFA
C     BLAS SAXPY,SDOT,SSCAL,SASUM
C     FORTRAN ABS,AMAX1,SIGN
C
C     INTERNAL VARIABLES
C
      REAL SDOT,EK,T,WK,WKM
      REAL ANORM,S,SASUM,SM,YNORM
      INTEGER INFO,J,K,KB,KP1,L
C
C
C     COMPUTE 1-NORM OF A
C
      ANORM = 0.0E0
      DO 10 J = 1, N
         ANORM = AMAX1(ANORM,SASUM(N,A(1,J),1))
   10 CONTINUE
C
C     FACTOR
C
      CALL SGEFA(A,LDA,N,IPVT,INFO)
C
C     RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) .
C     ESTIMATE = NORM(Z)/NORM(Y) WHERE  A*Z = Y  AND  TRANS(A)*Y = E .
C     TRANS(A)  IS THE TRANSPOSE OF A .  THE COMPONENTS OF  E  ARE
C     CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W  WHERE
C     TRANS(U)*W = E .  THE VECTORS ARE FREQUENTLY RESCALED TO AVOID
C     OVERFLOW.
C
C     SOLVE TRANS(U)*W = E
C
      EK = 1.0E0
      DO 20 J = 1, N
         Z(J) = 0.0E0
   20 CONTINUE
      DO 100 K = 1, N
         IF (Z(K) .NE. 0.0E0) EK = SIGN(EK,-Z(K))
         IF (ABS(EK-Z(K)) .LE. ABS(A(K,K))) GO TO 30
            S = ABS(A(K,K))/ABS(EK-Z(K))
            CALL SSCAL(N,S,Z,1)
            EK = S*EK
   30    CONTINUE
         WK = EK - Z(K)
         WKM = -EK - Z(K)
         S = ABS(WK)
         SM = ABS(WKM)
         IF (A(K,K) .EQ. 0.0E0) GO TO 40
            WK = WK/A(K,K)
            WKM = WKM/A(K,K)
         GO TO 50
   40    CONTINUE
            WK = 1.0E0
            WKM = 1.0E0
   50    CONTINUE
         KP1 = K + 1
         IF (KP1 .GT. N) GO TO 90
            DO 60 J = KP1, N
               SM = SM + ABS(Z(J)+WKM*A(K,J))
               Z(J) = Z(J) + WK*A(K,J)
               S = S + ABS(Z(J))
   60       CONTINUE
            IF (S .GE. SM) GO TO 80
               T = WKM - WK
               WK = WKM
               DO 70 J = KP1, N
                  Z(J) = Z(J) + T*A(K,J)
   70          CONTINUE
   80       CONTINUE
   90    CONTINUE
         Z(K) = WK
  100 CONTINUE
      S = 1.0E0/SASUM(N,Z,1)
      CALL SSCAL(N,S,Z,1)
C
C     SOLVE TRANS(L)*Y = W
C
      DO 120 KB = 1, N
         K = N + 1 - KB
         IF (K .LT. N) Z(K) = Z(K) + SDOT(N-K,A(K+1,K),1,Z(K+1),1)
         IF (ABS(Z(K)) .LE. 1.0E0) GO TO 110
            S = 1.0E0/ABS(Z(K))
            CALL SSCAL(N,S,Z,1)
  110    CONTINUE
         L = IPVT(K)
         T = Z(L)
         Z(L) = Z(K)
         Z(K) = T
  120 CONTINUE
      S = 1.0E0/SASUM(N,Z,1)
      CALL SSCAL(N,S,Z,1)
C
      YNORM = 1.0E0
C
C     SOLVE L*V = Y
C
      DO 140 K = 1, N
         L = IPVT(K)
         T = Z(L)
         Z(L) = Z(K)
         Z(K) = T
         IF (K .LT. N) CALL SAXPY(N-K,T,A(K+1,K),1,Z(K+1),1)
         IF (ABS(Z(K)) .LE. 1.0E0) GO TO 130
            S = 1.0E0/ABS(Z(K))
            CALL SSCAL(N,S,Z,1)
            YNORM = S*YNORM
  130    CONTINUE
  140 CONTINUE
      S = 1.0E0/SASUM(N,Z,1)
      CALL SSCAL(N,S,Z,1)
      YNORM = S*YNORM
C
C     SOLVE  U*Z = V
C
      DO 160 KB = 1, N
         K = N + 1 - KB
         IF (ABS(Z(K)) .LE. ABS(A(K,K))) GO TO 150
            S = ABS(A(K,K))/ABS(Z(K))
            CALL SSCAL(N,S,Z,1)
            YNORM = S*YNORM
  150    CONTINUE
         IF (A(K,K) .NE. 0.0E0) Z(K) = Z(K)/A(K,K)
         IF (A(K,K) .EQ. 0.0E0) Z(K) = 1.0E0
         T = -Z(K)
         CALL SAXPY(K-1,T,A(1,K),1,Z(1),1)
  160 CONTINUE
C     MAKE ZNORM = 1.0
      S = 1.0E0/SASUM(N,Z,1)
      CALL SSCAL(N,S,Z,1)
      YNORM = S*YNORM
C
      IF (ANORM .NE. 0.0E0) RCOND = YNORM/ANORM
      IF (ANORM .EQ. 0.0E0) RCOND = 0.0E0
      RETURN
      END
c///////////////////////////////////////////////////////////////////////////
c---------------------------------------------------------------------------
      SUBROUTINE SGESL(A,LDA,N,IPVT,B,JOB)
c---------------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER LDA,N,IPVT(1),JOB
      REAL A(LDA,1),B(1)
C
C     SGESL SOLVES THE REAL SYSTEM
C     A * X = B  OR  TRANS(A) * X = B
C     USING THE FACTORS COMPUTED BY SGECO OR SGEFA.
C
C     ON ENTRY
C
C        A       REAL(LDA, N)
C                THE OUTPUT FROM SGECO OR SGEFA.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY  A .
C
C        N       INTEGER
C                THE ORDER OF THE MATRIX  A .
C
C        IPVT    INTEGER(N)
C                THE PIVOT VECTOR FROM SGECO OR SGEFA.
C
C        B       REAL(N)
C                THE RIGHT HAND SIDE VECTOR.
C
C        JOB     INTEGER
C                = 0         TO SOLVE  A*X = B ,
C                = NONZERO   TO SOLVE  TRANS(A)*X = B  WHERE
C                            TRANS(A)  IS THE TRANSPOSE.
C
C     ON RETURN
C
C        B       THE SOLUTION VECTOR  X .
C
C     ERROR CONDITION
C
C        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
C        ZERO ON THE DIAGONAL.  TECHNICALLY THIS INDICATES SINGULARITY
C        BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
C        SETTING OF LDA .  IT WILL NOT OCCUR IF THE SUBROUTINES ARE
C        CALLED CORRECTLY AND IF SGECO HAS SET RCOND .GT. 0.0
C        OR SGEFA HAS SET INFO .EQ. 0 .
C
C     TO COMPUTE  INVERSE(A) * C  WHERE  C  IS A MATRIX
C     WITH  P  COLUMNS
C           CALL SGECO(A,LDA,N,IPVT,RCOND,Z)
C           IF (RCOND IS TOO SMALL) GO TO ...
C           DO 10 J = 1, P
C              CALL SGESL(A,LDA,N,IPVT,C(1,J),0)
C        10 CONTINUE
C
C     LINPACK. THIS VERSION DATED 08/14/78 .
C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C     SUBROUTINES AND FUNCTIONS
C
C     BLAS SAXPY,SDOT
C
C     INTERNAL VARIABLES
C
      REAL SDOT,T
      INTEGER K,KB,L,NM1
C
      NM1 = N - 1
      IF (JOB .NE. 0) GO TO 50
C
C        JOB = 0 , SOLVE  A * X = B
C        FIRST SOLVE  L*Y = B
C
         IF (NM1 .LT. 1) GO TO 30
         DO 20 K = 1, NM1
            L = IPVT(K)
            T = B(L)
            IF (L .EQ. K) GO TO 10
               B(L) = B(K)
               B(K) = T
   10       CONTINUE
            CALL SAXPY(N-K,T,A(K+1,K),1,B(K+1),1)
   20    CONTINUE
   30    CONTINUE
C
C        NOW SOLVE  U*X = Y
C
         DO 40 KB = 1, N
            K = N + 1 - KB
            B(K) = B(K)/A(K,K)
            T = -B(K)
            CALL SAXPY(K-1,T,A(1,K),1,B(1),1)
   40    CONTINUE
      GO TO 100
   50 CONTINUE
C
C        JOB = NONZERO, SOLVE  TRANS(A) * X = B
C        FIRST SOLVE  TRANS(U)*Y = B
C
         DO 60 K = 1, N
            T = SDOT(K-1,A(1,K),1,B(1),1)
            B(K) = (B(K) - T)/A(K,K)
   60    CONTINUE
C
C        NOW SOLVE TRANS(L)*X = Y
C
         IF (NM1 .LT. 1) GO TO 90
         DO 80 KB = 1, NM1
            K = N - KB
            B(K) = B(K) + SDOT(N-K,A(K+1,K),1,B(K+1),1)
            L = IPVT(K)
            IF (L .EQ. K) GO TO 70
               T = B(L)
               B(L) = B(K)
               B(K) = T
   70       CONTINUE
   80    CONTINUE
   90    CONTINUE
  100 CONTINUE
      RETURN
      END
      SUBROUTINE SAXPY(n,sa,sx,incx,sy,incy)
      IMPLICIT NONE
c
c     constant times a vector plus a vector.
c     uses unrolled loop for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      real sx(1),sy(1),sa
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if (sa .eq. 0.0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        sy(iy) = sy(iy) + sa*sx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        sy(i) = sy(i) + sa*sx(i)
   30 continue
      if( n .lt. 4 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,4
        sy(i) = sy(i) + sa*sx(i)
        sy(i + 1) = sy(i + 1) + sa*sx(i + 1)
        sy(i + 2) = sy(i + 2) + sa*sx(i + 2)
        sy(i + 3) = sy(i + 3) + sa*sx(i + 3)
   50 continue
      return
      end
      real function sdot(n,sx,incx,sy,incy)
      IMPLICIT NONE
c
c     forms the dot product of two vectors.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      real sx(1),sy(1),stemp
      integer i,incx,incy,ix,iy,m,mp1,n
c
      stemp = 0.0e0
      sdot = 0.0e0
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        stemp = stemp + sx(ix)*sy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      sdot = stemp
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        stemp = stemp + sx(i)*sy(i)
   30 continue
      if( n .lt. 5 ) go to 60
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) +
     *   sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4)
   50 continue
   60 sdot = stemp
      return
      end
      SUBROUTINE SSCAL(n,sa,sx,incx)
      IMPLICIT NONE
c
c     scales a vector by a constant.
c     uses unrolled loops for increment equal to 1.
c     jack dongarra, linpack, 3/11/78.
c
      real sa,sx(1)
      integer i,incx,m,mp1,n,nincx
c
      if(n.le.0)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      nincx = n*incx
      do 10 i = 1,nincx,incx
        sx(i) = sa*sx(i)
   10 continue
      return
c
c        code for increment equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        sx(i) = sa*sx(i)
   30 continue
      if( n .lt. 5 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        sx(i) = sa*sx(i)
        sx(i + 1) = sa*sx(i + 1)
        sx(i + 2) = sa*sx(i + 2)
        sx(i + 3) = sa*sx(i + 3)
        sx(i + 4) = sa*sx(i + 4)
   50 continue
      return
      end
      real function sasum(n,sx,incx)
      IMPLICIT NONE
c
c     takes the sum of the absolute values.
c     uses unrolled loops for increment equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      real sx(1),stemp
      integer i,incx,m,mp1,n,nincx
c
      sasum = 0.0e0
      stemp = 0.0e0
      if(n.le.0)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      nincx = n*incx
      do 10 i = 1,nincx,incx
        stemp = stemp + abs(sx(i))
   10 continue
      sasum = stemp
      return
c
c        code for increment equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,6)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        stemp = stemp + abs(sx(i))
   30 continue
      if( n .lt. 6 ) go to 60
   40 mp1 = m + 1
      do 50 i = mp1,n,6
        stemp = stemp + abs(sx(i)) + abs(sx(i + 1)) + abs(sx(i + 2))
     *  + abs(sx(i + 3)) + abs(sx(i + 4)) + abs(sx(i + 5))
   50 continue
   60 sasum = stemp
      return
      end
      SUBROUTINE SGEFA(A,LDA,N,IPVT,INFO)
      IMPLICIT NONE
      INTEGER LDA,N,IPVT(1),INFO
      REAL A(LDA,1)
C
C     SGEFA FACTORS A REAL MATRIX BY GAUSSIAN ELIMINATION.
C
C     SGEFA IS USUALLY CALLED BY SGECO, BUT IT CAN BE CALLED
C     DIRECTLY WITH A SAVING IN TIME IF  RCOND  IS NOT NEEDED.
C     (TIME FOR SGECO) = (1 + 9/N)*(TIME FOR SGEFA) .
C
C     ON ENTRY
C
C        A       REAL(LDA, N)
C                THE MATRIX TO BE FACTORED.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY  A .
C
C        N       INTEGER
C                THE ORDER OF THE MATRIX  A .
C
C     ON RETURN
C
C        A       AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
C                WHICH WERE USED TO OBTAIN IT.
C                THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE
C                L  IS A PRODUCT OF PERMUTATION AND UNIT LOWER
C                TRIANGULAR MATRICES AND  U  IS UPPER TRIANGULAR.
C
C        IPVT    INTEGER(N)
C                AN INTEGER VECTOR OF PIVOT INDICES.
C
C        INFO    INTEGER
C                = 0  NORMAL VALUE.
C                = K  IF  U(K,K) .EQ. 0.0 .  THIS IS NOT AN ERROR
C                     CONDITION FOR THIS SUBROUTINE, BUT IT DOES
C                     INDICATE THAT SGESL OR SGEDI WILL DIVIDE BY ZERO
C                     IF CALLED.  USE  RCOND  IN SGECO FOR A RELIABLE
C                     INDICATION OF SINGULARITY.
C
C     LINPACK. THIS VERSION DATED 08/14/78 .
C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C     SUBROUTINES AND FUNCTIONS
C
C     BLAS SAXPY,SSCAL,ISAMAX
C
C     INTERNAL VARIABLES
C
      REAL T
      INTEGER ISAMAX,J,K,KP1,L,NM1
C
C
C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
C
      INFO = 0
      NM1 = N - 1
      IF (NM1 .LT. 1) GO TO 70
      DO 60 K = 1, NM1
         KP1 = K + 1
C
C        FIND L = PIVOT INDEX
C
         L = ISAMAX(N-K+1,A(K,K),1) + K - 1
         IPVT(K) = L
C
C        ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
C
         IF (A(L,K) .EQ. 0.0E0) GO TO 40
C
C           INTERCHANGE IF NECESSARY
C
            IF (L .EQ. K) GO TO 10
               T = A(L,K)
               A(L,K) = A(K,K)
               A(K,K) = T
   10       CONTINUE
C
C           COMPUTE MULTIPLIERS
C
            T = -1.0E0/A(K,K)
            CALL SSCAL(N-K,T,A(K+1,K),1)
C
C           ROW ELIMINATION WITH COLUMN INDEXING
C
            DO 30 J = KP1, N
               T = A(L,J)
               IF (L .EQ. K) GO TO 20
                  A(L,J) = A(K,J)
                  A(K,J) = T
   20          CONTINUE
               CALL SAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
   30       CONTINUE
         GO TO 50
   40    CONTINUE
            INFO = K
   50    CONTINUE
   60 CONTINUE
   70 CONTINUE
      IPVT(N) = N
      IF (A(N,N) .EQ. 0.0E0) INFO = N
      RETURN
      END
      integer function isamax(n,sx,incx)
      IMPLICIT NONE
c
c     finds the index of element having max. absolute value.
c     jack dongarra, linpack, 3/11/78.
c
      real sx(1),smax
      integer i,incx,ix,n
c
      isamax = 0
      if( n .lt. 1 ) return
      isamax = 1
      if(n.eq.1)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      ix = 1
      smax = abs(sx(1))
      ix = ix + incx
      do 10 i = 2,n
         if(abs(sx(ix)).le.smax) go to 5
         isamax = i
         smax = abs(sx(ix))
    5    ix = ix + incx
   10 continue
      return
c
c        code for increment equal to 1
c
   20 smax = abs(sx(1))
      do 30 i = 2,n
         if(abs(sx(i)).le.smax) go to 30
         isamax = i
         smax = abs(sx(i))
   30 continue
      return
      end

max_stars_repo_name: NoiseCIEI/RayTomo

max_issues_repo_name: NoiseCIEI/RayTomo

max_forks_repo_name: NoiseCIEI/RayTomo


text:
module heatbalance

  use slf_cgs
  use globals, only: kappabp, kappsct, miu
  use iso_fortran_env, only: r64 => real64
  implicit none

  private :: ramp

contains

  !----------------------------------------------------------------------------!
  ! THIS IS DEPRECATED! do not use! use heatbil2 instead

  pure subroutine heatbil(T, Trad, Pgas, heat)

    real(r64), intent(inout) :: T
    real(r64), intent(in) :: Trad, Pgas, heat
    real(r64), dimension(3,2) :: kap
    real(r64), dimension(2,2) :: F
    real(r64) :: rho, bil, bil_T, rhotemp
    integer :: i
    integer, parameter :: niter = 64

    rhotemp = Pgas * miu * cgs_mhydr / cgs_boltz

    converge: do i = 1,niter

      rho = rhotemp / T

      call kappabp(rho, T, kap(1,1), kap(2,1), kap(3,1))
      call kappsct(rho, T, kap(1,2), kap(2,2), kap(3,2))

      F(1,:) = kap(1,:)
      F(2,:) = kap(3,:) - rho / T * kap(2,:)

      bil = (-T*heat + 4*cgs_stef*rhotemp*(T - Trad)*(4*Trad**4* &
      cgs_k_over_mec2*F(1, 2) + (T + Trad)*(T**2 + Trad**2)*F(1, 1)))/( &
      T*Trad**4)
      bil_T = 4*cgs_stef*rhotemp*(T**5*F(2, 1) + 3*T**4*F(1, 1) + 4*T**2*Trad &
            **4*cgs_k_over_mec2*F(2, 2) - 4*T*Trad**5*cgs_k_over_mec2*F(2, 2 &
            ) - T*Trad**4*F(2, 1) + 4*Trad**5*cgs_k_over_mec2*F(1, 2) + Trad &
            **4*F(1, 1))/(T**2*Trad**4)

      T = T - ramp(i,niter) * bil / bil_T

    end do converge

  end subroutine

  !----------------------------------------------------------------------------!
  ! solves the heating-cooling balance equation

  subroutine heatbil2(rho, tgas, trad, heat, isobar)
    real(r64), intent(in) :: trad, heat
    real(r64), intent(inout) :: rho, tgas
    logical, intent(in) :: isobar
    real(r64) :: rhotemp, kabpv(3), ksctv(3), cool, cool_dT, cool_dr
    integer :: i
    integer, parameter :: niter = 24

    rhotemp = rho * tgas

    tgas = heat  &
          & / (16 * cgs_k_over_mec2 * cgs_kapes * rho &
          & * cgs_stef * trad**4 )
    tgas = sqrt(trad**2 + tgas**2)

    do i = 1,niter

      call kappsct(rho, tgas, ksctv(1), ksctv(2), ksctv(3))
      call kappabp(rho, tgas, kabpv(1), kabpv(2), kabpv(3))

      cool = 4*cgs_stef*rho*(4*cgs_k_over_mec2*trad**4*(tgas - trad)*ksctv(1) &
            + (tgas**4 - trad**4)*kabpv(1))

      cool_dT = 4*cgs_stef*rho*(4*cgs_k_over_mec2*trad**4*(tgas - trad)*ksctv( &
            3) + 4*cgs_k_over_mec2*trad**4*ksctv(1) + 4*tgas**3*kabpv(1) + ( &
            tgas**4 - trad**4)*kabpv(3))

      if (isobar) then
        cool_dr = 4*cgs_stef*(4*cgs_k_over_mec2*trad**4*(tgas - trad)*ksctv(1) &
            + rho*(4*cgs_k_over_mec2*trad**4*(tgas - trad)*ksctv(2) + (tgas**4 &
            - trad**4)*kabpv(2)) + (tgas**4 - trad**4)*kabpv(1))
        cool_dT = cool_dT - rho / tgas * cool_dr
      end if

      tgas = tgas - (cool - heat) / cool_dT * ramp(i,niter)
      if (isobar) rho = rhotemp / tgas

    end do

  end subroutine heatbil2

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

  elemental function ramp(i,n) result(y)
    integer, intent(in) :: i,n
    real(r64) :: t,y
    t = merge(real(i) / n, 1.0, i .le. n)
    y = (3*t**2 - 8*t + 6) * t**2
  end function

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

end module

max_stars_repo_name: gronki/pydiskvert

max_issues_repo_name: gronki/pydiskvert

max_forks_repo_name: gronki/pydiskvert


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! Initialization File
!
! Generated by KPP-2.1 symbolic chemistry Kinetics PreProcessor
!       (http://www.cs.vt.edu/~asandu/Software/KPP)
! KPP is distributed under GPL, the general public licence
!       (http://www.gnu.org/copyleft/gpl.html)
! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa
! (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech
!     With important contributions from:
!        M. Damian, Villanova University, USA
!        R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany
!
! File                 : racm_esrlsorg_Initialize.f90
! Time                 : Mon Feb 27 11:08:48 2017
! Working directory    : /data/ksetigui/setigui/WRFV3_171015/chem/KPP/mechanisms/racm_esrlsorg
! Equation file        : racm_esrlsorg.kpp
! Output root filename : racm_esrlsorg
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

MODULE racm_esrlsorg_Initialize

  IMPLICIT NONE

CONTAINS

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! Initialize - function to initialize concentrations
!   Arguments :
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE Initialize ( )

  USE racm_esrlsorg_Global

  INTEGER :: i
  REAL(kind=dp) :: x

  CFACTOR = 1.000000e+00_dp

  x = (0.)*CFACTOR
  DO i = 1, NVAR
    VAR(i) = x
  END DO

  x = (0.)*CFACTOR
  DO i = 1, NFIX
    FIX(i) = x
  END DO

! constant rate coefficients
! END constant rate coefficients

! INLINED initializations

! End INLINED initializations

END SUBROUTINE Initialize

! End of Initialize function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

END MODULE racm_esrlsorg_Initialize


text:
Dos Pinos Housing Cooperative is a limited equity housing cooperative.  They are governed by an elected board of directors and use a management company for daytoday operations.  Owned by the residents, they are shareholders in the notforprofit corporation that holds a blanket mortgage over the entire property. The notforprofit cooperative corporation  in turn grants leases to each shareholder, giving them the right to live in a home.

There are 60 housing units at Dos Pinos, which include 1 bedroom, 2 bedroom, and 3 bedroom units.  All units are attached townhomes, organized into 9 buildings on a shared property. Member households are diverse, including single people, couples, and families with a wide range of ages and backgrounds. The property is located at 2550 Sycamore Lane, and is surrounded by residential housing on 3 sides and the greenbelt on the north side.

Dos Pinos Housing Cooperative offers a form of affordable homeownership in a cooperative community setting. It can be understood as a hybrid between homeownership and renting. Members join the coop by purchasing a share in the cooperative corporation, which gives them the right to live in a unit as well as voting rights within the cooperative corporation.  The cooperative corporation then leases a unit to the member. Members must pay a monthly carrying charge or assessment which includes overall community costs such as: water, sewage, management fees, maintenence, landscaping, and a share of the overall blanket mortgage.  This monthly assesssment is considerably less than rent at a comparablesized condo or apartment in Davis. The amount of the monthly assessment is set by the coop board of directors (which is elected from the membership of the coop). When a member decides to leave, they market and sell their unit for a set price.  This set price is based on the original purchase price, plus the value of interest and improvements. This provides a steady return on investment for member households.

Historically, there has been a one to 3 year waiting list, which is the case as of May 2012.

A flyer for a 2 bedroom unit, which was on sale in February 2011, is attached here File(DosPinos3aflier.pdf).

History

The Davis City Council first approved the construction in 1984, with the first occupants moving in during the fall of 1985.  The members of the cooperative elected the first board of directors and took over ownership of the initial 60 units in March 1986.

To learn more about rental housing in Davis, check out our Housing Guide as well as our Apartments pages.
To learn more about homeownership in Davis, check out the Real Estate page.
To learn more about cooperative housing in Davis, check out Cooperative Housing.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program testsaxpy_omp45_f
  implicit none
  integer, parameter :: n = ishft(1,21)
  integer :: i

  real,allocatable :: x(:), y(:)
  real :: a

  allocate( x(n), y(n) )
  x = 1.0
  y = 2.0
  a = 2.0

  !$omp target data map(to:x) map(tofrom:y) use_device_ptr(x,y)
    call testsaxpy_cudafortran(x,y,a,n)
  !$omp end target data

  write(*,*) "Ran FORTRAN OMP45 kernel. Max error: ", maxval(abs(y-4.0))
end program testsaxpy_omp45_f

max_stars_repo_name: bblakeley/FGPU

max_issues_repo_name: bblakeley/FGPU

max_forks_repo_name: bblakeley/FGPU


!===-- module/iso_c_binding.f90 --------------------------------------------===!
!
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!
!------------------------------------------------------------------------------!

! See Fortran 2018, clause 18.2

module iso_c_binding

  use __Fortran_builtins, only: &
    c_f_pointer => __builtin_c_f_pointer, &
    c_ptr => __builtin_c_ptr, &
    c_funptr => __builtin_c_funptr

  type(c_ptr), parameter :: c_null_ptr = c_ptr()
  type(c_funptr), parameter :: c_null_funptr = c_funptr()

  ! Table 18.2 (in clause 18.3.1)
  ! TODO: Specialize (via macros?) for alternative targets
  integer, parameter :: &
    c_int8_t = 1, &
    c_int16_t = 2, &
    c_int32_t = 4, &
    c_int64_t = 8, &
    c_int128_t = 16 ! anticipating future addition
  integer, parameter :: &
    c_int = c_int32_t, &
    c_short = c_int16_t, &
    c_long = c_int64_t, &
    c_long_long = c_int64_t, &
    c_signed_char = c_int8_t, &
    c_size_t = c_long_long, &
    c_intmax_t = c_int128_t, &
    c_intptr_t = c_size_t, &
    c_ptrdiff_t = c_size_t
  integer, parameter :: &
    c_int_least8_t = c_int8_t, &
    c_int_fast8_t = c_int8_t, &
    c_int_least16_t = c_int16_t, &
    c_int_fast16_t = c_int16_t, &
    c_int_least32_t = c_int32_t, &
    c_int_fast32_t = c_int32_t, &
    c_int_least64_t = c_int64_t, &
    c_int_fast64_t = c_int64_t, &
    c_int_least128_t = c_int128_t, &
    c_int_fast128_t = c_int128_t

  integer, parameter :: &
    c_float = 4, &
    c_double = 8, &
#if __x86_64__
    c_long_double = 10
#else
    c_long_double = 16
#endif

  integer, parameter :: &
    c_float_complex = c_float, &
    c_double_complex = c_double, &
    c_long_double_complex = c_long_double

  integer, parameter :: c_bool = 1 ! TODO: or default LOGICAL?
  integer, parameter :: c_char = 1

  ! C characters with special semantics
  character(kind=c_char, len=1), parameter :: c_null_char = achar(0)
  character(kind=c_char, len=1), parameter :: c_alert = achar(7)
  character(kind=c_char, len=1), parameter :: c_backspace = achar(8)
  character(kind=c_char, len=1), parameter :: c_form_feed = achar(12)
  character(kind=c_char, len=1), parameter :: c_new_line = achar(10)
  character(kind=c_char, len=1), parameter :: c_carriage_return = achar(13)
  character(kind=c_char, len=1), parameter :: c_horizontal_tab = achar(9)
  character(kind=c_char, len=1), parameter :: c_vertical_tab =  achar(11)

 contains

  logical function c_associated(c_ptr_1, c_ptr_2)
    type(c_ptr), intent(in) :: c_ptr_1
    type(c_ptr), intent(in), optional :: c_ptr_2
    if (c_ptr_1%__address == c_null_ptr%__address) then
      c_associated = .false.
    else if (present(c_ptr_2)) then
      c_associated = c_ptr_1%__address == c_ptr_2%__address
    else
      c_associated = .true.
    end if
  end function c_associated

  function c_loc(x)
    type(c_ptr) :: c_loc
    type(*), intent(in) :: x
    c_loc = c_ptr(loc(x))
  end function c_loc

  function c_funloc(x)
    type(c_funptr) :: c_funloc
    external :: x
    c_funloc = c_funptr(loc(x))
  end function c_funloc

  ! TODO c_f_procpointer
  ! TODO c_sizeof

end module iso_c_binding


text:
      subroutine readlut(csat_id,csat_typ,maxch,nch,
     &chtype,nx,ny,ri,rj,istatus)
c
      implicit none

      integer nx,ny
      integer maxch
      integer nch
      integer istatus
      integer i,n
      integer lend,lenf
      integer ispec
      integer istat

      logical   lgot_lut(maxch)

      real    ri(nx,ny,maxch)
      real    rj(nx,ny,maxch)
c     real    ri_in(nx,ny)
c     real    rj_in(nx,ny)
c     real    rdummy(nx,ny)

      character*100 cpath
      character*255 file
      character*6   csat_id
      character*3   csat_typ
      character*3   chtype(maxch)
      character*3   ct
c
c-----------------------------------------------------------------------
c
      istatus = 1

      call get_directory('static',cpath,lend)
      cpath=cpath(1:lend)//'lvd/'
      lend=index(cpath,' ')-1

c JS: 8-31-06. No longer reading a mapping look-up-table (lut).
c
      do i=1,maxch
         lgot_lut(i)=.true.   !!!.false.
      enddo

      do i=1,nch

         call lvd_file_specifier(chtype(i),ispec,istat)
         if(istat.eq.0)then
c           if(.not.lgot_lut(ispec))then
c              lgot_lut(ispec)=.true.
c              ct=chtype(i)
               if(ispec.eq.2.or.ispec.eq.4.or.ispec.eq.5)then  !all ir channels have the same ri/rj values
                  call move(rj(1,1,2),rj(1,1,ispec),nx,ny)
                  call move(ri(1,1,2),ri(1,1,ispec),nx,ny)
c                 ct='ir'
               endif

c              n=index(ct,' ')-1
c              if(n.le.0)n=3
c              file=cpath(1:lend)//csat_id//'-llij-'
c              lenf=index(file,' ')-1
c              file=file(1:lenf)//ct(1:n)//'-'//csat_typ//'.lut'
c              n=index(file,' ')-1
c              open(12,file=file,
c    &form='unformatted',status='old',err=101)
c              write(6,*)'Reading ',file(1:n)
c              read(12,err=23,end=23) rdummy
c              read(12,err=23,end=23) rdummy
c              read(12,err=23,end=23) ri_in
c              read(12,err=23,end=23) rj_in
c              close (12)
c              call move(ri_in,ri(1,1,ispec),nx,ny)
c              call move(rj_in,rj(1,1,ispec),nx,ny)
c           else
c              print*,'Not reading mapping look-up-table'
c              print*,'Return without ri/rj values'
c           endif

         endif
      enddo

      istatus = 0
      goto 1000

23    write(6,*)'Error reading or eof ll/ij lookup table'
      close(12)
      goto 1000

101   write(6,*)'Error opening file ',file(1:n)

1000  return
      end
c
c================================================
c
      subroutine check_luts(cfname_cur,isat,jtype,
     &chtype,maxchannels,nchannels,l_lut_flag,istatus)
c
c most applications will relocalize which will delete any
c existing lut's and new ones are automatically generated.
c For WFO and AFWA, the satellite data files can possibly
c change. This is particularly true of AFWA satellite data.
c
      implicit none

      integer        maxchannels
      integer        nchannels

      character      cfname_cur*(*)
      character      chtype(maxchannels)*3
      character      cname*100
      character      c_afwa_fname*100
      character      cfname*200

      integer        i,il
      integer        isat,jtype,kch
      integer        istatus
      integer        nlin,npix
      integer        nx,ny
      integer        ispec
      integer        ilat00,ilon00
      integer        i_la1,i_lo1
      integer        idx,idy,iresx,iresy
      integer        nw_line,nw_pix
      integer        nxaf,nyaf
      integer        lenf

      integer        strpix,strscnl,stppix,stpscnl
      integer        reqobstm,imgtype
      integer        iwidth,idepth
      integer        istrbdy1,istrbdy2,istpbdy1,istpbdy2
      integer        bepixfc,bescnfc,fsci,idecimat

      integer        istrtline,istrtpix
      integer        nw_vis_line_gwc
      integer        nw_vis_pix_gwc

      logical        l_lut_flag

      real resx,resy
      real centerlat,centerlon
      real rlat00,rlon00
      real rlatnxny,rlonnxny
      real rlatdxdy,rlondxdy
      real dx,dy
      real golatsbp,golonsbp,goalpha

      include 'satellite_dims_lvd.inc'
      include 'satellite_common_lvd.inc'

      istatus = 1
      l_lut_flag=.false.

      do i = 1,nchannels
      call lvd_file_specifier(chtype(i),ispec,istatus)

      if(ispec.eq.1)then

         resx=r_resolution_x_vis(jtype,isat)
         resy=r_resolution_y_vis(jtype,isat)
         nlin=n_lines_vis(jtype,isat)
         npix=n_pixels_vis(jtype,isat)
         nw_line=i_nwline_vis(jtype,isat)
         nw_pix=i_nwpix_vis(jtype,isat)
      elseif(ispec.eq.2.or.ispec.eq.4.or.ispec.eq.5)then

         resx=r_resolution_x_ir(jtype,isat)
         resy=r_resolution_y_ir(jtype,isat)
         nlin=n_lines_ir(jtype,isat)
         npix=n_pixels_ir(jtype,isat)
         nw_line=i_nwline_ir(jtype,isat)
         nw_pix=i_nwpix_ir(jtype,isat)

      elseif(ispec.eq.3)then

         resx=r_resolution_x_wv(jtype,isat)
         resy=r_resolution_y_wv(jtype,isat)
         nlin=n_lines_wv(jtype,isat)
         npix=n_pixels_wv(jtype,isat)
         nw_line=i_nwline_wv(jtype,isat)
         nw_pix=i_nwpix_wv(jtype,isat)

      endif

c check if namelist parameters are current

c
c --- WFO ---
c
      if(c_sat_types(jtype,isat).eq.'wfo')then

         call get_wfo_nav_parms(path_to_raw_sat(ispec,jtype,isat)
     &,chtype(i),centerlat,centerlon,rlat00,rlon00,rlatnxny,rlonnxny,
     &rlatdxdy,rlondxdy,dx,dy,nx,ny,istatus)

         if(istatus.eq.0)then

           ilat00=nint(rlat00*1000.)
           ilon00=nint(rlon00*1000.)
           i_la1 =nint(r_la1(jtype,isat)*1000.)
           i_lo1 =nint(r_lo1(jtype,isat)*1000.)
           iresx =nint(resx)
           iresy =nint(resy)
           idx   =nint(dx)
           idy   =nint(dy)

           if((ilat00.ne.i_la1).or.
     &        (ilon00.ne.i_lo1) )then
             write(6,*)'rlat00/rlon00/r_la1/rlo1 ',rlat00,rlon00,
     &r_la1(jtype,isat),r_lo1(jtype,isat)
              l_lut_flag=.true.
           endif
           if( idx.ne.iresx .or. idy.ne.iresy )then
              write(6,*)'dx/dy/resx/resy/ ',dx,dy,resx,resy
              l_lut_flag=.true.
           endif
           if( nlin.ne.ny .or. npix.ne.nx )then
              write(6,*)'nx/ny/npix/nlin ',nx,ny,npix,nlin
              l_lut_flag=.true.
           endif

         endif

      elseif(c_sat_types(jtype,isat).eq.'gwc')then
c
c --- AFWA ---
c for GOES data only
         cname=c_afwa_fname(c_sat_id(isat),chtype(i))
         call s_len(cname,il)
         lenf=index(path_to_raw_sat(ispec,jtype,isat),' ')-1
         cfname=path_to_raw_sat(ispec,jtype,isat)(1:lenf)//cname(1:il)

         call read_gwc_header(cfname,l_cell_afwa,
     &strpix,strscnl,stppix,stpscnl,
     &reqobstm,imgtype,golatsbp,golonsbp,iwidth,idepth,goalpha,istrbdy1,
     &istrbdy2,istpbdy1,istpbdy2,bepixfc,bescnfc,fsci,idecimat,istatus)

         if(istatus.eq.0)then
            if(l_cell_afwa)then
               if(iwidth*256.ne.npix.or.idepth*64.ne.nlin)then
                  l_lut_flag=.true.
               endif
            else
               nxaf=stppix-strpix+1
               nyaf=stpscnl-strscnl+1
               if(nxaf.ne.npix.or.nyaf.ne.nlin)then
                  l_lut_flag=.true.
               endif
            endif
c
c NOTE: nw_line and nw_pix are hardwired by src/include/satdata_lvd.f
c       however, gen_lut_gvar.f uses the new values for these variables.
c       this will always trip the lut regeneration unless nw_line/pix match.
c
            istrtline=nw_vis_line_gwc(chtype(i),idecimat,bescnfc,fsci)
            istrtpix=nw_vis_pix_gwc(chtype(i),idecimat,bepixfc,goalpha)
            if(istrtline.ne.nw_line.or.istrtpix.ne.nw_pix)then
               l_lut_flag=.true.
            endif

         else
            print*,'gwc header not read. No lut update'
         endif

      elseif(c_sat_types(jtype,isat).eq.'gvr')then

         if(resx.eq.0.0.or.resy.eq.0.0.or.
     &      nlin.eq.0.0.or.npix.eq.0.0)then

            l_lut_flag = .true.
         endif

      endif

      enddo
c
c this check is for gvar type data only. once a day using
c a current O&A block.
c
      if(cfname_cur(6:8).eq.'000'.and.
     &   (c_sat_types(jtype,isat).eq.'gvr'   .or.
     &    c_sat_types(jtype,isat).eq.'gwc')       )then
         l_lut_flag=.true.
c        write(6,*)'Auto-update the gvar navigation'
      endif

c 8-31-06: JS. forces an a-ok return and no regeneration of
c mapping. All mapping now done on the fly.

      l_lut_flag=.false.

      return
      end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


program pathlib_test

use, intrinsic :: iso_fortran_env, only : stderr=>error_unit
use pathlib, only : get_filename, mkdir, expanduser, is_absolute, make_absolute, directory_exists

implicit none (type, external)

call test_get_filename()

call test_expanduser_absolute()

call test_directory_exists()

contains

subroutine test_directory_exists()

integer :: i

if(.not.(directory_exists('.'))) error stop "did not detect '.' as directory"

open(newunit=i, file='test-pathlib.h5', status='replace')
close(i)
if((directory_exists('test-pathlib.h5'))) error stop "detected file as directory"
call unlink('test-pathlib.h5')

print *," OK: pathlib: directory_exists"
end subroutine test_directory_exists

subroutine test_get_filename()

character(:), allocatable:: fn
integer :: i
logical :: e

if(get_filename(' ') /= '') error stop 'empty 1'
if(get_filename(' ',' ') /= '') error stop 'empty 2'
!! " " instead of "" to avoid compile-time glitch error with GCC-10 with -Og

call unlink('test-pathlib.h5')

if(len(get_filename('test-pathlib.h5')) > 0) error stop 'not exist full 1'

fn = get_filename('test-pathlib')
if(len(fn) > 0) then
  write(stderr,*) 'ERROR: ',fn, len(fn)
  error stop 'not exist stem 1'
endif

!> touch empty file
open(newunit=i, file='test-pathlib.h5', status='replace')
close(i)
inquire(file='test-pathlib.h5', exist=e)
if(.not.e) error stop 'could not create test-pathlib.h5'

if(get_filename('test-pathlib.h5') /= 'test-pathlib.h5') error stop 'exist full 1'
if(get_filename('test-pathlib') /= 'test-pathlib.h5') error stop 'exist stem 1'

fn = get_filename('.', 'test-pathlib')
if(fn /= './test-pathlib.h5') error stop 'exist stem 2: ' // fn

fn = get_filename('./test-pathlib', 'test-pathlib')
if(fn /= './test-pathlib.h5') error stop 'exist parts 2: ' // fn

fn = get_filename('./test-pathlib.h5', 'test-pathlib')
if(fn /= './test-pathlib.h5') error stop 'exist full 2: ' // fn

call unlink('test-pathlib.h5')

open(newunit=i, file='test-pathlib.nc', status='replace')
close(i)

if(get_filename('test-pathlib.nc') /= 'test-pathlib.nc') error stop 'exist full 1a'
if(get_filename('test-pathlib') /= 'test-pathlib.nc') error stop 'exist stem 1a'
if(get_filename('.', 'test-pathlib') /= './test-pathlib.nc') error stop 'exist stem 2a'
if(get_filename('./test-pathlib', 'test-pathlib') /= './test-pathlib.nc') error stop 'exist parts 2a'

call unlink('test-pathlib.nc')

call mkdir('temp1/temp2')
call unlink('temp1/temp2/test-pathlib.h5')
fn = get_filename('temp1/temp2', 'test-pathlib')
if (fn /= '') error stop 'non-exist dir'
open(newunit=i, file='temp1/temp2/test-pathlib.h5', status='replace')
close(i)
fn = get_filename('temp1/temp2', 'test-pathlib')
if (fn /= 'temp1/temp2/test-pathlib.h5') error stop 'exist dir full 2'

fn = get_filename('./temp1/temp2', 'test-pathlib')
if (fn /= './temp1/temp2/test-pathlib.h5') error stop 'exist dir full 2a'

end subroutine test_get_filename

subroutine test_expanduser_absolute()

character(:), allocatable:: fn
character(16) :: fn2

fn = expanduser("~")
if (fn(1:1) == "/") then
  if (.not.is_absolute("/")) error stop "is_absolute('/') on Unix should be true"
  if (is_absolute("c:/")) error stop "is_absolute('c:/') on Unix should be false"

  fn2 = make_absolute("rel", "/foo")
  if (fn2 /= "/foo/rel") error stop "did not make_absolute Unix /foo/rel, got: " // fn2
else
  if (.not.is_absolute("J:/")) error stop "is_absolute('J:/') on Windows should be true"
  if (.not.is_absolute("j:/")) error stop "is_absolute('j:/') on Windows should be true"
  if (is_absolute("/")) error stop "is_absolute('/') on Windows should be false"

  fn2 = make_absolute("rel", "j:/foo")
  if (fn2 /= "j:/foo/rel") error stop "did not make_absolute Windows j:/foo/rel, got: " // fn2
endif

print *, "OK: pathlib: expanduser,is_absolute"
end subroutine test_expanduser_absolute

subroutine unlink(path)
character(*), intent(in) :: path
integer :: i
logical :: e

inquire(file=path, exist=e)
if (.not.e) return

open(newunit=i, file=path, status='old')
close(i, status='delete')
end subroutine unlink

end program


text:
cdis    Forecast Systems Laboratory
cdis    NOAA/OAR/ERL/FSL
cdis    325 Broadway
cdis    Boulder, CO     80303
cdis
cdis    Forecast Research Division
cdis    Local Analysis and Prediction Branch
cdis    LAPS
cdis
cdis    This software and its documentation are in the public domain and
cdis    are furnished "as is."  The United States government, its
cdis    instrumentalities, officers, employees, and agents make no
cdis    warranty, express or implied, as to the usefulness of the software
cdis    and documentation for any purpose.  They assume no responsibility
cdis    (1) for the use of the software and documentation; or (2) to provide
cdis     technical support to users.
cdis
cdis    Permission to use, copy, modify, and distribute this software is
cdis    hereby granted, provided that the entire disclaimer notice appears
cdis    in all copies.  All modifications to this software must be clearly
cdis    documented, and are solely the responsibility of the agent making
cdis    the modifications.  If significant modifications or enhancements
cdis    are made to this software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis
cdis
cdis
cdis
cdis
cdis

       subroutine NOWRADWSI_to_LAPS(ctype,
     &                    filename,
     &                    nlines,nelems,
     &                    imax,jmax,
     &                    lat,lon,
     &                    validTime,
     &                    rdbz,
     &                    maxradars,
     &                    radar_dist_min,
     &                    istatus)
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccc
c
c Routine reads netCDF nowrad (WSI) data using subroutine read_wsi_cdf.
c Nowrad data is then remapped to LAPS domain given lat/lon of domain.
c Routine automatically moves the boundary for any domain with the nowrad
c confines.

      integer maxradars
      integer nsites_present
      integer nsites_absent
      real    present_site_loc_i(maxradars)
      real    present_site_loc_j(maxradars)
      real    absent_site_loc_i(maxradars)
      real    absent_site_loc_j(maxradars)
      real    present_site_lat(maxradars)
      real    present_site_lon(maxradars)
      real    absent_site_lat(maxradars)
      real    absent_site_lon(maxradars)
      real    radar_lat(maxradars)
      real    radar_lon(maxradars)
      real    radar_elev(maxradars)
      real    radar_ri(maxradars)
      real    radar_rj(maxradars)

      real lat(imax,jmax)
      real lon(imax,jmax)
      real height_grid(imax,jmax)
      real rdbz(imax,jmax)
      real radar_dist_min(imax,jmax)
      real ri(imax,jmax)
      real rj(imax,jmax)
      real distmin,dist
      real rdistmax
      real xgrddismx
      real ygrddismx
      real grid_spacing_m
      real dlat,dlon
      real LoV, Latin, La1,La2,Lo1,Lo2
      real centerLon,topLat,dx,dy
      real nw(2),se(2)
      real pi,rad2dg,rii,rjj
      real r_missing_data
      real heigth_mean_grid
      real height_mean
      real height_sum
      real ridis,rjdis

      integer image_to_dbz(0:15)
      integer validTime
      integer istatus
      integer jstatus
      integer status
      integer i_base_value
      integer increment

      character filename*200
      character ctype*3

      integer image(nelems,nlines)
      integer indxofthosein(maxradars)

      common /cegrid/nx,ny,nz,nw,se,rlatc,rlonc

c ----------------------------------------------------------------
c Read data and get navigation info to build lat/lon to i/j table
c ----------------------------------------------------------------
      istatus=1

      nradars_in = 0 ! initialize this

      call read_nowrad_cdf(ctype,filename,nlines,nelems,
     + dlat,dlon,La1,Lo1,La2,Lo2,centerlon,topLat,validTime,
     + dx,dy,Lov, Latin, image, istatus)
c
c this routine converts the raw bytes to integers scaled 0 - 64
c it also determines the wsi grid i/j locations for absent
c or presently reporting radars.

      call cvt_wsi_nowrad(ctype,nelems,nlines,image
     +,nsites_present,present_site_loc_i,present_site_loc_j
     +,nsites_absent,absent_site_loc_i,absent_site_loc_j
     +,maxradars,istatus)
      if(istatus .ne. 1)then
          write(6,*)'Bad status returned from cvt_wsi_nowrad'
          return
      endif

c     call cvt_wsi_nowrad(ctype,nelems,nlines,image,istatus)

      pi=acos(-1.)
      rad2dg=180.0/pi

      if(ctype.eq.'wfo')then

         Lo1=-Lo1
         Lo2=-Lo2
         if(Lov.lt.-180.0)Lov=Lov+360.
         if(Lo1.lt.-180.0)Lo1=Lo1+360.
         if(Lo2.lt.-180.0)Lo2=Lo2+360.

         if(Lov.gt.180.0)Lov=Lov-360.
         if(Lo1.gt.180.0)Lo1=Lo1-360.
         if(Lo2.gt.180.0)Lo2=Lo2-360.

         call gen_rirj_lam(imax,jmax,lat,lon,nelems,nlines
     &        ,La1,Lo1,La2,Lo2,Latin,Lov,dx,dy,ri,rj
     &        ,istatus)
c

      elseif(ctype.eq.'wsi')then
c
c load common block for cyclindrical equidistant
c
           nx=nelems
           ny=nlines
           nz=1
           nw(1)=la1
           nw(2)=lo1
           se(1)=la2
           se(2)=lo2
           rlatc=(topLat - (dlat*(nlines-1)*0.5))*rad2dg
           rlonc=centerLon*rad2dg

c--------------------------------------------------
c compute the wsi grid i/j values for domain lat/lon
c--------------------------------------------------
           call latlon_2_ceij(imax*jmax,lat,lon,ri,rj)

           call check_domain_vrc(imax,jmax,ri,rj,nx,ny)

c--------------------------------------------------
c compute radar site lat/lon locations both present
c and missing for the wsi mosaic.
c--------------------------------------------------
           call ceij_2_latlon(nsites_present,present_site_loc_i
     .,present_site_loc_j,present_site_lat,present_site_lon)

           call ceij_2_latlon(nsites_absent,absent_site_loc_i
     .,absent_site_loc_j,absent_site_lat,absent_site_lon)

      endif

      write(*,*)' NOWRAD data prepared for requested time '
      write(*,*)
      write(*,*)'Valid Time ',validTime

c ---------------------------------------------------
c    Build radar count level to dbz look up
c ---------------------------------------------------
      increment=5
      i_base_value=7
      do i = 1,15
         image_to_dbz(i)=(i-1)*increment+i_base_value
      end do
      image_to_dbz(0)=0
c ---------------------------------------------------
c compute grid ratio
c ---------------------------------------------------
      call get_grid_spacing(grid_spacing_m,istatus)
      r_grid_ratio=sqrt(Dx*Dx + Dy*Dy)/grid_spacing_m
c ---------------------------------------------------
c  Remap NOWRAD-WSI data to LAPS grid.
c ---------------------------------------------------
      Call process_nowrad_z(imax,jmax,
     &                  r_grid_ratio,
     &                  image_to_dbz,
     &                  image,
     &                  ri,
     &                  rj,
     &                  nlines,nelems, ! input array dimensions
     &                  rdbz,
     &                  istatus)

c----------------------------------------------------
c Determine which radars in wsi grid are in this domain
c and compute/save the minimum distance to radar.
c----------------------------------------------------
      call get_r_missing_data(r_missing_data,istatus)
      call  read_static_grid(imax,jmax,'AVG',height_grid,istatus)

      if(ctype.eq.'wsi')then

         height_sum=0
         do j=1,jmax
         do i=1,imax
            height_sum=height_sum+height_grid(i,j)
         enddo
         enddo
         height_mean_grid=height_sum/(imax*jmax)
         nradars_tot=0
         rdistmax=480000.
         xgrddismx=imax+rdistmax/grid_spacing_m
         ygrddismx=jmax+rdistmax/grid_spacing_m
         xgrddismn=-rdistmax/grid_spacing_m
         ygrddismn=-rdistmax/grid_spacing_m
         height_sum=0.0

         print*,'mx num grid points for search x: ',xgrddismx
         print*,'mx num grid points for search y: ',ygrddismx
         print*,'mn num grid points for search x: ',xgrddismn
         print*,'mn num grid points for search y: ',ygrddismn

         do i=1,nsites_present
c
c fudge factor on site lat/lon due to unknown systematic
c offset.
           present_site_lat(i)=present_site_lat(i)-0.07
           present_site_lon(i)=present_site_lon(i)+0.02

           call latlon_to_rlapsgrid(present_site_lat(i),
     &                              present_site_lon(i),
     &                              lat,lon,            !LAPS lat/lon arrays
     &                              imax,jmax,          !LAPS horiz domain
     &                              rii,rjj,            !Output: real i,j, scalars
     &                              jstatus)

           ridis=abs(rii)
           rjdis=abs(rjj)

c          if(ridis.le.xgrddismx .and. rjdis.le.ygrddismx)then

           if( (rii.le.xgrddismx .and. rii.ge.xgrddismn) .and.
     &         (rjj.le.ygrddismx .and. rjj.ge.ygrddismn)  )then

              nradars_tot = nradars_tot + 1
              radar_lat(nradars_tot)=present_site_lat(i)
              radar_lon(nradars_tot)=present_site_lon(i)
              radar_ri(nradars_tot)=rii
              radar_rj(nradars_tot)=rjj
              radar_elev(nradars_tot)=r_missing_data

              if( (rii.gt.0.0  .and.  rii.lt.imax) .and.
     &            (rjj.gt.0.0  .and.  rjj.lt.jmax)  )then

                 call bilinear_interp_extrap(rii,rjj,imax,jmax
     &                    ,height_grid,result,istatus)

                 height_sum=height_sum+result
                 radar_elev(nradars_tot)=result
                 nradars_in=nradars_in+1
                 indxofthosein(nradars_in)=nradars_tot
              endif

           endif

        enddo

        if(nradars_in.gt.0)then
           print*
           print*,'Found ',nradars_in,' in this domain'
           height_mean=height_sum/nradars_in
        else
           print*
           print*,'Interesting! No radars in domain'
           print*,'set mean radar height to mean grid elev'
           print*,' = ',height_mean_grid
           height_mean = height_mean_grid
        endif
        print*
        print*,'Found ',nradars_tot,' to consider in min dist array'
        print*,'Mean height of radars = ',height_mean
        print*

        do i=1,nradars_tot
        do j=1,nradars_in
         ii=indxofthosein(j)
         if(ii.eq.i)then
          print*,'#/ri/rj/lat/lon/elev: ',ii,radar_ri(ii),radar_rj(ii)
     .,radar_lat(ii),radar_lon(ii),radar_elev(ii)
         endif
        enddo
        enddo

        where(radar_elev.eq.r_missing_data)radar_elev=height_mean_grid

        print*,'Determine nearest radar distance array'
        print*

        do j=1,jmax
        do i=1,imax

           distmin=r_missing_data
           do k=1,nradars_tot
              if(radar_elev(k).lt.10000.)then
                call latlon_to_radar(lat(i,j),lon(i,j),height_grid(i,j)
     1,azimuth,slant_range,elev,radar_lat(k),radar_lon(k),radar_elev(k))
                if(slant_range.lt.distmin.and.slant_range.lt.rdistmax)
     &             distmin=slant_range
              endif
           enddo
           radar_dist_min(i,j)=distmin
       enddo
       enddo

      endif

      goto 16
19    print*,'Error in nowrad_2_laps, terminating'
      goto 16
14    print*,'NOWRAD data not found for given time or'
      print*,'Data could be bad. No vrc produced'

16    print*,'Finished in nowrad_2_laps'
      return
      end


text:
module mod_monolis_linalg_com
  use mod_monolis_prm
  use mod_monolis_com
  use mod_monolis_mat
  implicit none

  integer(kint), parameter :: tagSum = 1
  integer(kint), parameter :: tagMax = 2
  integer(kint), parameter :: tagMin = 3

contains

  subroutine monolis_allreduce_R1(val, tag, comm)
    implicit none
    integer(kint), intent(in) :: tag, comm
    integer(kint) :: N, ierr
    real(kdouble) :: val, temp

#ifdef WITH_MPI
    temp = 0.0d0
    N = 1
    if(tag == tagSum)then
      call MPI_allreduce(val, temp, N, MPI_DOUBLE_PRECISION, MPI_SUM, comm, ierr)
    elseif(tag == tagMax)then
      call MPI_allreduce(val, temp, N, MPI_DOUBLE_PRECISION, MPI_MAX, comm, ierr)
    elseif(tag == tagMin)then
      call MPI_allreduce(val, temp, N, MPI_DOUBLE_PRECISION, MPI_MIN, comm, ierr)
    endif
    val = temp
#endif
  end subroutine monolis_allreduce_R1

  subroutine monolis_allreduce_R(N, val, tag, comm)
    implicit none
    integer(kint), intent(in) :: tag, comm
    integer(kint) :: N, ierr
    real(kdouble) :: val(N), temp(N)

#ifdef WITH_MPI
    temp = 0.0d0
    if(tag == tagSum)then
      call MPI_allreduce(val, temp, N, MPI_DOUBLE_PRECISION, MPI_SUM, comm, ierr)
    elseif(tag == tagMax)then
      call MPI_allreduce(val, temp, N, MPI_DOUBLE_PRECISION, MPI_MAX, comm, ierr)
    elseif(tag == tagMin)then
      call MPI_allreduce(val, temp, N, MPI_DOUBLE_PRECISION, MPI_MIN, comm, ierr)
    endif
    val = temp
#endif
  end subroutine monolis_allreduce_R

  subroutine monolis_allreduce_I1(val, tag, comm)
    implicit none
    integer(kint), intent(in) :: tag, comm
    integer(kint)  :: N, ierr, val, temp

#ifdef WITH_MPI
    temp = 0
    N = 1
    if(tag == tagSum)then
      call MPI_allreduce(val, temp, N, MPI_INTEGER, MPI_SUM, comm, ierr)
    elseif(tag == tagMax)then
      call MPI_allreduce(val, temp, N, MPI_INTEGER, MPI_MAX, comm, ierr)
    elseif(tag == tagMin)then
      call MPI_allreduce(val, temp, N, MPI_INTEGER, MPI_MIN, comm, ierr)
    endif
    val = temp
#endif
  end subroutine monolis_allreduce_I1

  subroutine monolis_allreduce_I(N, val, tag, comm)
    implicit none
    integer(kint), intent(in) :: tag, comm
    integer(kint)  :: N, ierr, val(N), temp(N)

#ifdef WITH_MPI
    temp = 0
    if(tag == tagSum)then
      call MPI_allreduce(val, temp, N, MPI_INTEGER, MPI_SUM, comm, ierr)
    elseif(tag == tagMax)then
      call MPI_allreduce(val, temp, N, MPI_INTEGER, MPI_MAX, comm, ierr)
    elseif(tag == tagMin)then
      call MPI_allreduce(val, temp, N, MPI_INTEGER, MPI_MIN, comm, ierr)
    endif
    val = temp
#endif
  end subroutine monolis_allreduce_I

  subroutine monolis_SendRecv_pre_R(send_n_neib, send_neib_pe, recv_n_neib, recv_neib_pe, &
    & send_index, send_item, recv_index, recv_item, &
    & ws, wr, val, ndof, comm)
    implicit none
    integer(kint) :: send_n_neib, recv_n_neib
    integer(kint) :: iS, in, j, k, ierr
    integer(kint) :: i, ndof, comm
    integer(kint), pointer :: send_neib_pe(:)
    integer(kint), pointer :: send_index(:)
    integer(kint), pointer :: send_item (:)
    integer(kint), pointer :: recv_neib_pe(:)
    integer(kint), pointer :: recv_index(:)
    integer(kint), pointer :: recv_item (:)
    integer(kint) :: sta1(monolis_status_size, send_n_neib)
    integer(kint) :: sta2(monolis_status_size, recv_n_neib)
    integer(kint) :: req1(send_n_neib)
    integer(kint) :: req2(recv_n_neib)
    real(kdouble) :: val(:), ws(:), wr(:)

#ifdef WITH_MPI
    do i = 1, send_n_neib
      iS = send_index(i-1)
      in = send_index(i  ) - iS
      if(in == 0) cycle
      do j = iS+1, iS+in
        do k = 1, ndof
          ws(ndof*(j-1)+k) = val(ndof*(send_item(j)-1)+k)
        enddo
      enddo
      call MPI_Isend(ws(ndof*iS+1), ndof*in, MPI_DOUBLE_PRECISION, send_neib_pe(i), 0, comm, req1(i), ierr)
    enddo

    do i = 1, recv_n_neib
      iS = recv_index(i-1)
      in = recv_index(i  ) - iS
      if(in == 0) cycle
      call MPI_Irecv(wr(ndof*iS+1), ndof*in, MPI_DOUBLE_PRECISION, recv_neib_pe(i), 0, comm, req2(i), ierr)
    enddo

    call MPI_waitall(recv_n_neib, req2, sta2, ierr)

    do i = 1, recv_n_neib
      iS = recv_index(i-1)
      in = recv_index(i  ) - iS
      do j = iS+1, iS+in
        do k = 1, ndof
          val(ndof*(recv_item(j)-1)+k) = wr(ndof*(j-1)+k)
        enddo
      enddo
    enddo

    call MPI_waitall(send_n_neib, req1, sta1, ierr)
#endif
  end subroutine monolis_SendRecv_pre_R

  subroutine monolis_SendRecv_post_R(n_neib, neib_pe, send_index, send_item, recv_index, recv_item, &
  & ws, wr, val, ndof, comm)
    implicit none
    integer(kint) :: n_neib
    integer(kint) :: iS, in, j, k, ierr
    integer(kint) :: i, ndof, comm
    integer(kint), pointer :: neib_pe(:)
    integer(kint), pointer :: send_index(:)
    integer(kint), pointer :: send_item (:)
    integer(kint), pointer :: recv_index(:)
    integer(kint), pointer :: recv_item (:)
    integer(kint) :: sta1(monolis_status_size, n_neib)
    integer(kint) :: sta2(monolis_status_size, n_neib)
    integer(kint) :: req1(n_neib)
    integer(kint) :: req2(n_neib)
    real(kdouble) :: val(:), ws(:), wr(:)

#ifdef WITH_MPI

#endif
  end subroutine monolis_SendRecv_post_R

  subroutine monolis_SendRecv_pre_I(send_n_neib, send_neib_pe, recv_n_neib, recv_neib_pe, &
    & send_index, send_item, recv_index, recv_item, &
    & ws, wr, val, ndof, comm)
    implicit none
    integer(kint) :: send_n_neib, recv_n_neib
    integer(kint) :: iS, in, j, k, ierr
    integer(kint) :: i, ndof, comm
    integer(kint), pointer :: send_neib_pe(:)
    integer(kint), pointer :: send_index(:)
    integer(kint), pointer :: send_item (:)
    integer(kint), pointer :: recv_neib_pe(:)
    integer(kint), pointer :: recv_index(:)
    integer(kint), pointer :: recv_item (:)
    integer(kint) :: sta1(monolis_status_size, send_n_neib)
    integer(kint) :: sta2(monolis_status_size, recv_n_neib)
    integer(kint) :: req1(send_n_neib)
    integer(kint) :: req2(recv_n_neib)
    integer(kint) :: val(:), ws(:), wr(:)

#ifdef WITH_MPI
    do i = 1, send_n_neib
      iS = send_index(i-1)
      in = send_index(i  ) - iS
      if(in == 0) cycle
      do j = iS+1, iS+in
        do k = 1, ndof
          ws(ndof*(j-1)+k) = val(ndof*(send_item(j)-1)+k)
        enddo
      enddo
      call MPI_Isend(ws(ndof*iS+1), ndof*in, MPI_INTEGER, send_neib_pe(i), 0, comm, req1(i), ierr)
    enddo

    do i = 1, recv_n_neib
      iS = recv_index(i-1)
      in = recv_index(i  ) - iS
      if(in == 0) cycle
      call MPI_Irecv(wr(ndof*iS+1), ndof*in, MPI_INTEGER, recv_neib_pe(i), 0, comm, req2(i), ierr)
    enddo

    call MPI_waitall(recv_n_neib, req2, sta2, ierr)

    do i = 1, recv_n_neib
      iS = recv_index(i-1)
      in = recv_index(i  ) - iS
      do j = iS+1, iS+in
        do k = 1, ndof
          val(ndof*(recv_item(j)-1)+k) = wr(ndof*(j-1)+k)
        enddo
      enddo
    enddo

    call MPI_waitall(send_n_neib, req1, sta1, ierr)
#endif
  end subroutine monolis_SendRecv_pre_I

  subroutine monolis_allgather_I1(sval, rbuf, comm)
    implicit none
    integer(kint) :: sval    !send buffer
    integer(kint) :: rbuf(:) !receive buffer
    integer(kint) :: comm
#ifdef WITH_MPI
    integer(kint) :: ierr
    call MPI_allgather(sval, 1, MPI_INTEGER, rbuf, 1, MPI_INTEGER, comm, ierr )
#else
    rbuf(1) = sval
#endif
  end subroutine monolis_allgather_I1

  subroutine monolis_gatherv_R(sbuf, sc, &
      &     rbuf, rcs, disp, root, comm)
    implicit none
    integer(kint) :: sc      !send counts
    integer(kint) :: rcs(:)  !receive counts
    integer(kint) :: disp(:) !displacement
    integer(kint) :: root
    integer(kint) :: comm
    real(kdouble) :: sbuf(:) !send buffer
    real(kdouble) :: rbuf(:) !receive buffer
#ifdef WITH_MPI
    integer(kint) :: ierr
    call MPI_gatherv( sbuf, sc, MPI_REAL8, &
      &     rbuf, rcs, disp, MPI_REAL8, root, comm, ierr )
#else
    rbuf(1:sc) = sbuf(1:sc)
#endif
  end subroutine monolis_gatherv_R

  subroutine monolis_scatterv_R(sbuf, scs, disp, &
      &     rbuf, rc, root, comm)
    implicit none
    integer(kint) :: scs(:)  !send counts
    integer(kint) :: disp(:) !displacement
    integer(kint) :: rc      !receive counts
    integer(kint) :: root
    integer(kint) :: comm
    real(kdouble) :: rbuf(:) !receive buffer
    real(kdouble) :: sbuf(:) !send buffer
#ifdef WITH_MPI
    integer(kint) :: ierr
    call MPI_scatterv( sbuf, scs, disp, MPI_REAL8, &
      &     rbuf, rc, MPI_REAL8, root, comm, ierr )
#else
    rbuf(1:rc) = sbuf(1:rc)
#endif
  end subroutine monolis_scatterv_R

end module mod_monolis_linalg_com

max_stars_repo_name: nqomorita/monolis

max_issues_repo_name: nqomorita/monolis

max_forks_repo_name: nqomorita/monolis