text:
The Office of External Affairs is a physical location on the third floor of the Memorial Union, housing three ASUCD units which deal with advocacy:

  ASUCD City/County Affairs Office
  ASUCD Lobby Corps
  ASUCD University Affairs Office

It was formerly the name of a consolidated unit that oversaw all ASUCD Advocacy operations, though it was unwieldy in ASUCDs flat organizational structure, and was disbanded.  Details are below.

 ASUCD Unit

The Office of External Affairs was a part of the ASUCD from 2004 to 2006.  The OEA served to educate and advocate for undergraduate students on city, county, state and national issues. The OEA involved students in this advocacy through voter registration and education, grassroots organizing, lobbying elected officials, educational forums, and other out reach projects.

The unit was part of the ASUCD President ASUCD Presidents cabinet and worked with him on related special projects and issues. The OEA also worked with the University of California Student Association (UCSA) and the United States Student Association (USSA) on UC system wide, state, and federal legislative issues related to higher education.

In FY20062007 the office was broken up into three units; ASUCD Lobby Corps, ASUCD City/County Affairs Office, and ASUCD State and National Affairs Office (SNAO). However, the SNAO became the University Affairs Office in 2007 in reaction to an advocacy shift cause by ASUCDs pull out from UCSA.  For more information, see the History section on this entry.

 Offices within External Affairs

  ASUCD City/County Affairs Office
  ASUCD Lobby Corps
  ASUCD State and National Affairs Office

 Staff
  20062007 (Constituent units)
   Michelle Shaffie: Director of ASUCD State and National Affairs Office
   Matt Shannon: Assistant Director of ASUCD State and National Affairs
   Anneka Euser: Campus Organizing

   Derick Lennox: External Director ASUCD Lobby Corps
   Danielle Vaticano: Internal Director ASUCD Lobby Corps
   Erica Rozetti: Field/Recruitment Organizer

   Jarvis Burston: Director of ASUCD City and County Affairs Office

  20052006
   James Schwab: Director of the Office of External Affairs 20052006
   Liz Burrell: State and National Affairs Director
   Marwa Rifahie: Assistant Director, Legislative Affairs
   Users/AliRamezanzadeh Ali Ramezanzadeh: Campus Outreach Director
   Greg Russell: ASUCD City/County Affairs Office City and County Affairs Director
   Laura Henderson: Assistant Director, Voter Registration
   Charles Palmares: Assistant Director, Education and Forums
   Users/WencheMolenaar Wenche Molenaar: ASUCD StudentPolice Relations Committee Chair
   Users/MichaelCaporusso Mike Caporusso ASUCD Lobby Corps Director
  20042005
   Brian McInnis: Director of the Office of External Affairs

 History

The ASUCD City/County Lobby evolved from the City Student Lobby.  The City Student Lobby was established during the 19731974 academic year. This was in response to the need for effective student representation within the city of Davis and the county of Yolo.  Students, as a distinct sector of the population, have special needs and interests which the City/County Lobby (now the City/Council Affairs Office) can meet.  Because the quality of students lives and education are affected by community decisions, students require a means of influencing the local decisionmaking process.

In 1993, due to the ruling by the California Supreme Court on Smith v. Regents, structural and programmatic changes were necessary for all of ASUCDs Lobby Units.  The decision limited the use of mandatory student fees for advocating ideological or political causes outside the U.C. system.  This ruling in effect required the cessation of all advocacy functions of this office, since the office was supported by mandatory ASUCD fees.  Thus, the ASUCD City/County Lobby was renamed ASUCD City/County Affairs Office to reflect the instituted changes and the offices new role as a purely educational and outreach natured unit.

The need for a bridge between students and their local elected officials is fulfilled by the City/County Affairs Office.  The Office analyzed the decisions and actions of the Davis City Council, the Yolo County Board of Supervisors, and the various commissions within the city and county.  The Office researched community issues, and analyzed their potential impact on Davis students.  The Office publicized pertinent information among the student population.  Students educated in this way are better able to make informed decisions about the community issues which affect them.  The Office then solicited student input to identify student positions, enabling the City/County Affairs Office to responsibly represent the student interest.

The Office dealt with housing issues as well as supplying apartment complex information.  At one time, a specialized branch of the Housing Department in the City/County Lobby was Coop Housing.  Previously, the City/County Affairs Office sponsored the ASUCD Housing Day which is held annually in order to provide needed information to students and apartment managers.

The Office had a Community Relations Department which had two major functions: to act as Community City/County Liaison between city and the county organizations and to conduct Community Projects.  The Lobby also had an Environmental Issues Coordinator.  Next to Project Recycle, this is the only department which specifically dealt with environmental issues within ASUCD.  This unit worked with EPPC to develop community wide environmental projects and influence environmental policy made locally.

As a distinct member of the ASUCD cabinet, the City/County Affairs Office served the students of Davis as an information source for the students within the Associated Students of the University of California at Davis.

In 1989, the ASUCD Executive Council created the ASUCD National Lobby in order to provide UCD Students with a voice on the federal government level. Similar to the ASUCD State Lobby, this advocacy unit deals with federal issues and advocates the viewpoint of UC Davis students to the US Congress, Department of Education, and other federal agencies dealing with higher education.

In 1993, again due to the ruling by the California Supreme Court on Smith v. Regents, structural and programmatic changes were necessary for all of ASUCDs Lobby units. The ASUCD National Lobby was renamed ASUCD National Affairs Office to reflect the instituted changes and the offices new role as a purely research, educational, and outreach natured unit. The ASUCD National Lobby and State Lobby coordinated extensively with both the UC Student Association and the United States Student Association to help secure the state and federal policy dreams students hold for higher education.

In spring of 2003, the ASUCD Lobby Corps was created to teach UCD undergraduates how to lobby and better represent the student body in state and local government.  UCD alumni and administrators come as guest speakers to teach students how to support then retitled ASUCD State and National Affairs Office as well as City/County Affairs Office at the city, state, and national level.

In 2004, the ASUCD Executive Office connected the three main advocative efforts of the association by creating the ASUCD Office of External Affairs. This new office incorporates advocacy in the city of Davis, the county of Yolo, the state of California, as well as in the Federal Government. Through research, analysis, and various lobbying tactics the External Affairs office will strive to represent UC Davis students and their concerns to the best of their abilities. A central role of the ASUCD Office of External Affairs will be to incorporate a broad base of UC Davis students not only in comprehension of these external issues, but also to empower students through lobby and direct action trainings and experience.

The ASUCD External Affairs Office is meant as a resource for UC Students and to act as that resource through the Executive Branch of ASUCD. The ASUCD Office of External Affairs is meant to be primarily directed by the Executive Branch and from APPROVED ASUCD Senate Resolution or Bills pertaining to external issues facing students at UC Davis.

In 2006, the Office of External Affairs was broken back up into its pre2004 units, City/County Affairs and State and National Affairs as well as Lobby Corps branching off to become its own unit. The position of Director of External Affairs was moved into the Executive Office as a cabinet position.  The following year, the Director of External Affairs was no longer needed.

 Gossipy History

Legend tells us that the OEA was only created because Kalen Gallagher had promised both Liz Burrell and Brian McInnis the job of Director of State and National Affairs.  The solution?  Create a position for McInnis above SNAO Director, and make everyone happy.  Not that it actually did.

20080930 07:30:00 nbsp Legend tells us?  Is that how it happened?  Interesting gossipy history... Users/EBurrell

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! RUN: %python %S/test_errors.py %s %flang_fc1
! Check errors on MAX/MIN with keywords, a weird case in Fortran
real :: x = 0.0 ! prevent folding
!ERROR: Argument keyword 'a1=' was repeated in call to 'max'
print *, max(a1=x,a1=1)
!ERROR: Keyword argument 'a1=' has already been specified positionally (#1) in this procedure reference
print *, max(x,a1=1)
print *, max(a1=x,a2=0,a4=0) ! ok
print *, max(x,0,a99=0) ! ok
!ERROR: Argument keyword 'a06=' is not known in call to 'max'
print *, max(a1=x,a2=0,a06=0)
end


text:
      subroutine wshrun
c
c     + + + PURPOSE + + +
c
c     SR WSHRUN routes the runoff from impoundments and
c     channels through the watershed.
c
c     Called from: SR WSHDRV
c     Author(s): Ascough II
c     Reference in User Guide:
c
c     Version:
c     Date recoded:
c     Recoded by:
c
c     + + + KEYWORDS + + +
c
c     + + + PARAMETERS + + +
c
      include 'pmxcsg.inc'
      include 'pmxelm.inc'
      include 'pmxhil.inc'
      include 'pmximp.inc'
      include 'pmxpln.inc'
      include 'pmxprt.inc'
c
c     + + + ARGUMENT DECLARATIONS + + +
c
c     + + + ARGUMENT DEFINITIONS + + +
c
c     + + + COMMON BLOCKS + + +
c
      include 'cdiss11.inc'
      include 'cirspri.inc'
      include 'coutchn.inc'
      include 'cstore1.inc'
      include 'cstruct.inc'

      include 'chydrol.inc'
      include 'cchpek.inc'
      include 'cstruc.inc'
c
c     + + + LOCAL VARIABLES + + +
c
c     + + + LOCAL DEFINITIONS + + +
c
c     + + + SAVES + + +
c
c     + + + SUBROUTINES CALLED + + +
c
c     wshpek
c
c     + + + DATA INITIALIZATIONS + + +
c
c     + + + END SPECIFICATIONS + + +
c
c
c     calculate the average runoff on the watershed at the outlet
c     of the current channel element
c
      rofave(ielmt) = runvol(ielmt) / wsarea(ielmt)
c
c     save monthly, annual, and total runoff from each channel
c     element
c
      trunm(ielmt) = trunm(ielmt) + runvol(ielmt)
      truny(ielmt) = truny(ielmt) + runvol(ielmt)
      trunt(ielmt) = trunt(ielmt) + runvol(ielmt)
cd    Added by S. Dun, Dec 20, 2007
c     for more output of annual summary of the chnnel
      tronvy(ielmt) = tronvy(ielmt) + rvolon(ielmt)
      tronvt(ielmt) = tronvt(ielmt) + rvolon(ielmt)
      sbrvty(ielmt) = sbrvty(ielmt) + sbrunf(ielmt-nhill)
     1               * charea(ielmt-nhill)
      sbrvtt(ielmt) = sbrvtt(ielmt) + sbrunf(ielmt-nhill)
     1               * charea(ielmt-nhill)
c     End adding
c
c     the duration for the channel element is the maximum
c     duration of the watershed elements that may feed it
c
cd    Modified by S. Dun 1/08/2004
cd    This duration value has given in WSHCQI already.
cd    Only the dur for channel could be changed after WSHIRS.for
cd
cd      watdur(ielmt) = amax1(watdur(nhleft(ielmt)),
cd     1    watdur(nhrght(ielmt)),watdur(nhtop(ielmt)),
cd     1    watdur(ncleft(ielmt)),watdur(ncrght(ielmt)),
cd     1    watdur(nctop(ielmt)),watdur(nileft(ielmt)),
cd     1    watdur(nirght(ielmt)),watdur(nitop(ielmt)),dur,irdur)
      watdur(ielmt) = amax1(watdur(ielmt),dur)
cd    End modifying
c
c
c     calculate peak runoff using modified Rational method
c     (ipeak = 1) or modified CREAMS method (ipeak = 2)
c
      call wshpek
c
      return
      end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text:
    module handles
    use CAMB
    use Precision
    use ModelParams
    use Transfer
    use iso_c_binding
    implicit none

    Type c_MatterTransferData
        integer   ::  num_q_trans   !    number of steps in k for transfer calculation
        type(c_ptr) :: q_trans, sigma_8, sigma2_vdelta_8, TransferData
        integer :: sigma_8_size(2)
        integer :: sigma2_vdelta_8_size(2)
        integer :: TransferData_size(3)
    end Type c_MatterTransferData

    Type c_ClTransferData
        integer :: NumSources
        integer :: q_size
        type(c_ptr) :: q
        integer delta_size(3)
        type(c_ptr) Delta_p_l_k
        integer l_size
        type(c_ptr) ls
        !skip limber for now...
    end Type c_ClTransferData

    contains

    !SPECIAL BRIDGE ROUTINES FOR PYTHONset

    subroutine CAMBdata_new(handle)
    type(c_ptr), intent(out) :: handle
    Type (CAMBdata), pointer :: pCAMBdata

    allocate(pCAMBdata)
    call CAMB_InitCAMBdata(pCAMBdata)
    call CAMB_setDefParams(pCAMBdata%Params)
    handle = c_loc(pCAMBdata)

    end subroutine CAMBdata_new

    subroutine CAMBdata_free(cptr)
    type(c_ptr)  :: cptr
    Type (CAMBdata), pointer :: pCAMBdata

    call c_f_pointer(cptr, pCAMBdata)
    call CAMB_FreeCAMBdata(pCAMBdata)
    deallocate(pCAMBdata)

    end subroutine CAMBdata_free

    subroutine CAMBdata_setParams(data, Params)
    Type (CAMBdata), target :: data
    type(CAMBparams) :: Params

    data%Params = Params

    end subroutine CAMBdata_setParams

    subroutine CAMBdata_getParams(data, handle)
    Type (CAMBdata), target :: data
    type(c_ptr), intent(out)  ::  handle

    handle = c_loc(data%Params)

    end subroutine CAMBdata_getParams

    function CAMBdata_GetTransfers(data, Params, onlytransfer) result(error)
    Type (CAMBdata):: data
    type(CAMBparams) :: Params, P
    logical(kind=c_bool)  :: onlytransfer
    integer :: error

    P = Params
    if (P%DoLensing .and. (P%NonLinear == NonLinear_lens .or. P%NonLinear == NonLinear_both)) then
        P%WantTransfer = .true.
        call Transfer_SetForNonlinearLensing(P%Transfer)
    end if
    call Transfer_SortAndIndexRedshifts(P%Transfer)
    error = 0

    P%OnlyTransfers = onlytransfer
    call CAMB_GetTransfers(P, data, error)

    end function CAMBdata_GetTransfers

    subroutine CAMBdata_SetParamsForBackground(data, P)
    Type (CAMBdata):: data
    type(CAMBparams) :: P

    global_error_flag = 0
    data%Params = P
    call CAMBParams_Set(data%Params)
    end subroutine CAMBdata_SetParamsForBackground

    function CAMBdata_CalcBackgroundTheory(data, P) result(error)
    use cambmain, only: initvars
    Type (CAMBdata):: data
    type(CAMBparams) :: P
    integer error

    global_error_flag = 0
    data%Params = P
    call CAMBParams_Set(data%Params)
    call InitVars !calculate thermal history, e.g. z_drag etc.
    error=global_error_flag

    end function CAMBdata_CalcBackgroundTheory

    subroutine CAMBdata_MatterTransferData(data, cData)
    Type(CAMBdata), target :: data
    Type(c_MatterTransferData) :: cData

    cData%num_q_trans = data%MTrans%num_q_trans
    cData%q_trans = c_loc(data%MTrans%q_trans)
    cData%sigma_8 = c_loc(data%MTrans%sigma_8)
    cData%sigma2_vdelta_8 = c_loc(data%MTrans%sigma2_vdelta_8)
    cData%TransferData = c_loc(data%MTrans%TransferData)
    cData%q_trans = c_loc(data%MTrans%q_trans)
    cData%sigma_8_size = shape(data%MTrans%sigma_8)
    cData%sigma2_vdelta_8_size = shape(data%MTrans%sigma2_vdelta_8)
    cData%TransferData_size = shape(data%MTrans%TransferData)

    end subroutine CAMBdata_MatterTransferData

    subroutine CAMBdata_ClTransferData(data, cData, i)
    Type(CAMBdata), target :: data
    Type(c_ClTransferData) :: cData
    integer, intent(in) :: i

    if (i==0) then
        call Convert_ClTransferData(data%ClTransScal, cData)
    else if (i==1) then
        call Convert_ClTransferData(data%ClTransVec, cData)
    else if (i==2) then
        call Convert_ClTransferData(data%ClTransTens, cData)
    else
        error stop 'Unknown ClTransferData index'
    end if

    end subroutine CAMBdata_ClTransferData

    subroutine Convert_ClTransferData(data, cData)
    Type(ClTransferData), target :: data
    Type(c_ClTransferData) :: cData

    cData%NumSources = data%NumSources
    if (associated(Data%q%points)) then
        cData%q_size = size(data%q%points)
        cData%q = c_loc(data%q%points)
    else
        cData%q_size = 0
    end if
    if (associated(data%Delta_p_l_k)) then
        cData%delta_size = shape(Data%Delta_p_l_k)
        cData%delta_p_l_k = c_loc(Data%Delta_p_l_k)
    else
        cData%delta_size = 0
    end if
    cdata%l_size = Data%ls%l0
    cdata%ls = c_loc(Data%ls%l)

    end subroutine Convert_ClTransferData

    subroutine CAMBdata_GetLinearMatterPower(data, PK, var1, var2, hubble_units)
    Type(CAMBdata) :: data
    real(dl) :: PK(data%MTrans%num_q_trans,data%Params%Transfer%PK_num_redshifts)
    integer, intent(in) :: var1, var2
    logical :: hubble_units

    call Transfer_GetUnsplinedPower(data%MTrans, PK, var1, var2, hubble_units)

    end subroutine CAMBdata_GetLinearMatterPower

    subroutine CAMBdata_GetNonLinearMatterPower(data, PK, var1, var2, hubble_units)
    Type(CAMBdata) :: data
    real(dl) :: PK(data%MTrans%num_q_trans,data%Params%Transfer%PK_num_redshifts)
    integer, intent(in) :: var1, var2
    logical :: hubble_units

    call Transfer_GetUnsplinedNonlinearPower(data%MTrans, PK, var1, var2, hubble_units)

    end subroutine CAMBdata_GetNonLinearMatterPower

    subroutine CAMBdata_GetMatterPower(data, outpower, minkh, dlnkh, npoints, var1, var2)
    Type(CAMBdata) :: data
    real(dl), intent(out) :: outpower(npoints,data%Params%Transfer%PK_num_redshifts)
    real(dl), intent(in) :: minkh, dlnkh
    integer, intent(in) :: npoints, var1, var2
    integer i

    do i=1,data%Params%Transfer%PK_num_redshifts
        call Transfer_GetMatterPowerD(data%MTrans,outpower(:,i), data%Params%Transfer%PK_num_redshifts-i+1, &
            & 1, minkh, dlnkh, npoints, var1, var2)
    end do

    end subroutine CAMBdata_GetMatterPower

    subroutine CAMB_setinitialpower(Params, P)
    type(CAMBparams) :: Params
    type(InitialPowerParams) :: P

    Params%InitPower = P

    end subroutine CAMB_setinitialpower

    subroutine CAMB_SetTotCls(lmax, tot_scalar_Cls, in)
    integer, intent(IN) :: lmax, in
    real(dl), intent(OUT) :: tot_scalar_cls(4, 0:lmax)
    integer l

    tot_scalar_cls = 0
    do l=lmin, lmax
        if (CP%WantScalars .and. l<= CP%Max_l) then
            if (CP%DoLensing) then
                if (l<=lmax_lensed) tot_scalar_cls(1:4,l) = Cl_lensed(l, in, CT_Temp:CT_Cross)
            else
                tot_scalar_cls(1:2,l) = Cl_scalar(l, in,  C_Temp:C_E)
                tot_scalar_cls(4,l) = Cl_scalar(l, in,  C_Cross)
            endif
        end if
        if (CP%WantTensors .and. l <= CP%Max_l_tensor) then
            tot_scalar_cls(1:4,l) = tot_scalar_cls(1:4,l) + Cl_tensor(l, in,  CT_Temp:CT_Cross)
        end if
    end do

    end subroutine CAMB_SetTotCls

    subroutine CAMB_SetUnlensedCls(lmax, unlensed_cls, in)
    integer, intent(IN) :: lmax, in
    real(dl), intent(OUT) :: unlensed_cls(4,0:lmax)
    integer l

    unlensed_cls = 0
    do l=lmin, lmax
        if (CP%WantScalars .and. l<= CP%Max_l) then
            unlensed_cls(1:2,l) = Cl_scalar(l, in,  C_Temp:C_E)
            unlensed_cls(4,l) = Cl_scalar(l, in,  C_Cross)
        end if
        if (CP%WantTensors .and. l <= CP%Max_l_tensor) then
            unlensed_cls(1:4,l) = unlensed_cls(1:4,l) + Cl_tensor(l, in,  CT_Temp:CT_Cross)
        end if
    end do

    end subroutine CAMB_SetUnlensedCls

    subroutine CAMB_SetLensPotentialCls(lmax, cls, in)
    use constants
    integer, intent(IN) :: lmax, in
    real(dl), intent(OUT) :: cls(3, 0:lmax) !phi-phi, phi-T, phi-E
    integer l

    cls = 0
    if (CP%WantScalars .and. CP%DoLensing) then
        do l=lmin, min(lmax,CP%Max_l)
            cls(1,l) = Cl_scalar(l,in,C_Phi) * (real(l+1)/l)**2/const_twopi
            cls(2:3,l) = Cl_scalar(l,in,C_PhiTemp:C_PhiE) * ((real(l+1)/l)**1.5/const_twopi)
        end do
    end if

    end subroutine CAMB_SetLensPotentialCls

    subroutine CAMB_SetUnlensedScalCls(lmax, scalar_Cls, in)
    integer, intent(IN) :: lmax, in
    real(dl), intent(OUT) :: scalar_Cls(4, 0:lmax)
    integer lmx

    scalar_Cls = 0
    if (CP%WantScalars) then
        lmx = min(CP%Max_l, lmax)
        scalar_Cls(1:2,lmin:lmx) = transpose(Cl_Scalar(lmin:lmx, in,C_Temp:C_E))
        scalar_Cls(4,lmin:lmx) = Cl_Scalar(lmin:lmx, in,C_Cross)
    end if

    end subroutine CAMB_SetUnlensedScalCls

    subroutine CAMB_SetlensedScalCls(lmax, lensed_Cls, in)
    integer, intent(IN) :: lmax, in
    real(dl), intent(OUT) :: lensed_Cls(4, 0:lmax)
    integer lmx

    lensed_Cls = 0
    if (CP%WantScalars .and. CP%DoLensing) then
        lmx = min(lmax,lmax_lensed)
        lensed_Cls(1:4,lmin:lmx) = transpose(Cl_lensed(lmin:lmx, in,CT_Temp:CT_Cross))
    end if

    end subroutine CAMB_SetlensedScalCls

    subroutine CAMB_SetTensorCls(lmax, tensor_Cls, in)
    integer, intent(IN) :: lmax, in
    real(dl), intent(OUT) :: tensor_Cls(4, 0:lmax)
    integer lmx

    tensor_Cls = 0
    if (CP%WantTensors) then
        lmx = min(lmax,CP%Max_l_tensor)
        tensor_Cls(1:3,lmin:lmx) = transpose(Cl_Tensor(lmin:lmx, in, CT_Temp:CT_Cross))
    end if

    end subroutine CAMB_SetTensorCls

    subroutine CAMB_SetUnlensedScalarArray(lmax, ScalarArray, in, n)
    integer, intent(IN) :: lmax, in, n
    real(dl), intent(OUT) :: ScalarArray(n, n, 0:lmax)
    integer l

    ScalarArray = 0
    if (CP%WantScalars) then
        do l=lmin, min(lmax,CP%Max_l)
            ScalarArray(1:n,1:n,l) = Cl_scalar_array(l, in, 1:n,1:n)
        end do
    end if

    end subroutine CAMB_SetUnlensedScalarArray

    subroutine CAMB_SetBackgroundOutputs_z(redshifts,n)
    integer, intent(in) :: n
    real(dl), intent(in) :: redshifts(n)

    if (associated(BackgroundOutputs%z_outputs)) deallocate(BackgroundOutputs%z_outputs)
    if (n>0) then
        allocate(BackgroundOutputs%z_outputs(n))
        BackgroundOutputs%z_outputs = redshifts
    else
        nullify(BackgroundOutputs%z_outputs)
    end if

    end subroutine CAMB_SetBackgroundOutputs_z

    function CAMB_GetNumBackgroundOutputs()
    integer CAMB_GetNumBackgroundOutputs

    if (.not. associated(BackgroundOutputs%z_outputs)) then
        CAMB_GetNumBackgroundOutputs = 0
    else
        CAMB_GetNumBackgroundOutputs = size(BackgroundOutputs%z_outputs)
    end if

    end function CAMB_GetNumBackgroundOutputs

    subroutine CAMB_GetBackgroundOutputs(outputs, n)
    use constants
    integer, intent(in) :: n
    real(dl), intent(out) :: outputs(4,n)
    integer i

    if (associated(BackgroundOutputs%z_outputs)) then
        do i=1, size(BackgroundOutputs%z_outputs)
            outputs(1,i) = BackgroundOutputs%rs_by_D_v(i)
            outputs(2,i) = BackgroundOutputs%H(i)*c/1e3_dl
            outputs(3,i) = BackgroundOutputs%DA(i)
            outputs(4,i) = (1+BackgroundOutputs%z_outputs(i))* &
                BackgroundOutputs%DA(i) * BackgroundOutputs%H(i) !F_AP parameter
        end do
    end if

    end subroutine CAMB_GetBackgroundOutputs

    subroutine set_cls_template(cls_template)
    character(len=*), intent(in) :: cls_template

    highL_unlensed_cl_template = trim(cls_template)

    end subroutine set_cls_template

    function CAMB_PrimordialPower(Params, k, powers, n,  i) result(err)
    use constants
    type(CAMBparams) :: Params
    integer, intent(in) :: i,n
    real(dl), intent(in) :: k(n)
    real(dl), intent(out) :: powers(n)
    real(dl) curv
    integer err,ix

    global_error_flag = 0
    curv =-Params%omegak/((c/1000)/Params%h0)**2
    call InitializePowers(Params%InitPower,curv)
    if (global_error_flag==0) then
        do ix =1, n
            if (i==0) then
                powers(ix) = ScalarPower(k(ix),1)
            elseif (i==2) then
                powers(ix) = TensorPower(k(ix),1)
            else
                error stop 'Unknown power type index'
            end if
            if (global_error_flag /= 0) exit
        end do
    end if
    err= global_error_flag

    end function CAMB_PrimordialPower

    subroutine GetOutputEvolutionFork(EV, times, outputs)
    use Transfer
    use CAMBmain
    implicit none
    type(EvolutionVars) EV
    real(dl), intent(in) :: times(:)
    real(dl), intent(out) :: outputs(:,:,:)
    real(dl) tau,tol1,tauend, taustart
    integer j,ind,itf
    real(dl) c(24),w(EV%nvar,9), y(EV%nvar)
    real(dl) yprime(EV%nvar), ddelta, delta, adotoa,dtauda, growth, x_e, a, x_e_recomb
    external dtauda
    real, target :: Arr(Transfer_max)

    w=0
    y=0
    taustart = GetTauStart(min(500._dl,EV%q))
    call initial(EV,y, taustart)

    tau=taustart
    ind=1
    tol1=tol/exp(AccuracyBoost-1)
    do j=1,size(times)
        tauend = times(j)
        if (tauend<taustart) cycle

        call GaugeInterface_EvolveScal(EV,tau,y,tauend,tol1,ind,c,w)
        yprime = 0
        EV%OutputTransfer =>  Arr
        call derivs(EV,EV%ScalEqsToPropagate,tau,y,yprime)
        nullify(EV%OutputTransfer)
        a = y(1)
        outputs(1:Transfer_Max, j, EV%q_ix) = Arr
        outputs(Transfer_Max+1, j, EV%q_ix) = a
        outputs(Transfer_Max+2, j, EV%q_ix) = y(2) !etak
        adotoa = 1/(y(1)*dtauda(y(1)))
        ddelta= (yprime(3)*grhoc+yprime(4)*grhob)/(grhob+grhoc)
        delta=(grhoc*y(3)+grhob*y(4))/(grhob+grhoc)
        growth= ddelta/delta/adotoa
        outputs(Transfer_Max+3, j, EV%q_ix) = adotoa !hubble
        outputs(Transfer_Max+4, j, EV%q_ix) = growth
        if (.not. EV%no_phot_multpoles) then
            outputs(Transfer_Max+5, j, EV%q_ix) = y(EV%g_ix+1) !v_g
            if (EV%TightCoupling) then
                outputs(Transfer_Max+6, j, EV%q_ix) = EV%pig
                outputs(Transfer_Max+7, j, EV%q_ix) = EV%pig/4 !just first order result
            else
                outputs(Transfer_Max+6, j, EV%q_ix) = y(EV%g_ix+2) !pi_g
                outputs(Transfer_Max+7, j, EV%q_ix) = y(EV%polind+2) !E_2
            end if
        end if
        if (.not. EV%no_nu_multpoles) then
            outputs(Transfer_Max+8, j, EV%q_ix) = y(EV%r_ix+1) !v_r
        end if

        if (global_error_flag/=0) return
    end do
    end subroutine GetOutputEvolutionFork

    function CAMB_TimeEvolution(nq, q, ntimes, times, noutputs, outputs) result(err)
    integer, intent(in) :: nq, ntimes, noutputs
    real(dl), intent(in) :: q(nq), times(ntimes)
    real(dl), intent(out) :: outputs(Transfer_Max+8, ntimes, nq)
    integer err
    integer q_ix
    Type(EvolutionVars) :: Ev

    global_error_flag = 0
    outputs = 0
    !$OMP PARALLEL DO DEFAUlT(SHARED),SCHEDUlE(DYNAMIC), PRIVATE(EV, q_ix)
    do q_ix= 1, nq
        if (global_error_flag==0) then
            EV%q_ix = q_ix
            EV%q = q(q_ix)
            EV%TransferOnly=.false.
            EV%q2=EV%q**2
            call GetNumEqns(EV)
            call GetOutputEvolutionFork(EV, times, outputs)
        end if
    end do
    !$OMP END PARALLEL DO
    err = global_error_flag
    end function CAMB_TimeEvolution

    ! END BRIDGE FOR PYTHON

    end module handles

max_stars_repo_name: Ji-Ping-Dai/CosmoCl

max_issues_repo_name: Ji-Ping-Dai/GeneralPPF

max_forks_repo_name: Ji-Ping-Dai/GeneralPPF


text:
C----------------------------------------------------------------------
        subroutine ftcpdt(iunit,ounit,status)

C       copies the data from the IUNIT CHDU to the data of the OUNIT CHDU.
C       This will overwrite any data already in the OUNIT CHDU.

C       iunit   i  fortran unit number of the input file to be copied
C       ounit   i  fortran unit number of the output file to be copied to
C       status  i  output error status
C
C       written by Wm Pence, HEASARC/GSFC, Aug 1993

        integer iunit,ounit,status

C       COMMON BLOCK DEFINITIONS:--------------------------------------------
        integer nb,ne
        parameter (nb = 20)
        parameter (ne = 200)
        integer bufnum,chdu,hdutyp,maxhdu,hdstrt,hdend,nxthdr,dtstrt
        integer nxtfld
        logical wrmode
        common/ft0001/bufnum(199),chdu(nb),hdutyp(nb),maxhdu(nb),
     &  wrmode(nb),hdstrt(nb,ne),hdend(nb),nxthdr(nb),dtstrt(nb),nxtfld
        character*1 cbuff(2880), xdummy(2880)
        common/ftheap/cbuff,xdummy
C       END OF COMMON BLOCK DEFINITIONS-----------------------------------

        integer ibuff,obuff,nblock,i

        if (status .gt. 0)return

        if (iunit .eq. ounit)then
                status=101
                return
        end if

        ibuff=bufnum(iunit)
        obuff=bufnum(ounit)

C       determine HDU structure as defined by keywords in output file
        call ftrdef(ounit,status)

C       Calculate the number of bytes to be copied.  By definition there
C       will be an integral number of 2880-byte logical blocks to be copied
        nblock=(hdstrt(ibuff,chdu(ibuff)+1)-dtstrt(ibuff))/2880

        if (nblock .gt. 0)then
C           move to the beginning of the data in the input and output files
            call ftmbyt(iunit,dtstrt(ibuff),.false.,status)
            call ftmbyt(ounit,dtstrt(obuff),.true.,status)

C           now copy the data one block at a time
            do 30 i=1,nblock
                call ftgcbf(iunit,0,2880,cbuff,status)
                call ftpcbf(ounit,0,2880,cbuff,status)
30          continue
        end if
        end

ext: f

max_stars_repo_name: iraf-community/stsdas

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
      subroutine wpot (nph, edens, imt, inrm,
     1                 rho, vclap, vcoul, vtot, ntitle, title)

c     Writes potentials to file name POTxx.DAT for each unique pot.

      implicit double precision (a-h, o-z)

      include '../HEADERS/const.h'
      include '../HEADERS/dim.h'

      dimension rho(251,0:nphx+1)
      dimension vcoul(251,0:nphx+1)
      dimension edens(251,0:nphx)
      dimension vclap(251,0:nphx)
      dimension vtot (251,0:nphx)
      dimension imt(0:nphx)
      dimension inrm(0:nphx)
      character*80 title(ntitle)

      character*30 fname
c#mn
       external rr

c     note units --
c     potentials in hartrees, so that v * 27.2 -> eV
c     density in #/(bohr)**3, so rho * e / (.529)**3 -> e/(Ang)**3

      do 180  iph = 0, nph
c        prepare file for unique potential data
         write(fname,172)  iph
  172    format('pot', i2.2, '.dat')
         open (unit=1, file=fname, status='unknown', iostat=ios)
         call chopen (ios, fname, 'wpot')
         call wthead(1, ntitle, title)
         write(1,173)  iph, imt(iph), inrm(iph)
  173    format (1x, 3i4, '  Unique potential, I_mt, I_norman.',
     1          '    Following data in atomic units.')
         write(1,*) ' iph ', iph
         write(1,174)
  174    format ('   i      r         vcoul        rho',
     1           '     ovrlp vcoul  ovrlp vtot  ovrlp rho')
c        need some limit here, 1250 points is silly.  Use
c        r <= 38, which gives 249 points with usual rgrid
         do 178  i = 1, 251
            if (rr(i) .gt. 38)  goto 179
            write(1,176) i, rr(i), vcoul(i,iph), rho(i,iph)/(4*pi),
     1                vclap(i,iph), vtot(i,iph), edens(i,iph)/(4*pi)
  176       format (1x, i4, 1p, 6e12.4)
  178    continue
  179    continue
         close(unit=1)
  180 continue

      return
      end


text:
      DOUBLE PRECISION FUNCTION DZSUM1( N, CX, INCX )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. Scalar Arguments ..
      INTEGER            INCX, N
*     ..
*     .. Array Arguments ..
      COMPLEX*16         CX( * )
*     ..
*
*  Purpose
*  =======
*
*  DZSUM1 takes the sum of the absolute values of a complex
*  vector and returns a double precision result.
*
*  Based on DZASUM from the Level 1 BLAS.
*  The change is to use the 'genuine' absolute value.
*
*  Contributed by Nick Higham for use with ZLACON.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The number of elements in the vector CX.
*
*  CX      (input) COMPLEX*16 array, dimension (N)
*          The vector whose elements will be summed.
*
*  INCX    (input) INTEGER
*          The spacing between successive values of CX.  INCX > 0.
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER            I, NINCX
      DOUBLE PRECISION   STEMP
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS
*     ..
*     .. Executable Statements ..
*
      DZSUM1 = 0.0D0
      STEMP = 0.0D0
      IF( N.LE.0 )
     $   RETURN
      IF( INCX.EQ.1 )
     $   GO TO 20
*
*     CODE FOR INCREMENT NOT EQUAL TO 1
*
      NINCX = N*INCX
      DO 10 I = 1, NINCX, INCX
*
*        NEXT LINE MODIFIED.
*
         STEMP = STEMP + ABS( CX( I ) )
   10 CONTINUE
      DZSUM1 = STEMP
      RETURN
*
*     CODE FOR INCREMENT EQUAL TO 1
*
   20 CONTINUE
      DO 30 I = 1, N
*
*        NEXT LINE MODIFIED.
*
         STEMP = STEMP + ABS( CX( I ) )
   30 CONTINUE
      DZSUM1 = STEMP
      RETURN
*
*     End of DZSUM1
*
      END

max_stars_repo_name: JamiePringle/ROMScode-PringleBBLpaper

max_issues_repo_name: JamiePringle/ROMScode-PringleBBLpaper

max_forks_repo_name: JamiePringle/ROMScode-PringleBBLpaper


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     ################################################################
c     ##                                                            ##
c     ##  subroutine field  --  get the potential energy functions  ##
c     ##                                                            ##
c     ################################################################
c
c
c     "field" sets the force field potential energy functions from
c     a parameter file and modifications specified in a keyfile
c
c
      subroutine field
      implicit none
      include 'sizes.i'
      include 'keys.i'
      include 'potent.i'
      integer i
      character*120 record
c
c
c     set the default values for the active potentials
c
      use_bond = .true.
      use_angle = .true.
      use_strbnd = .true.
      use_urey = .true.
      use_angang = .true.
      use_opbend = .true.
      use_opdist = .true.
      use_improp = .true.
      use_imptor = .true.
      use_tors = .true.
      use_pitors = .true.
      use_strtor = .true.
      use_tortor = .true.
      use_vdw = .true.
      use_charge = .true.
      use_chgdpl = .true.
      use_dipole = .true.
      use_mpole = .true.
      use_polar = .true.
      use_rxnfld = .false.
      use_solv = .true.
      use_metal = .false.
      use_geom = .true.
      use_extra = .true.
c
c     read the potential energy force field parameter file
c
      call getprm
c
c     check keywords for potential function control parameters
c
      do i = 1, nkey
         record = keyline(i)
         call prmkey (record)
      end do
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


text:
!=======================================================================
!  OpSoilPi, Subroutine, U. Singh and C.H.Porter
!-----------------------------------------------------------------------
!  Inorganic Phosphorus
!-----------------------------------------------------------------------
!  REVISION HISTORY
!  01/24/2005 CHP Written
!  06/17/2006 AJG Renamed all root-zone and no-root-zone variables with
!                 Rts and noRts.
!-----------------------------------------------------------------------
!  Called: SoilPi
!  Calls: GETLUN, HEADER, YR_DOY, SUMVALS
!=====================================================================

      SUBROUTINE OpSoilPi(CONTROL, ISWITCH,
     &  CImmobP, CMinerP, FertData, PUptake, !FracRts,
!    &  ProfLab2Act, ProfAct2Lab, ProfAct2Sta, ProfSta2Act,
     &  SOILPROP, SPi_AVAIL, SPiAvlProf, PiLabile,
     &  SPiSolProf, SPiLabProf, SPiActProf, SPiStaProf, SPiTotProf)
!    &  SPiSolRtsProf, SPiSolNoRtsProf, SPiLabRtsProf, SPiLabNoRtsProf)
!    &  PiSolRts, PiSolNoRts, PiLabRts, PiLabNoRts,
!    &  SPiLabRts, SPiLabNoRts, SPiSolRts, SPiSolNoRts)

C-----------------------------------------------------------------------
      USE ModuleDefs     !Definitions of constructed variable types,
                         ! which contain control information, soil
                         ! parameters, hourly weather data.
!     VSH
      USE CsvOutput
      USE Linklist
      IMPLICIT NONE
      SAVE

      CHARACTER*1  IDETL, IDETP, ISWPHO, RNMODE
      CHARACTER*8  LayerText(11)
      CHARACTER*10, PARAMETER :: OUTSP = 'SoilPi.OUT'

      INTEGER DAS, DOY, DYNAMIC, ERRNUM, FROP
      INTEGER L, LUN, NLAYR, REPNO, RUN, YEAR, YRDOY

      REAL SPiAvlProf, SPiSolProf, SPiLabProf
      REAL SPiActProf, SPiStaProf, SPiTotProf
!     REAL ProfLab2Act, ProfAct2Lab, ProfAct2Sta, ProfSta2Act
      Real CumFertP, CMinerP, CImmobP, CumUptakeP

      REAL, DIMENSION(NL) :: PUptake, SPi_AVAIL !, FracRts
      REAL, DIMENSION(NL) :: PiLabile
!     REAL, DIMENSION(NL) :: PiSolRts, PiSolNoRts
!     REAL, DIMENSION(NL) :: PiLabRts, PiLabNoRts
!     REAL, DIMENSION(NL) :: SPiSolRts, SPiSolNoRts
!     REAL, DIMENSION(NL) :: SPiLabRts, SPiLabNoRts

!      REAL SPiSolRtsProf, SPiSolNoRtsProf
!      REAL SPiLabRtsProf, SPiLabNoRtsProf

      LOGICAL FEXIST, DOPRINT

!     Arrays which contain data for printing in SUMMARY.OUT file
      INTEGER, PARAMETER :: SUMNUM = 4
      CHARACTER*4, DIMENSION(SUMNUM) :: LABEL
      REAL, DIMENSION(SUMNUM) :: VALUE

!-----------------------------------------------------------------------
!     Define constructed variable types based on definitions in
!     ModuleDefs.for.
      TYPE (ControlType)  CONTROL
      TYPE (FertType)     FERTDATA
      TYPE (SwitchType)   ISWITCH
      TYPE (SoilType)     SOILPROP

!     Transfer values from constructed data types into local variables.
      DAS     = CONTROL % DAS
      DYNAMIC = CONTROL % DYNAMIC
      FROP    = CONTROL % FROP
      REPNO   = CONTROL % REPNO
      RNMODE  = CONTROL % RNMODE
      RUN     = CONTROL % RUN
      YRDOY   = CONTROL % YRDOY

      CumFertP = FertData % AMTFER(P)

      NLAYR = SOILPROP % NLAYR

      FMOPT = ISWITCH % FMOPT   ! VSH

!***********************************************************************
!***********************************************************************
!     Seasonal initialization - run once per season
!***********************************************************************
      IF (DYNAMIC .EQ. SEASINIT) THEN
!-----------------------------------------------------------------------
      CMinerP    = 0.0
      CImmobP    = 0.0
      CumUptakeP = 0.0

      IDETP = ISWITCH % IDETP
      ISWPHO= ISWITCH % ISWPHO
      IDETL = ISWITCH % IDETL
      IF (IDETP .NE. 'Y' .OR. ISWPHO .NE. 'Y' .OR. IDETL == '0') RETURN

      IF (FMOPT == 'A' .OR. FMOPT == ' ') THEN   ! VSH
      CALL GETLUN(OUTSP, LUN)

      INQUIRE (FILE = OUTSP, EXIST = FEXIST)
      IF (FEXIST) THEN
        OPEN (UNIT = LUN, FILE = OUTSP, STATUS = 'OLD',
     &    IOSTAT = ERRNUM, POSITION = 'APPEND')
      ELSE
        OPEN (UNIT = LUN, FILE = OUTSP, STATUS = 'NEW',
     &    IOSTAT = ERRNUM)
        WRITE(LUN,'("*Soil Inorganic Phosphorus daily output file")')
      ENDIF
      END IF   ! VSH

      IF (RUN .EQ. 1 .OR. INDEX('QF',RNMODE) .LE. 0) THEN
        IF (FMOPT == 'A' .OR. FMOPT == ' ') THEN   ! VSH
        !For sequenced run, use replicate
        ! number instead of run number in header.
        IF (RNMODE .EQ. 'Q') THEN
          CALL HEADER(SEASINIT, LUN, REPNO)
        ELSE
          CALL HEADER(SEASINIT, LUN, RUN)
        ENDIF
        END IF   ! VSH

!       Text describing soil layer depths
        LayerText = SoilProp % LayerText
!       LayerText(11) is for layers 5 thru NLAYR
        IF (NLAYR > 4) THEN
          LayerText(5) = SOILPROP % LayerText(11)
        ENDIF

        IF (FMOPT == 'A' .OR. FMOPT == ' ') THEN   ! VSH
        WRITE (LUN,'("!",T112,"P avail (kg/ha) by soil depth (cm):",
     &    T157,"P uptake (kg/ha) by soil depth (cm):",
     &    T202,"Labile P (ppm) by soil depth (cm):")')
        WRITE (LUN, '("!",T106,15(1X,A8))')
     &    (LayerText(L),L=1,5),(LayerText(L),L=1,5),(LayerText(L),L=1,5)
        WRITE (LUN,120)
  120   FORMAT('@YEAR DOY   DAS',
     &  '     PIAD    PAVLD    PSOLD    PLABD    PACTD    PSTAD',
!    &  '  SPLabRt  SPLabNR  SPSolRt  SPSolNR',
     &  '     PAPC     PMNC     PIMC     PUPC',
     &  '    PAV1D    PAV2D    PAV3D    PAV4D    PAV5D',
     &  '    PUP1D    PUP2D    PUP3D    PUP4D    PUP5D', !,
     &  '    PLAB1    PLAB2    PLAB3    PLAB4    PLAB5')
!!       Temp variables
!     &  '    ST2AC    AC2ST    AC2LA    LA2AC',
!     &  '    FRRT1    FRRT2    FRRT3    FRRT4    FRRT5',
!     &  '   PSOLR1   PSOLR2   PSOLR3   PSOLR4   PSOLR5',
!     &  '   PSOLN1   PSOLN2   PSOLN3   PSOLN4   PSOLN5',
!     &  '   PLABR1   PLABR2   PLABR3   PLABR4   PLABR5',
!     &  '   PLABN1   PLABN2   PLABN3   PLABN4   PLABN5',
!     &  '  SPSOLR1  SPSOLR2  SPSOLR3  SPSOLR4  SPSOLR5',
!     &  '  SPSOLN1  SPSOLN2  SPSOLN3  SPSOLN4  SPSOLN5',
!     &  '  SPLABR1  SPLABR2  SPLABR3  SPLABR4  SPLABR5',
!     &  '  SPLABN1  SPLABN2  SPLABN3  SPLABN4  SPLABN5')
         END IF   ! VSH
      ENDIF

      ENDIF

!!********************************************************************
!!     Temporary output file
!      CALL GETLUN("DetailP.OUT", LUN2)
!
!      INQUIRE (FILE = "DetailP.OUT", EXIST = FEXIST)
!      IF (FEXIST) THEN
!        OPEN (UNIT = LUN2, FILE = "DetailP.OUT", STATUS = 'OLD',
!     &    IOSTAT = ERRNUM, POSITION = 'APPEND')
!      ELSE
!        OPEN (UNIT = LUN2, FILE = "DetailP.OUT", STATUS = 'NEW',
!     &    IOSTAT = ERRNUM)
!        WRITE(LUN2,'("*Soil Inorganic Phosphorus daily DETAIL file")')
!      ENDIF
!
!      IF (RUN .EQ. 1 .OR. INDEX('QF',RNMODE) .LE. 0) THEN
!
!        !For sequenced run, use replicate
!        ! number instead of run number in header.
!        IF (RNMODE .EQ. 'Q') THEN
!          CALL HEADER(SEASINIT, LUN2, REPNO)
!        ELSE
!          CALL HEADER(SEASINIT, LUN2, RUN)
!        ENDIF
!
!        WRITE (LUN2,120)
!  120   FORMAT('@YEAR DOY   DAS',
!     &  '   SPiLR1   SPiLR2   SPiLR3   SPiLR4   SPiLR5',
!     &  '   SPiLN1   SPiLN2   SPiLN3   SPiLN4   SPiLN5',
!     &  '   SPiSR1   SPiSR2   SPiSR3   SPiSR4   SPiSR5',
!     &  '   SPiSN1   SPiSN2   SPiSN3   SPiSN4   SPiSN5',
!     &  '    PiSR1    PiSR2    PiSR3    PiSR4    PiSR5',
!     &  '    PiSN1    PiSN2    PiSN3    PiSN4    PiSN5',
!     &  '      SW1      SW2      SW3      SW4      SW5',
!     &  '    FRRT1    FRRT2    FRRT3    FRRT4    FRRT5')
!      ENDIF

!      WRITE(,) (SPiLabRts(L),L=1,5), (SPiLabNoRts(L),L=1,5), SPiSolRts(L), SPiSolNoRts(L),
!     &        PiSolRts(L), PiSolNoRts(L),
!***********************************************************************
!***********************************************************************
!     OUTPUT
!***********************************************************************
      IF (DYNAMIC == OUTPUT .OR. DYNAMIC == SEASINIT) THEN
!-----------------------------------------------------------------------

      DO L = 1, NLAYR
        CumUptakeP = CumUptakeP + PUptake(L)
      ENDDO

!     Daily printout
      IF (IDETP .NE. 'Y' .OR. ISWPHO .NE. 'Y' .OR. IDETL == '0') RETURN

      DOPRINT = .FALSE.
      SELECT CASE(DYNAMIC)
      CASE (SEASINIT)
        IF (RUN == 1 .OR. INDEX('QF',RNMODE) <= 0) DOPRINT = .TRUE.
      CASE (OUTPUT)
        IF (MOD(DAS, FROP) == 0) DOPRINT = .TRUE.
      END SELECT

      IF (DOPRINT) THEN
        CALL YR_DOY(YRDOY, YEAR, DOY)
        IF (FMOPT == 'A' .OR. FMOPT == ' ') THEN   ! VSH
        WRITE (LUN,300) YEAR, DOY, DAS,
     &    SPiTotProf, SPiAvlProf, SPiSolProf,
     &    SPiLabProf, SPiActProf, SPiStaProf,
!    &    SPiLabRtsProf, SPiLabNoRtsProf, SPiSolRtsProf,SPiSolNoRtsProf,
     &    CumFertP, CMinerP, CImmobP, CumUptakeP,
     &    SPi_AVAIL(1:5), PUptake(1:5),   !,
     &    PiLabile(1:5)
!!       Temp variables:
!     &    ProfSta2Act, ProfAct2Sta, ProfAct2Lab, ProfLab2Act,
!     &    (FracRts(L),L=1,5), (PiSolRts(L),L=1,5), (PiSolNoRts(L),L=1,5)
!     &   ,(PiLabRts(L),L=1,5), (PiLabNoRts(L),L=1,5),
!     &    (SPiSolRts(L),L=1,5), (SPiSolNoRts(L),L=1,5)
!     &   ,(SPiLabRts(L),L=1,5), (SPiLabNoRts(L),L=1,5)
  300   FORMAT(1X,I4,1X,I3.3,1X,I5,
     &    F9.1, 2F9.3, 3F9.1, F9.1, 3F9.2, 10F9.3, 5F9.2)
        END IF   ! VSH
      ENDIF

!     VSH
      IF (FMOPT == 'C') THEN
         CALL CsvOutSoilPi(EXPNAME, CONTROL%RUN, CONTROL%TRTNUM,
     &CONTROL%ROTNUM, CONTROL%REPNO, YEAR, DOY, DAS,
     &SPiTotProf, SPiAvlProf, SPiSolProf, SPiLabProf, SPiActProf,
     &SPiStaProf, CumFertP, CMinerP, CImmobP, CumUptakeP, SPi_AVAIL,
     &PUptake, PiLabile,
     &vCsvlineSoilPi, vpCsvlineSoilPi, vlngthSoilPi)

         CALL LinklstSoilPi(vCsvlineSoilPi)
      END IF

!***********************************************************************
!***********************************************************************
!     SEASEND
!***********************************************************************
      ELSEIF (DYNAMIC .EQ. SEASEND) THEN
!-----------------------------------------------------------------------
!     Write end of season summary info for SUMMARY.OUT file
!     Scratch file has been opened by subroutine OPSUM, so
!     just need to retrieve correct unit number.

!     Store Summary.out labels and values in arrays to send to
!     OPSUM routines for printing.  Integers are temporarily
!     saved as real numbers for placement in real array.

      LABEL(1) = 'PI#M'; VALUE(1) = FLOAT(FertData % NAPFER(P))
      LABEL(2) = 'PICM'; VALUE(2) = CumFertP      !P applied (kg/ha)
      LABEL(3) = 'PUPC'; VALUE(3) = CumUptakeP    !Cumul P uptake(kg/ha)
      LABEL(4) = 'SPAM'; VALUE(4) = SPiTotProf    !Soil P at mat (kg/ha)

!     Send labels and values to OPSUM
      CALL SUMVALS (SUMNUM, LABEL, VALUE)

      IF (IDETP .NE. 'Y' .OR. ISWPHO .NE. 'Y' .OR. IDETL == '0') RETURN

!       Also print on last day if not already done.
        IF (MOD(DAS, FROP) .NE. 0) THEN
          CALL YR_DOY(YRDOY, YEAR, DOY)

          IF (FMOPT == 'A' .OR. FMOPT == ' ') THEN   ! VSH
          WRITE (LUN,300) YEAR, DOY, DAS,
     &    SPiTotProf, SPiAvlProf, SPiSolProf,
     &    SPiLabProf, SPiActProf, SPiStaProf,
!    &    SPiLabRtsProf, SPiLabNoRtsProf, SPiSolRtsProf,SPiSolNoRtsProf,
     &    CumFertP, CMinerP, CImmobP, CumUptakeP,
     &    SPi_AVAIL(1:5), PUptake(1:5)  !,
!!       Temp variables:
!     &    ProfSta2Act, ProfAct2Sta, ProfAct2Lab, ProfLab2Act,
!     &    (FracRts(L),L=1,5), (PiSolRts(L),L=1,5), (PiSolNoRts(L),L=1,5)
!     &   ,(PiLabRts(L),L=1,5), (PiLabNoRts(L),L=1,5),
!     &    (SPiSolRts(L),L=1,5), (SPiSolNoRts(L),L=1,5)
!     &   ,(SPiLabRts(L),L=1,5), (SPiLabNoRts(L),L=1,5)
          END IF   ! VSH
        ENDIF

!       Close daily output files.
        CLOSE (LUN)

!***********************************************************************
!***********************************************************************
!     END OF DYNAMIC IF CONSTRUCT
!***********************************************************************
      ENDIF
!***********************************************************************

      RETURN
      END SUBROUTINE OpSoilPi

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


text:
C=======================================================================
C  FOR_SEEDDM, Subroutine
C  Calculates pest damage in seed and shells.
C-----------------------------------------------------------------------
C  REVISION       HISTORY
C  01/01/1990 WDB Written
C  02/23/1998 CHP Modified for PEST Module
C  01/12/1999 GH  Incorporated into CROPGRO
C  06/06/2001 CHP Correction for 0 divide
C-----------------------------------------------------------------------
C  Called by: PEST
C  Calls:     None
C=======================================================================
      SUBROUTINE FOR_SEEDDM(
     &    DAS, LAGSD, LNGPEG, NR2, PHTIM, PHTHRS8,        !Input
     &    PSDDL, PSDDM, PSDDS, PSHDL, PSHDM, PSHDS,       !Input
     &    NSDDL, NSDDM, NSDDS, NSHDL, NSHDM, NSHDS,       !Input/Output
     &    SDNO, SHELN, SWIDOT,                            !Input/Output
     &    WSDDL, WSDDM, WSDDS,                            !Input/Output
     &    WSHDL, WSHDM, WSHDS, WTSD, WTSHE, WSHIDT,       !Input/Output
     &    CSDM, CSDN, CSHM, CSHN, SDDES, SDIDOT,          !Output
     &    SHIDOT, TSDNOL, TSDNOM, TSDNOS, TSDWTL,         !Output
     &    TSDWTM, TSDWTS, TSHNOL, TSHNOM, TSHNOS,         !Output
     &    TSHWTL, TSHWTM, TSHWTS,                         !Output
     &    DYNAMIC, WSDD, PSDD,SDWT)                       !Control

!-----------------------------------------------------------------------
      USE ModuleDefs     !Definitions of constructed variable types,
        ! which contain control information, soil
        ! parameters, hourly weather data.
!     NCOHORTS defined in ModuleDefs.for

      IMPLICIT NONE
      SAVE

      INTEGER I, NR2, NPP, DAS, DYNAMIC

      REAL PHTIM(365)
      REAL PAGE,DEST,WPS
      REAL WTSD(NCOHORTS) ,WTSHE(NCOHORTS)
      REAL SDNO(NCOHORTS), SHELN(NCOHORTS)
      REAL SWIDOT, WSHIDT, SDIDOT, SHIDOT, LAGSD, LNGPEG
C     Seed Variables
      REAL NSDDS, NSDDL, NSDDM, PSDDS, PSDDL, PSDDM, WSDDS, WSDDL, WSDDM
      REAL TSDNO, TSDWT, CSDM, CSDN
      REAL TSDNOS, TSDNOL, TSDNOM, TSDWTS, TSDWTL, TSDWTM
      REAL SDDES(NCOHORTS), SDWDES(NCOHORTS), SDNDES(NCOHORTS)
      REAL WSDD,PSDD,SDWT
C     Shell Variables
      REAL NSHDS,NSHDL,NSHDM,PSHDS,PSHDL,PSHDM,WSHDS,WSHDL,WSHDM
      REAL TSHNO, TSHWT, CSHM, CSHN, PHTHRS8
      REAL TSHNOS,TSHNOL,TSHNOM,TSHWTS,TSHWTL,TSHWTM
      REAL SHWDES(NCOHORTS), SHNDES(NCOHORTS)

C      SAVE TSDNO, TSDWT, TSHNO, TSHWT

!***********************************************************************
!***********************************************************************
!     Seasonal initialization - run once per season
!***********************************************************************
      IF (DYNAMIC .EQ. SEASINIT) THEN
!-----------------------------------------------------------------------
      DO I=1, NCOHORTS
        SDWDES(I) = 0.0
        SDNDES(I) = 0.0
        SHWDES(I) = 0.0
        SHNDES(I) = 0.0
        SDDES(I)  = 0.0
      ENDDO

      SDIDOT = 0.0
      SHIDOT = 0.0
      SWIDOT = 0.0
      WSHIDT = 0.0

      TSDWTS = 0.0
      TSDNOS = 0.0
      TSHWTS = 0.0
      TSHNOS = 0.0

      TSDWTL = 0.0
      TSDNOL = 0.0
      TSHWTL = 0.0
      TSHNOL = 0.0

      TSDWTM = 0.0
      TSDNOM = 0.0
      TSHWTM = 0.0
      TSHNOM = 0.0

      CSDM = 0.
      CSDN = 0.
      CSHM = 0.
      CSHN = 0.

      TSDNO = 0.
      TSDWT = 0.
      TSHNO = 0.
      TSHWT = 0.

!***********************************************************************
!***********************************************************************
!     Daily rate calculations
!***********************************************************************
      ELSEIF (DYNAMIC .EQ. RATE) THEN
!-----------------------------------------------------------------------
      SDIDOT = 0.0
      SHIDOT = 0.0
      SWIDOT = 0.0
      WSHIDT = 0.0

      SDWDES = 0.0    !
      SDNDES = 0.0    !Arrays
      SHWDES = 0.0    !
      SHNDES = 0.0    !

      TSDWTS = 0.0
      TSDNOS = 0.0
      TSHWTS = 0.0
      TSHNOS = 0.0

      TSDWTL = 0.0
      TSDNOL = 0.0
      TSHWTL = 0.0
      TSHNOL = 0.0

      TSDWTM = 0.0
      TSDNOM = 0.0
      TSHWTM = 0.0
      TSHNOM = 0.0

!-----------------------------------------------------------------------
!     Compute total seed and shell mass and number for small, large and
!         mature pod classes.
!-----------------------------------------------------------------------
      DO NPP = 1, DAS-NR2-1
        PAGE = PHTIM(DAS-NR2) - PHTIM(NPP)
!-----------------------------------------------------------------------
C     Count small seed cohorts less than LAGSD (R3) old
C-----------------------------------------------------------------------
        IF (WTSD(NPP) .LT. 0.001 .AND. SDNO(NPP) .GT. 0.0) THEN
        TSDWTS = TSDWTS + WTSD(NPP)
        TSDNOS = TSDNOS + SDNO(NPP)
        TSHWTS = TSHWTS + WTSHE(NPP)
        TSHNOS = TSHNOS + SHELN(NPP)
        ENDIF
C-----------------------------------------------------------------------
C     Count large seed cohorts between small and full pod (R6)
C-----------------------------------------------------------------------
        IF (WTSD(NPP) .GT. 0.0 .AND. PAGE .GE. LAGSD .AND.
     &    PAGE .LT. PHTHRS8) THEN
        TSDWTL = TSDWTL + WTSD(NPP)
        TSDNOL = TSDNOL + SDNO(NPP)
        TSHWTL = TSHWTL + WTSHE(NPP)
        TSHNOL = TSHNOL + SHELN(NPP)
        ENDIF
C-----------------------------------------------------------------------
C    Count mature cohorts that are older than full pod (R6)
C-----------------------------------------------------------------------
        IF (PAGE .GT. PHTHRS8) THEN
        TSDWTM = TSDWTM + WTSD(NPP)
        TSDNOM = TSDNOM + SDNO(NPP)
        TSHWTM = TSHWTM + WTSHE(NPP)
        TSHNOM = TSHNOM + SHELN(NPP)
        ENDIF
      ENDDO

      TSDNOS = MAX(0.,TSDNOS)
      TSDWTS = MAX(0.,TSDWTS)
      TSHNOS = MAX(0.,TSHNOS)
      TSHWTS = MAX(0.,TSHWTS)

      TSDNOL = MAX(0.,TSDNOL)
      TSDWTL = MAX(0.,TSDWTL)
      TSHNOL = MAX(0.,TSHNOL)
      TSHWTL = MAX(0.,TSHWTL)

      TSDNOM = MAX(0.,TSDNOM)
      TSDWTM = MAX(0.,TSDWTM)
      TSHNOM = MAX(0.,TSHNOM)
      TSHWTM = MAX(0.,TSHWTM)
!-----------------------------------------------------------------------
C      Compute total seed and shell number and mass in the damaged
C      size class (These values not used - chp)
C---------------------------------------------------------------------
      TSDNO = TSDNOS + TSDNOL + TSDNOM
      TSDWT = TSDWTS + TSDWTL + TSDWTM
      TSHNO = TSHNOS + TSHNOL + TSHNOM
      TSHWT = TSHWTS + TSHWTL + TSHWTM

!-----------------------------------------------------------------------
!     Loop through cohorts and apply damage
!-----------------------------------------------------------------------
      DO NPP=1,DAS-NR2-1
        PAGE = PHTIM(DAS-NR2) - PHTIM(NPP)

C-----------------------------------------------------------------------
C    Apply damage as percent of seed number or seed mass
C-----------------------------------------------------------------------
        IF (PSDDS .GT. 0.0 .AND. WTSD(NPP) .LT. 0.0001 .AND.
     &    SDNO(NPP) .GT. 0.0) THEN
        SDWDES(NPP) = SDWDES(NPP) + PSDDS * WTSD(NPP) / 100.0
        SDNDES(NPP) = SDNDES(NPP) + PSDDS * SDNO(NPP) / 100.0
        IF (PSHDS .LT. 0.001 .AND. NSHDS .LT. 0.001) THEN
        SDDES(NPP) = SDDES(NPP) + PSDDS * SDNO(NPP) / 100.0
        ENDIF
        ENDIF

        IF (PSDDL .GT. 0.0 .AND. PAGE .GE. LAGSD .AND.
     &    PAGE .LT. PHTHRS8) THEN
        SDWDES(NPP) = SDWDES(NPP) + PSDDL * WTSD(NPP) / 100.0
        SDNDES(NPP) = SDNDES(NPP) + PSDDL * SDNO(NPP) / 100.0
        IF (PSHDL .LT. 0.001 .AND. NSHDL .LT. 0.001) THEN
        SDDES(NPP) = SDDES(NPP) + PSDDL * SDNO(NPP)/100.0
        ENDIF
        ENDIF

        IF (PSDDM .GT. 0.0 .AND. PAGE .GT. PHTHRS8) THEN
        SDWDES(NPP) = SDWDES(NPP) + PSDDM * WTSD(NPP) / 100.0
        SDNDES(NPP) = SDNDES(NPP) + PSDDM * SDNO(NPP) / 100.0
        IF (PSHDM .LT. 0.001 .AND. NSHDM .LT. 0.001) THEN
        SDDES(NPP) = SDDES(NPP)+ PSDDM * SDNO(NPP)/100.0
        ENDIF
        ENDIF

C-----------------------------------------------------------------------
C     Apply damage as number of seed/m2
C-----------------------------------------------------------------------
        IF (NSDDS .GT. 0.0 .AND. WTSD(NPP) .LT. 0.0001 .AND.
     &  SDNO(NPP) .GT. 0.0 .AND. TSDNOS .GT. 0.0) THEN
        NSDDS = MIN(NSDDS,TSDNOS)
        SDWDES(NPP) = SDWDES(NPP) + WTSD(NPP) * NSDDS/TSDNOS
        SDNDES(NPP) = SDNDES(NPP) + SDNO(NPP) * NSDDS/TSDNOS
        IF (PSHDS .LT. 0.001 .AND. NSHDS .LT. 0.001) THEN
        SDDES(NPP) = SDDES(NPP) +  NSDDS * SDNO(NPP)/TSDNOS
        ENDIF
        ENDIF

        IF (NSDDL .GT. 0.0 .AND. PAGE .GE. LAGSD .AND.
     &   PAGE .LT. PHTHRS8 .AND. TSDNOL .GT. 0.0) THEN
        NSDDL = MIN(NSDDL,TSDNOL)
        SDWDES(NPP) = SDWDES(NPP) + WTSD(NPP) * NSDDL/TSDNOL
        SDNDES(NPP) = SDNDES(NPP) + SDNO(NPP) * NSDDL/TSDNOL
        IF (PSHDL .LT. 0.001 .AND. NSHDM .LT. 0.001) THEN
        SDDES(NPP) = SDDES(NPP) +  NSDDL * SDNO(NPP)/TSDNOL
        ENDIF
        ENDIF

        IF (NSDDM .GT. 0.0 .AND. PAGE .GT. PHTHRS8 .AND.
     &    TSDNOM .GT. 0.0) THEN
        NSDDM = MIN(NSDDM,TSDNOM)
        SDWDES(NPP) = SDWDES(NPP) + WTSD(NPP) * NSDDM/TSDNOM
        SDNDES(NPP) = SDNDES(NPP) + SDNO(NPP) * NSDDM/TSDNOM
        IF (PSHDM .LT. 0.001 .AND. NSHDM .LT. 0.001) THEN
        SDDES(NPP) = SDDES(NPP) + NSDDM * SDNO(NPP)/TSDNOM
        ENDIF
        ENDIF

C-----------------------------------------------------------------------
C     Apply damage as mass of seed
C-----------------------------------------------------------------------
        IF (WTSD(NPP) .LE. 0.001) GOTO 1000
        WPS = SDNO(NPP)/WTSD(NPP)

        IF (WSDDS .GT. 0.0 .AND. PAGE .LT. LAGSD .AND.
     &  PAGE .GT. LNGPEG .AND. TSDNOS .GT. 0.0) THEN
        WSDDS = MIN(WSDDS,TSDWTS)
        SDWDES(NPP) = SDWDES(NPP) + WSDDS*SDNO(NPP)/TSDNOS
        SDNDES(NPP) = SDNDES(NPP) + WSDDS*SDNO(NPP)*WPS/TSDNOS
        IF (PSHDS .GT. 0.0 .AND. NSHDS .GT. 0.0 .AND.
     &    TSDWTS .GT. 0.0001) THEN
        SDDES(NPP) = SDDES(NPP) + WPS*WSDDS*SDNO(NPP)/TSDNOS
        ENDIF
        ENDIF

        IF (WSDDL .GT. 0.0 .AND. PAGE .GE. LAGSD .AND.
     &  PAGE .LT. PHTHRS8 .AND. TSDNOL .GT. 0.0) THEN
        WSDDL = MIN(WSDDL,TSDWTL)
        SDWDES(NPP) = SDWDES(NPP) + WSDDL*SDNO(NPP)/TSDNOL
        SDNDES(NPP) = SDNDES(NPP) + WSDDL*SDNO(NPP)*WPS/TSDNOL
        IF (PSHDL .GT. 0.0 .AND. NSHDL .GT. 0.0) THEN
        SDDES(NPP) = SDDES(NPP) + WPS*WSDDL*SDNO(NPP)/TSDNOL
        ENDIF
        ENDIF

        IF (WSDDM .GT. 0.0 .AND. PAGE .GT. PHTHRS8 .AND.
     &  TSDNOM .GT. 0.0) THEN
        WSDDM = MIN(WSDDM,TSDWTM)
        SDWDES(NPP) = SDWDES(NPP) + WSDDM*SDNO(NPP)/TSDNOM
        SDNDES(NPP) = SDNDES(NPP) + WSDDM*SDNO(NPP)*WPS/TSDNOM
        IF (PSHDM .GT. 0.0 .AND. NSHDM .GT. 0.0) THEN
        SDDES(NPP) = SDDES(NPP) + WPS*WSDDM*SDNO(NPP)/TSDNOM
        ENDIF
        ENDIF
 1000   CONTINUE

C-----------------------------------------------------------------------
C        Apply damage as percent shell mass
C-----------------------------------------------------------------------
        IF (PSHDS .GT. 0.0 .AND. PAGE .GT. LNGPEG .AND.
     &  PAGE .LT. LAGSD) THEN
        SHWDES(NPP) = SHWDES(NPP) + PSHDS * WTSHE(NPP) / 100.0
        SHNDES(NPP) = SHNDES(NPP) + PSHDS * SHELN(NPP) / 100.0
        ENDIF

        IF(PSHDL .GT. 0.0 .AND. PAGE .GE. LAGSD .AND.
     &  PAGE .LT. PHTHRS8) THEN
        SHWDES(NPP) = SHWDES(NPP) + PSHDL * WTSHE(NPP) / 100.0
        SHNDES(NPP) = SHNDES(NPP) + PSHDL * SHELN(NPP) / 100.0
        ENDIF

        IF(PSHDM .GT. 0.0 .AND. PAGE .GE. PHTHRS8) THEN
        SHWDES(NPP) = SHWDES(NPP) + PSHDM * WTSHE(NPP) / 100.0
        SHNDES(NPP) = SHNDES(NPP) + PSHDM * SHELN(NPP) / 100.0
        ENDIF

C-----------------------------------------------------------------------
C     Apply damage as number of shells /m2
C-----------------------------------------------------------------------
        IF (NSHDS .GT. 0.0 .AND. PAGE .GT. LNGPEG .AND.
     &  PAGE .LT. LAGSD .AND. TSHNOS .GT. 0.0) THEN
        NSHDS = MIN(NSHDS,TSHNOS)
        SHWDES(NPP) = SHWDES(NPP) + WTSHE(NPP)*NSHDS/TSHNOS
        SHNDES(NPP) = SHNDES(NPP) + SHELN(NPP)*NSHDS/TSHNOS
        ENDIF

        IF (NSHDL .GT. 0.0 .AND. PAGE .GT. LAGSD .AND.
     &  PAGE .LT. PHTHRS8 .AND. TSHNOL .GT. 0) THEN
        NSHDL = MIN(NSHDL,TSHNOL)
        SHWDES(NPP) = SHWDES(NPP) + WTSHE(NPP)*NSHDL/TSHNOL
        SHNDES(NPP) = SHNDES(NPP) + SHELN(NPP)*NSHDL/TSHNOL
        ENDIF

        IF (NSHDM .GT. 0.0 .AND. PAGE .GT. PHTHRS8 .AND.
     &  TSHNOM .GT. 0.0) THEN
        NSHDM = MIN(NSHDM,TSHNOM)
        SHWDES(NPP) = SHWDES(NPP) + WTSHE(NPP)*NSHDM/TSHNOM
        SHNDES(NPP) = SHNDES(NPP) + SHELN(NPP)*NSHDM/TSHNOM
        ENDIF

C-----------------------------------------------------------------------
C     Apply damage as weight of shell /m2
C-----------------------------------------------------------------------
        IF(WSHDS .GT. 0.0 .AND. PAGE .GT. LNGPEG .AND.
     &  PAGE .LT. LAGSD .AND. TSHNOS .GT. 0.0) THEN
        WSHDS = MIN(WSHDS,TSHWTS)
        SHWDES(NPP) = SHWDES(NPP) + WSHDS*WTSHE(NPP)/TSHWTS
        DEST = WSHDS * WTSHE(NPP)/TSHWTS
        SHNDES(NPP) = SHNDES(NPP) + DEST*SHELN(NPP)/WTSHE(NPP)
        ENDIF

        IF (WSHDL .GT. 0.0 .AND. PAGE .GT. LAGSD .AND.
     &  PAGE .LT. PHTHRS8 .AND. TSHNOL .GT. 0.0) THEN
        WSHDL = MIN(WSHDL,TSHWTL)
        SHWDES(NPP) = SHWDES(NPP) + WSHDL*WTSHE(NPP)/TSHWTL
        DEST = WSHDL * WTSHE(NPP)/TSHWTL
        SHNDES(NPP) = SHNDES(NPP) + DEST*SHELN(NPP)/WTSHE(NPP)
        ENDIF

        IF (WSHDM .GT. 0.0 .AND. PAGE .GT. PHTHRS8 .AND.
     &  TSHNOM .GT. 0.0) THEN
        WSHDM = MIN(WSHDM,TSHWTM)
        SHWDES(NPP) = SHWDES(NPP) + WSHDM*WTSHE(NPP)/TSHWTM
        DEST = WSHDM * WTSHE(NPP)/TSHWTM
        SHNDES(NPP) = SHNDES(NPP) + DEST*SHELN(NPP)/WTSHE(NPP)
        ENDIF

!       Check that mass and number destroyed are not greater than
!             actual mass and number.
        SDWDES(NPP) = MIN(SDWDES(NPP), WTSD(NPP))
        SDNDES(NPP) = MIN(SDNDES(NPP), SDNO(NPP))
        SHWDES(NPP) = MIN(SHWDES(NPP), WTSHE(NPP))
        SHNDES(NPP) = MIN(SHNDES(NPP), SHELN(NPP))

!       Update daily total mass and number of seed and shell destroyed.
        SWIDOT = SWIDOT + SDWDES(NPP)
        SDIDOT = SDIDOT + SDNDES(NPP)
        WSHIDT = WSHIDT + SHWDES(NPP)
        SHIDOT = SHIDOT + SHNDES(NPP)

      ENDDO       !End of cohort loop

      !----------------------------------------------------------------
      ! This section was added to compute seed damage in CERES MAIZE.
      ! The main difference is that CERES assumes all seeds are the
      ! same age, thus there are no seed chorts.
      !
      ! Apply damage as absolute (g seed/m2/d)
      IF(WSDD.GT.0) THEN
        SWIDOT = SWIDOT + WSDD
      ENDIF

      ! Apply damage as percent (%/d)
      IF (PSDD .GT. 0) THEN
      SWIDOT = SWIDOT + SDWT*(PSDD/100)
      ENDIF

!***********************************************************************
!***********************************************************************
!     Daily integration
!***********************************************************************
      ELSEIF (DYNAMIC .EQ. INTEGR) THEN
C-----------------------------------------------------------------------
      CSDM = CSDM + SWIDOT
      CSDN = CSDN + SDIDOT
      CSHM = CSHM + WSHIDT
      CSHN = CSHN + SHIDOT

      DO NPP = 1, DAS - NR2 - 1
        PAGE = PHTIM(DAS - NR2) - PHTIM(NPP)
        WTSD(NPP)  = MAX(0., WTSD(NPP)  - SDWDES(NPP))
        SDNO(NPP)  = MAX(0., SDNO(NPP)  - SDNDES(NPP))
        WTSHE(NPP) = MAX(0., WTSHE(NPP) - SHWDES(NPP))
        SHELN(NPP) = MAX(0., SHELN(NPP) - SHNDES(NPP))
      ENDDO

!***********************************************************************
      ENDIF
!***********************************************************************
!***********************************************************************
!     END OF DYNAMIC IF CONSTRUCT
!***********************************************************************
      END   ! SUBROUTINE FOR_SEEDDM

!-----------------------------------------------------------------------
!     Variable Definitions
!-----------------------------------------------------------------------
! CSDM     Cumulative seed mass destroyed (g/m2)
! CSDN     Cumulative number of seeds destroyed (#/m2)
! CSHM     Cumulative shell mass destroyed (g/m2)
! CSHN     Cumulative number of shells destroyed (#/m2)
! DAS      Days after start of simulation (d)
! DEST     Percentage of shells destroyed
! DYNAMIC  Module control variable; =RUNINIT, SEASINIT, RATE, EMERG,
!            INTEGR, OUTPUT, or SEASEND
! LAGSD    Time required between shell growth and seed growth, per cohort
!            (Photo-thermal days)
! LNGPEG   Time between start of peg (full flower) and shell formation (for
!            peanuts only).  Defines slow growth period. (Photo-thermal days)
! NPP      Cohort number used as index in loops
! NR2      Day when 50% of plants have one peg (peanuts only) (d)
! NSDDL    Daily number of large seeds damaged (#/m2/d)
! NSDDM    Daily number of mature seeds damaged (#/m2/d)
! NSDDS    Daily number of small seeds damaged (#/m2/d)
! NSHDL    Daily number of large shells damaged (#/m2/d)
! NSHDM    Daily number of mature shells damaged (#/m2/d)
! NSHDS    Daily number of small shells damaged (#/m2/d)
! PAGE     Photothermal age of each cohort  (Photo-thermal days)
! PHTHRS8  Threshold time that must accumulate in phase 8 for the next
!            stage to occur.  Equivalent to PHTHRS(8) in Subroutine PHENOLOG.
!            (photothermal days)
! PHTIM    Cumulative photothermal time ages of seeds and shells
! PSDDL    Percent large seed mass damage (%/d)
! PSDDM    Percent mature seed mass damage (%/d)
! PSDDS    Percent small seed mass damage (%/d)
! PSHDL    Percent large shell mass damage (%/d)
! PSHDM    Percent mature shell mass damage (%/d)
! PSHDS    Percent small shell mass damage (%/d)
! SDDES(J) Number of seeds destroyed today in cohort J when shells are
!            not destroyed (#/m2/d)
! SDIDOT   Number of seeds destroyed on the current day (#/m2/d)
! SDNDES(J)Number of seeds destroyed in cohort J by pests or disease
!            (#/m2/day)
! SDNO(J)  Number of seeds for cohort J (#/m2)
! SDWT     Seed weight, g/m2
! SDWDES(J)Mass of seeds in cohort J destroyed today by pest or disease
!            (g/m2/d)
! SHELN(J) Number of shells for cohort J (#/m2)
! SHIDOT   Number of shells destroyed on current day (#/m2/d)
! SWIDOT   Daily seed mass damage (g/m2/day)
! TSDNO    Total seed number in all cohorts (#/m2)
! TSDNOL   Total number of large seeds (#/m2)
! TSDNOM   Total number of mature seeds (#/m2)
! TSDNOS   Total number of small seeds (#/m2)
! TSDWT    Total seed mass in all cohorts (g/m2)
! TSDWTL   Seed mass for large seeds (g/m2)
! TSDWTM   Seed mass for mature seeds (g/m2)
! TSDWTS   Seed mass for small seeds (g/m2)
! TSHNO    Total shell number in all cohorts (#/m2)
! TSHNOL   Number shells with large seeds (#/m2)
! TSHNOM   Number shells with mature seeds (#/m2)
! TSHNOS   Number shells with small seeds (#/m2)
! TSHWT    Total shell mass in all cohorts (g/m2)
! TSHWTL   Shell mass with large seeds (g/m2)
! TSHWTM   Shell mass with mature seeds (g/m2)
! TSHWTS   Shell mass with small seeds (g/m2)
! WPS      Number seeds per gram (#/g)
! WSDDL    Daily mass of large seed damaged (g/m2/day)
! WSDDM    Daily mass of mature seed damaged (g/m2/day)
! WSDDS    Daily mass of small seed damaged (g/m2/day)
! WSHDL    Daily mass of large shell damaged (g/m2/day)
! WSHDM    Daily mass of mature shell damaged (g/m2/day)
! WSHDS    Daily mass of small shell damaged (g/m2/day)
! WSHIDT   Weight of shell tissue consumed by pests today (g[shell]/m2-d)
! WTSD(J)  Seed mass  for cohort J (g/m2)
! WTSHE(J) Shell mass  for cohort J (g/m2)
!-----------------------------------------------------------------------
!     End Subroutine FOR_SEEDDM
!-----------------------------------------------------------------------

max_stars_repo_name: vstocca/dssat-csm-os

max_issues_repo_name: vstocca/dssat-csm-os

max_forks_repo_name: vstocca/dssat-csm-os


text:
c---------------------------------------------------------------------
c---------------------------------------------------------------------

       subroutine x_solve(rho_i, us, speed, rhs, nx, nxmax, ny, nz)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c this function performs the solution of the approximate factorization
c step in the x-direction for all five matrix components
c simultaneously. The Thomas algorithm is employed to solve the
c systems for the x-lines. Boundary conditions are non-periodic
c---------------------------------------------------------------------

       include 'header.h'

       integer nx, nxmax, ny, nz
       double precision rho_i(  0:nxmax-1,0:ny-1,0:nz-1),
     $                  us   (  0:nxmax-1,0:ny-1,0:nz-1),
     $                  speed(  0:nxmax-1,0:ny-1,0:nz-1),
     $                  rhs  (5,0:nxmax-1,0:ny-1,0:nz-1)

       integer i, j, k, i1, i2, m
       double precision  ru1, fac1, fac2

c---------------------------------------------------------------------
c---------------------------------------------------------------------

       if (timeron) call timer_start(t_xsolve)
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(fac2,m,fac1,i2,i1,ru1,i,j,k)
!$OMP&  SHARED(comz6,comz1,comz4,comz5,c2dttx1,dttx1,dttx2,dx1,dxmax,
!$OMP& c1c5,dx5,con43,dx2,c3c4,nx,ny,nz)
       do  k = 1, nz-2

c---------------------------------------------------------------------
c Computes the left hand side for the three x-factors
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     zap the whole left hand side for starters
c---------------------------------------------------------------------
          do  j = 1, ny-2
             do   m = 1, 5
                lhs (m,0,j) = 0.0d0
                lhsp(m,0,j) = 0.0d0
                lhsm(m,0,j) = 0.0d0
                lhs (m,nx-1,j) = 0.0d0
                lhsp(m,nx-1,j) = 0.0d0
                lhsm(m,nx-1,j) = 0.0d0
             end do

c---------------------------------------------------------------------
c      next, set all diagonal values to 1. This is overkill, but
c      convenient
c---------------------------------------------------------------------
             lhs (3,0,j) = 1.0d0
             lhsp(3,0,j) = 1.0d0
             lhsm(3,0,j) = 1.0d0
             lhs (3,nx-1,j) = 1.0d0
             lhsp(3,nx-1,j) = 1.0d0
             lhsm(3,nx-1,j) = 1.0d0
          end do

c---------------------------------------------------------------------
c      first fill the lhs for the u-eigenvalue
c---------------------------------------------------------------------
          do  j = 1, ny-2
             do  i = 0, nx-1
                ru1 = c3c4*rho_i(i,j,k)
                cv(i) = us(i,j,k)
                rhon(i) = dmax1(dx2+con43*ru1,
     >                          dx5+c1c5*ru1,
     >                          dxmax+ru1,
     >                          dx1)
             end do

             do  i = 1, nx-2
                lhs(1,i,j) =   0.0d0
                lhs(2,i,j) = - dttx2 * cv(i-1) - dttx1 * rhon(i-1)
                lhs(3,i,j) =   1.0d0 + c2dttx1 * rhon(i)
                lhs(4,i,j) =   dttx2 * cv(i+1) - dttx1 * rhon(i+1)
                lhs(5,i,j) =   0.0d0
             end do
          end do

c---------------------------------------------------------------------
c      add fourth order dissipation
c---------------------------------------------------------------------

          do  j = 1, ny-2
             i = 1
             lhs(3,i,j) = lhs(3,i,j) + comz5
             lhs(4,i,j) = lhs(4,i,j) - comz4
             lhs(5,i,j) = lhs(5,i,j) + comz1

             i = 2
             lhs(2,i,j) = lhs(2,i,j) - comz4
             lhs(3,i,j) = lhs(3,i,j) + comz6
             lhs(4,i,j) = lhs(4,i,j) - comz4
             lhs(5,i,j) = lhs(5,i,j) + comz1
          end do

          do  j = 1, ny-2
             do   i=3, nx-4
                lhs(1,i,j) = lhs(1,i,j) + comz1
                lhs(2,i,j) = lhs(2,i,j) - comz4
                lhs(3,i,j) = lhs(3,i,j) + comz6
                lhs(4,i,j) = lhs(4,i,j) - comz4
                lhs(5,i,j) = lhs(5,i,j) + comz1
             end do
          end do

          do  j = 1, ny-2
             i = nx-3
             lhs(1,i,j) = lhs(1,i,j) + comz1
             lhs(2,i,j) = lhs(2,i,j) - comz4
             lhs(3,i,j) = lhs(3,i,j) + comz6
             lhs(4,i,j) = lhs(4,i,j) - comz4

             i = nx-2
             lhs(1,i,j) = lhs(1,i,j) + comz1
             lhs(2,i,j) = lhs(2,i,j) - comz4
             lhs(3,i,j) = lhs(3,i,j) + comz5
          end do

c---------------------------------------------------------------------
c      subsequently, fill the other factors (u+c), (u-c) by adding to
c      the first
c---------------------------------------------------------------------
          do  j = 1, ny-2
             do   i = 1, nx-2
                lhsp(1,i,j) = lhs(1,i,j)
                lhsp(2,i,j) = lhs(2,i,j) -
     >                            dttx2 * speed(i-1,j,k)
                lhsp(3,i,j) = lhs(3,i,j)
                lhsp(4,i,j) = lhs(4,i,j) +
     >                            dttx2 * speed(i+1,j,k)
                lhsp(5,i,j) = lhs(5,i,j)
                lhsm(1,i,j) = lhs(1,i,j)
                lhsm(2,i,j) = lhs(2,i,j) +
     >                            dttx2 * speed(i-1,j,k)
                lhsm(3,i,j) = lhs(3,i,j)
                lhsm(4,i,j) = lhs(4,i,j) -
     >                            dttx2 * speed(i+1,j,k)
                lhsm(5,i,j) = lhs(5,i,j)
             end do
          end do

c---------------------------------------------------------------------
c                          FORWARD ELIMINATION
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c      perform the Thomas algorithm; first, FORWARD ELIMINATION
c---------------------------------------------------------------------

          do  j = 1, ny-2
             do    i = 0, nx-3
                i1 = i  + 1
                i2 = i  + 2
                fac1      = 1.d0/lhs(3,i,j)
                lhs(4,i,j)  = fac1*lhs(4,i,j)
                lhs(5,i,j)  = fac1*lhs(5,i,j)
                do    m = 1, 3
                   rhs(m,i,j,k) = fac1*rhs(m,i,j,k)
                end do
                lhs(3,i1,j) = lhs(3,i1,j) -
     >                         lhs(2,i1,j)*lhs(4,i,j)
                lhs(4,i1,j) = lhs(4,i1,j) -
     >                         lhs(2,i1,j)*lhs(5,i,j)
                do    m = 1, 3
                   rhs(m,i1,j,k) = rhs(m,i1,j,k) -
     >                         lhs(2,i1,j)*rhs(m,i,j,k)
                end do
                lhs(2,i2,j) = lhs(2,i2,j) -
     >                         lhs(1,i2,j)*lhs(4,i,j)
                lhs(3,i2,j) = lhs(3,i2,j) -
     >                         lhs(1,i2,j)*lhs(5,i,j)
                do    m = 1, 3
                   rhs(m,i2,j,k) = rhs(m,i2,j,k) -
     >                         lhs(1,i2,j)*rhs(m,i,j,k)
                end do
             end do
          end do

c---------------------------------------------------------------------
c      The last two rows in this zone are a bit different,
c      since they do not have two more rows available for the
c      elimination of off-diagonal entries
c---------------------------------------------------------------------

          do  j = 1, ny-2
             i  = nx-2
             i1 = nx-1
             fac1      = 1.d0/lhs(3,i,j)
             lhs(4,i,j)  = fac1*lhs(4,i,j)
             lhs(5,i,j)  = fac1*lhs(5,i,j)
             do    m = 1, 3
                rhs(m,i,j,k) = fac1*rhs(m,i,j,k)
             end do
             lhs(3,i1,j) = lhs(3,i1,j) -
     >                      lhs(2,i1,j)*lhs(4,i,j)
             lhs(4,i1,j) = lhs(4,i1,j) -
     >                      lhs(2,i1,j)*lhs(5,i,j)
             do    m = 1, 3
                rhs(m,i1,j,k) = rhs(m,i1,j,k) -
     >                      lhs(2,i1,j)*rhs(m,i,j,k)
             end do
c---------------------------------------------------------------------
c            scale the last row immediately
c---------------------------------------------------------------------
             fac2             = 1.d0/lhs(3,i1,j)
             do    m = 1, 3
                rhs(m,i1,j,k) = fac2*rhs(m,i1,j,k)
             end do
          end do

c---------------------------------------------------------------------
c      do the u+c and the u-c factors
c---------------------------------------------------------------------

          do  j = 1, ny-2
             do    i = 0, nx-3
                i1 = i  + 1
                i2 = i  + 2
                m = 4
                fac1       = 1.d0/lhsp(3,i,j)
                lhsp(4,i,j)  = fac1*lhsp(4,i,j)
                lhsp(5,i,j)  = fac1*lhsp(5,i,j)
                rhs(m,i,j,k) = fac1*rhs(m,i,j,k)
                lhsp(3,i1,j) = lhsp(3,i1,j) -
     >                        lhsp(2,i1,j)*lhsp(4,i,j)
                lhsp(4,i1,j) = lhsp(4,i1,j) -
     >                        lhsp(2,i1,j)*lhsp(5,i,j)
                rhs(m,i1,j,k) = rhs(m,i1,j,k) -
     >                        lhsp(2,i1,j)*rhs(m,i,j,k)
                lhsp(2,i2,j) = lhsp(2,i2,j) -
     >                        lhsp(1,i2,j)*lhsp(4,i,j)
                lhsp(3,i2,j) = lhsp(3,i2,j) -
     >                        lhsp(1,i2,j)*lhsp(5,i,j)
                rhs(m,i2,j,k) = rhs(m,i2,j,k) -
     >                        lhsp(1,i2,j)*rhs(m,i,j,k)
                m = 5
                fac1       = 1.d0/lhsm(3,i,j)
                lhsm(4,i,j)  = fac1*lhsm(4,i,j)
                lhsm(5,i,j)  = fac1*lhsm(5,i,j)
                rhs(m,i,j,k) = fac1*rhs(m,i,j,k)
                lhsm(3,i1,j) = lhsm(3,i1,j) -
     >                        lhsm(2,i1,j)*lhsm(4,i,j)
                lhsm(4,i1,j) = lhsm(4,i1,j) -
     >                        lhsm(2,i1,j)*lhsm(5,i,j)
                rhs(m,i1,j,k) = rhs(m,i1,j,k) -
     >                        lhsm(2,i1,j)*rhs(m,i,j,k)
                lhsm(2,i2,j) = lhsm(2,i2,j) -
     >                        lhsm(1,i2,j)*lhsm(4,i,j)
                lhsm(3,i2,j) = lhsm(3,i2,j) -
     >                        lhsm(1,i2,j)*lhsm(5,i,j)
                rhs(m,i2,j,k) = rhs(m,i2,j,k) -
     >                        lhsm(1,i2,j)*rhs(m,i,j,k)
             end do
          end do

c---------------------------------------------------------------------
c         And again the last two rows separately
c---------------------------------------------------------------------
          do  j = 1, ny-2
             i  = nx-2
             i1 = nx-1
             m = 4
             fac1       = 1.d0/lhsp(3,i,j)
             lhsp(4,i,j)  = fac1*lhsp(4,i,j)
             lhsp(5,i,j)  = fac1*lhsp(5,i,j)
             rhs(m,i,j,k) = fac1*rhs(m,i,j,k)
             lhsp(3,i1,j) = lhsp(3,i1,j) -
     >                      lhsp(2,i1,j)*lhsp(4,i,j)
             lhsp(4,i1,j) = lhsp(4,i1,j) -
     >                      lhsp(2,i1,j)*lhsp(5,i,j)
             rhs(m,i1,j,k) = rhs(m,i1,j,k) -
     >                      lhsp(2,i1,j)*rhs(m,i,j,k)
             m = 5
             fac1       = 1.d0/lhsm(3,i,j)
             lhsm(4,i,j)  = fac1*lhsm(4,i,j)
             lhsm(5,i,j)  = fac1*lhsm(5,i,j)
             rhs(m,i,j,k) = fac1*rhs(m,i,j,k)
             lhsm(3,i1,j) = lhsm(3,i1,j) -
     >                      lhsm(2,i1,j)*lhsm(4,i,j)
             lhsm(4,i1,j) = lhsm(4,i1,j) -
     >                      lhsm(2,i1,j)*lhsm(5,i,j)
             rhs(m,i1,j,k) = rhs(m,i1,j,k) -
     >                      lhsm(2,i1,j)*rhs(m,i,j,k)
c---------------------------------------------------------------------
c               Scale the last row immediately
c---------------------------------------------------------------------
             rhs(4,i1,j,k) = rhs(4,i1,j,k)/lhsp(3,i1,j)
             rhs(5,i1,j,k) = rhs(5,i1,j,k)/lhsm(3,i1,j)
          end do

c---------------------------------------------------------------------
c                         BACKSUBSTITUTION
c---------------------------------------------------------------------

          do  j = 1, ny-2
             i  = nx-2
             i1 = nx-1
             do   m = 1, 3
                rhs(m,i,j,k) = rhs(m,i,j,k) -
     >                             lhs(4,i,j)*rhs(m,i1,j,k)
             end do

             rhs(4,i,j,k) = rhs(4,i,j,k) -
     >                          lhsp(4,i,j)*rhs(4,i1,j,k)
             rhs(5,i,j,k) = rhs(5,i,j,k) -
     >                          lhsm(4,i,j)*rhs(5,i1,j,k)
          end do

c---------------------------------------------------------------------
c      The first three factors
c---------------------------------------------------------------------
          do  j = 1, ny-2
             do    i = nx-3, 0, -1
                i1 = i  + 1
                i2 = i  + 2
                do   m = 1, 3
                   rhs(m,i,j,k) = rhs(m,i,j,k) -
     >                          lhs(4,i,j)*rhs(m,i1,j,k) -
     >                          lhs(5,i,j)*rhs(m,i2,j,k)
                end do

c---------------------------------------------------------------------
c      And the remaining two
c---------------------------------------------------------------------
                rhs(4,i,j,k) = rhs(4,i,j,k) -
     >                          lhsp(4,i,j)*rhs(4,i1,j,k) -
     >                          lhsp(5,i,j)*rhs(4,i2,j,k)
                rhs(5,i,j,k) = rhs(5,i,j,k) -
     >                          lhsm(4,i,j)*rhs(5,i1,j,k) -
     >                          lhsm(5,i,j)*rhs(5,i2,j,k)
             end do
          end do

       end do
!$OMP END PARALLEL DO
       if (timeron) call timer_stop(t_xsolve)

c---------------------------------------------------------------------
c      Do the block-diagonal inversion
c---------------------------------------------------------------------
       if (timeron) call timer_start(t_ninvr)
       call ninvr(rhs, nx, nxmax, ny, nz)
       if (timeron) call timer_stop(t_ninvr)

       return
       end

max_stars_repo_name: robert-mijakovic/readex-ptf

max_issues_repo_name: robert-mijakovic/readex-ptf

max_forks_repo_name: robert-mijakovic/readex-ptf


text:
C     LAST UPDATE 16/10/96
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      PROGRAM SAMPLE
      IMPLICIT NONE
C
C Purpose: Re-sample CCP13  intensity output data by fitting either to
C          the Shannon sampling points or (more efficiently) to the
C          Fourier-Bessel sampling positions.
C
C Calls   2: RDCOM  , ROOTJ
C Called by:
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C Parameters:
C
       INTEGER NLLPX,MAXRTS,MAXLPT
       PARAMETER(NLLPX=4096,MAXRTS=128,MAXLPT=512)
       REAL PI,TOL
       PARAMETER(PI=3.141592654,TOL=1.0E-05)
C
C Local arrays:
C
       REAL VALS(10),RINT(NLLPX),SIG(NLLPX),RL(2,NLLPX),D(NLLPX),
     &      ROOTS(MAXRTS),Y(MAXLPT),FB(MAXLPT,MAXRTS),V(MAXRTS,MAXRTS),
     &      W(MAXRTS),CVM(MAXRTS,MAXRTS),X(MAXRTS)
       INTEGER ITEM(20),INDEX(MAXLPT)
       CHARACTER*40 WORD(10)
C
C Local variables:
C
       REAL RESMIN,RESLIM,DMAX,DR,SIGO,RMIN,RMAX,CELL,A,SIGTMP,
     &      TMP,CHISQ,WMAX,THRESH,XMAX,XSCALE
       INTEGER NW,NV,IRC,IREF,I,J,M,L,LMIN,LMAX,NDATA,NSVREJ,NSTART,NMAX
       INTEGER KUNIT,N,K,NR,UHELIX,THELIX,LACTUAL,LSCALE,NSTACK
       LOGICAL NEXT,FEX
       CHARACTER*40 OFNAM,INFIL
C
C External functions:
C
       INTEGER NMIN
       REAL KCOF
       EXTERNAL NMIN,KCOF
C
C Data:
C
       DATA KUNIT /12/
C-----------------------------------------------------------------------
C
C========Open log file
C
      OPEN(UNIT=4,FILE='SAMPLE.LOG',STATUS='UNKNOWN')
      WRITE(6,1000)
      WRITE(4,1000)
C
C========Set up defaults
C
      NEXT = .TRUE.
      RESLIM = 0.0
      RESMIN = 0.0
      DMAX = 0.0
      SIGO = 0.0
      LMIN = 0
      LMAX = 0
      UHELIX = 1
      THELIX = 1
      NSTART = 1
      NSTACK = 1
      LSCALE = 1
      XSCALE = 5.0
      A = 100.0
C
C=======Read keyworded input
C
 10   WRITE(6,'(A8,$)')'SAMPLE> '
      CALL RDCOM(5,6,WORD,NW,VALS,NV,ITEM,4,IRC)
      IF(IRC.EQ.2)GOTO 10
      IF(IRC.EQ.1)GOTO 9999
      IF(ITEM(1).EQ.0)GOTO 10
      IF(ITEM(1).NE.2)THEN
         WRITE(6,1010)
         WRITE(4,1010)
         GOTO 10
      ENDIF
C
      IF(WORD(1)(1:3).EQ.'RUN')THEN
         NEXT = .FALSE.
      ELSEIF(WORD(1)(1:4).EQ.'FILE')THEN
C
C========Read input filename
C
         WRITE(6,1030)
         WRITE(4,1030)
         READ(5,'(A)',END=9999)INFIL
         WRITE(4,'(A)')INFIL
         INQUIRE(FILE=INFIL,EXIST=FEX)
         IF(.NOT.FEX)THEN
            WRITE(6,1036)INFIL
            WRITE(4,1036)INFIL
            GOTO 10
         ENDIF
         WRITE(6,1035)
         WRITE(4,1035)
         READ(5,'(A)',END=9999)OFNAM
         WRITE(4,'(A)')OFNAM
         OPEN(UNIT=KUNIT,FILE=OFNAM,STATUS='UNKNOWN',IOSTAT=IRC)
         IF(IRC.NE.0)THEN
            WRITE(6,1046)OFNAM
            WRITE(4,1046)OFNAM
            GOTO 9999
         ENDIF
      ELSEIF(WORD(1)(1:4).EQ.'RESO')THEN
         IF(ITEM(2).NE.1)THEN
            WRITE(6,1020)
            WRITE(4,1020)
         ELSE
            RESMIN = VALS(1)
            RESLIM = VALS(2)
         ENDIF
      ELSEIF(WORD(1)(1:4).EQ.'SIGM')THEN
         IF(ITEM(2).NE.1)THEN
            WRITE(6,1020)
            WRITE(4,1020)
         ELSE
            SIGO = VALS(1)
         ENDIF
      ELSEIF(WORD(1)(1:4).EQ.'EXTR')THEN
         IF(ITEM(2).NE.1)THEN
            WRITE(6,1020)
            WRITE(4,1020)
         ELSE
            XSCALE = VALS(1)
         ENDIF
      ELSEIF(WORD(1)(1:4).EQ.'HELI')THEN
         IF(ITEM(2).NE.1)THEN
            WRITE(6,1020)
            WRITE(4,1020)
         ELSE
            IF(NV.GE.1)UHELIX = NINT(VALS(1))
            IF(NV.GE.2)THELIX = NINT(VALS(2))
         ENDIF
      ELSEIF(WORD(1)(1:4).EQ.'STAR')THEN
         IF(ITEM(2).NE.1)THEN
            WRITE(6,1020)
            WRITE(4,1020)
         ELSE
            NSTART = NINT(VALS(1))
         ENDIF
      ELSEIF(WORD(1)(1:4).EQ.'STAC')THEN
         IF(ITEM(2).NE.1)THEN
            WRITE(6,1020)
            WRITE(4,1020)
         ELSE
            NSTACK = NINT(VALS(1))
         ENDIF
      ELSEIF(WORD(1)(1:4).EQ.'RADI')THEN
         IF(ITEM(2).NE.1)THEN
            WRITE(6,1020)
            WRITE(4,1020)
         ELSE
            A = VALS(1)
         ENDIF
      ELSEIF(WORD(1)(1:4).EQ.'SCAL')THEN
         IF(ITEM(2).NE.1)THEN
            WRITE(6,1020)
            WRITE(4,1020)
         ELSE
            LSCALE = NINT(VALS(1))
         ENDIF
      ELSE
         WRITE(6,1025)
         WRITE(4,1025)
      ENDIF
C
      IF(NEXT)GOTO 10
C
C========Open input file to determine Bragg/continuous etc
C
      IREF = 0
      OPEN(UNIT=10,FILE=INFIL,STATUS='OLD',IOSTAT=IRC)
      IF(IRC.NE.0)THEN
         WRITE(6,1045)INFIL
         WRITE(4,1045)INFIL
         GOTO 9999
      ENDIF
 30   CALL RDCOM(10,0,WORD,NW,VALS,NV,ITEM,0,IRC)
      IF(IRC.EQ.2)THEN
         WRITE(6,1050)
         WRITE(4,1050)
      ELSEIF(IRC.EQ.1)THEN
         GOTO 45
      ENDIF
      IF(ITEM(1).EQ.2)THEN
         IF(WORD(1)(1:4).EQ.'BRAG')THEN
            WRITE(6,2000)
            WRITE(4,2000)
            STOP
         ELSEIF(WORD(1)(1:4).EQ.'CELL')THEN
            CELL = VALS(1)
            WRITE(KUNIT,1047)CELL
         ELSEIF(WORD(1)(1:4).EQ.'DELT')THEN
            DR = VALS(1)
         ENDIF
      ELSEIF(ITEM(1).EQ.1)THEN
         IREF = IREF + 1
         IF(IREF.GT.NLLPX)THEN
            WRITE(6,1055)NLLPX
            WRITE(4,1055)NLLPX
            GOTO 45
         ENDIF
         RL(2,IREF) = VALS(3)
         RL(1,IREF) = VALS(4)
         RINT(IREF) = VALS(6)
         SIG(IREF) = VALS(7)
         D(IREF) = SQRT((RL(2,IREF)/CELL)**2
     &                  +RL(1,IREF)**2)
         IF(D(IREF).GT.DMAX)DMAX = D(IREF)
         IF(NINT(RL(2,IREF)).LT.LMIN)LMIN = NINT(RL(2,IREF))
         IF(NINT(RL(2,IREF)).GT.LMAX)LMAX = NINT(RL(2,IREF))
      ENDIF
      GOTO 30
 45   WRITE(6,1060)IREF
      WRITE(4,1060)IREF
      CLOSE(10)
      WRITE(6,1070)DMAX
      WRITE(4,1070)DMAX
      IF(RESLIM.EQ.0.0)RESLIM = DMAX
      DO 140 L=LMIN,LMAX
         J = 0
         RMIN = 1.0
         RMAX = -1.0
         DO 50 I=1,IREF
            IF(NINT(RL(2,I)).EQ.L)THEN
               IF(SIG(I).GT.0.0.AND.RINT(I).GT.0.0)THEN
                  IF(RINT(I)/SIG(I).GT.SIGO.AND.
     &               D(I).LT.RESLIM.AND.D(I).GT.RESMIN)THEN
                     J = J + 1
                     IF(RL(1,I).LT.RMIN)RMIN = RL(1,I)
                     IF(RL(1,I).GT.RMAX)RMAX = RL(1,I)
                     INDEX(J) = I
                  ENDIF
               ENDIF
            ENDIF
 50      CONTINUE
         NMAX = NINT(2.0*PI*A*RMAX+2.0)
         WRITE(6,1080)L,J,RMIN,RMAX,NMAX
         WRITE(4,1080)L,J,RMIN,RMAX,NMAX
C
C========Find minimum bessel function order contributing to this
C========layerline
C
         LACTUAL = LSCALE*L
         N = NMIN(LACTUAL,UHELIX,THELIX,NSTART,NSTACK,NMAX)
         IF(N.LT.0)THEN
            WRITE(6,1085)NMAX
            WRITE(4,1085)NMAX
            GOTO 140
         ENDIF
C
C========Determine the largest root required and find roots
C
         XMAX = PI*(XSCALE+4.0*A*RMAX)
         CALL ROOTJ(2*N,XMAX,ROOTS,MAXRTS,NR)
         WRITE(6,1090)N,NR
         WRITE(4,1090)N,NR
C
C========Warn user about under-sampling
C
         IF(J.LT.NR)THEN
            WRITE(6,2010)
            WRITE(4,2010)
         ENDIF
C
C========Form design matrix for Fourier-Bessel smoothing
C
         NDATA = MAX(NR,J)
         DO 65 K=1,NDATA
            IF(K.LE.J)THEN
               M = INDEX(K)
               TMP = RL(1,M)*4.0*PI*A
               SIGTMP = SIG(M)
               DO 55 I=1,NR
                  FB(K,I) = KCOF(TMP,ROOTS(I),N)/SIGTMP
 55            CONTINUE
               Y(K) = RINT(M)/SIGTMP
            ELSE
               DO 60 I=1,NR
                  FB(K,I) = 0.0
 60            CONTINUE
               Y(K) = 0.0
            ENDIF
 65      CONTINUE
C
C========Decompose design matrix
C
         CALL SVDCMP(FB,NDATA,NR,MAXLPT,MAXRTS,W,V)
C
C========Filter singular values
C
         WMAX = 0.0
         DO 70 I=1,NR
            IF(W(I).GT.WMAX)WMAX = W(J)
 70      CONTINUE
         THRESH = TOL*WMAX
         NSVREJ = 0
         DO 80 I=1,NR
            IF(W(I).LT.THRESH)THEN
               W(I) = 0.0
               NSVREJ = NSVREJ + 1
            ENDIF
 80      CONTINUE
         WRITE(6,1100)NSVREJ
         WRITE(4,1100)NSVREJ
C
C========Back-substitute to get least-squares solution
C
         CALL SVBKSB(FB,W,V,NDATA,NR,MAXLPT,MAXRTS,Y,X)
C
C========Form covariance matrix
C
         CALL SVDVAR(V,NR,MAXRTS,W,CVM,MAXRTS)
C
C========Write out reduced data
C
         DO 90 I=1,NR
            TMP = ROOTS(I)/(4.0*PI*A)
            IF(TMP.GE.RMIN.AND.TMP.LE.RMAX)
     &         WRITE(KUNIT,2000,ERR=9999)L,TMP,X(I),SQRT(CVM(I,I))
 90      CONTINUE
C
C========Compute chi-squared
C
         CHISQ = 0.0
         DO 130 K=1,NDATA
            IF(K.LE.J)THEN
               M = INDEX(K)
               TMP = RL(1,M)*4.0*PI*A
               SIGTMP = SIG(M)
               DO 100 I=1,NR
                  FB(K,I) = KCOF(TMP,ROOTS(I),N)/SIGTMP
 100           CONTINUE
            ELSE
               DO 110 I=1,NR
                  FB(K,I) = 0.0
 110           CONTINUE
            ENDIF
            TMP = 0.0
            DO 120 I=1,NR
               TMP = TMP + FB(K,I)*X(I)
 120        CONTINUE
            TMP = TMP - Y(K)
            CHISQ = CHISQ + TMP*TMP
 130     CONTINUE
         WRITE(6,1110)CHISQ
         WRITE(4,1110)CHISQ
 140  CONTINUE
      CLOSE(KUNIT)
      STOP 'Normal stop'
 9999 STOP 'USER STOP'
C
 1000 FORMAT(1X,'SAMPLE: Re-sampling of continuous layer-line',
     &       1X,'intensity'/
     &       1X,'Last update 16/10/96'/)
 1010 FORMAT(1X,'***Keyword expected')
 1020 FORMAT(1X,'***Value expected')
 1025 FORMAT(1X,'***Unrecognised keyword')
 1030 FORMAT(1X,'Enter input filename')
 1035 FORMAT(1X,'Enter output filename')
 1036 FORMAT(1X,'***Input file does not exist:  ',A40)
 1045 FORMAT(1X,'***Error opening input file:   ',A40)
 1046 FORMAT(1X,'***Error opening output file:  ',A40)
 1047 FORMAT(1X,'CONTINUOUS'/1X,'CELL ',F10.4)
 1050 FORMAT(1X,'SAMPLE: Error reading input file')
 1055 FORMAT(1X,'Number of reflections exceeds ',I4)
 1060 FORMAT(1X,'Number of reflections read from file:',I5)
 1070 FORMAT(1X,'Maximum d*  ',F8.6)
 1080 FORMAT(/1X,'Layer line ',I3,5X,'No. of intensities read ', I4/
     &       1X,'Min. rec. radius ',F8.6,5X,'Max. rec. radius ',F8.6/
     &       1X,'Maximum possible |n| ',I4)
 1085 FORMAT(1X,'Warning: No orders contribute for |n| <=',I3/
     &       1X,'=== Skipping layerline ===')
 1090 FORMAT(1X,'Minimum |n| ',I3/
     &       1X,'No. of roots of 2|n| determined ',I4)
 1100 FORMAT(1X,'No. of singular values filtered ',I4)
 1110 FORMAT(1X,'Chi-squared for this layerline = ',E12.4)
 2000 FORMAT(2(1X,'999'),I4,E12.4,3X,'1',2E12.4)
 2010 FORMAT(1X,'Warning: No. of roots > No. of intensities')
 1140 FORMAT(1X,'Number of reflections written  ',I8)
      END


! { dg-do compile }
subroutine foo ()
   real(4), allocatable, save :: tmp (:, :)
   real(4), pointer, save :: arr (:, :, :)
   integer :: l, m, n
   tmp = (cshift(cshift(arr (:,:,l),m,2),n,1))
end subroutine foo