text:
!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
!> \brief  This module provides functions to deal with time and increment of stress analysis

module m_fstr_TimeInc
  use m_fstr
  use m_step
  use m_timepoint

  implicit none

  real(kind=kreal), private :: current_time = 0.d0 !< current time
  real(kind=kreal), private :: time_inc = 0.d0 !< current increment of time
  real(kind=kreal), private :: time_inc_base = 0.d0 !< base increment of time

contains

  !! basic functions to get time/timeinc
  real(kind=kreal) function fstr_get_time()
    fstr_get_time = current_time
  end function

  real(kind=kreal) function fstr_get_timeinc()
    fstr_get_timeinc = time_inc
  end function

  real(kind=kreal) function fstr_get_timeinc_base()
    fstr_get_timeinc_base = time_inc_base
  end function

  !! basic functions to set/update time/timeinc
  subroutine fstr_set_time( time )
    real(kind=kreal), intent(in) :: time  !< input time
    current_time = time
  end subroutine

  subroutine fstr_set_timeinc( dtime )
    real(kind=kreal), intent(in) :: dtime  !< input increment of time
    time_inc = dtime
  end subroutine

  subroutine fstr_set_timeinc_base( dtime_base )
    real(kind=kreal), intent(in) :: dtime_base  !< input increment of time
    time_inc_base = dtime_base
  end subroutine

  subroutine fstr_proceed_time()
    current_time = current_time + time_inc
  end subroutine

  !! functions to print status
  subroutine fstr_TimeInc_PrintSTATUS_init
    write(ISTA,'(A10,"-+-",A60,"-+-",A40)') repeat("-",10),repeat("-",60),repeat("-",40)
    write(ISTA,'(2A5," | ",2A5,2A7,3A12," | ",A40)') '     ', '     ', ' ', ' # of',  'MAX #', 'TOT #',  '','', '', ''
    write(ISTA,'(2A5," | ",2A5,2A7,3A12," | ",A7,A33)') &
      & 'STEP', 'SUB', 'STAT', ' CONT', 'NEWTON', 'NEWTON', 'START', 'TIME','END', 'MESSAGE',''
    write(ISTA,'(2A5," | ",2A5,2A7,3A12," | ",A40)') &
      & '', 'STEP', '  ', 'ITER',  'ITER', 'ITER',  'TIME', 'INC', 'TIME',''
    write(ISTA,'(A10,"-+-",A60,"-+-",A40)') repeat("-",10),repeat("-",60),repeat("-",40)
  end subroutine

  subroutine fstr_TimeInc_PrintSTATUS( stepinfo, fstrPARAM, totstep, substep, NRstatI, NRstatR, AutoINC_stat, Cutback_stat )
    type(step_info), intent(in)       ::  stepinfo
    type(fstr_param), intent(in)      ::  fstrPARAM
    integer(kind=kint), intent(in)    ::  totstep
    integer(kind=kint), intent(in)    ::  substep
    integer(kind=kint), intent(in)    ::  NRstatI(:)  !previous Newton Raphson Iteration Log
    real(kind=kreal), intent(in)      ::  NRstatR(:)  !previous Newton Raphson Iteration Log
    integer(kind=kint), intent(inout) ::  AutoINC_stat
    integer(kind=kint), intent(in)    ::  Cutback_stat

    character(len=5) :: cstep, csstep, cstate, ccont
    character(len=7) :: cmaxn, ctotn
    character(len=12) :: ctime, cdtime, etime
    character(len=40) :: message
    type(tParamAutoInc) :: pAinc

    pAinc = fstrPARAM%ainc(stepinfo%AincParam_id)

    write(cstep,'(I5)') totstep
    write(csstep,'(I5)') substep
    write(ccont,'(I5)') NRstatI(knstCITER)
    write(cmaxn,'(I7)') NRstatI(knstMAXIT)
    write(ctotn,'(I7)') NRstatI(knstSUMIT)
    write(ctime,'(1pE12.4)') current_time
    write(cdtime,'(1pE12.4)') time_inc

    if( Cutback_stat > 0 ) then
      write(etime,'(1pE12.4)') current_time
      write(cstate,'(I4,A)') Cutback_stat,'F'
      if( NRstatI(knstDRESN) == 1 ) write(message,'(A)') 'Failed to converge due to MAXITER.'
      if( NRstatI(knstDRESN) == 2 ) write(message,'(A)') 'Failed to converge due to MAXRES.'
      if( NRstatI(knstDRESN) == 3 ) write(message,'(A)') 'Failed to converge due to MAXCONTITER.'
      if( NRstatI(knstDRESN) == 4 ) write(message,'(A)') 'Failed to converge due to MatSolveError.'
      if( Cutback_stat == pAinc%CBbound ) write(message,'(A)') '# of successive cutback reached max.'
    else
      write(etime,'(1pE12.4)') current_time+time_inc
      write(cstate,'(A5)') 'S'
      write(message,'(A)') ''
    endif

    write(ISTA,'(2A5," | ",2A5,2A7,3A12," | ",A40)') cstep, csstep, cstate, ccont, &
      &  cmaxn, ctotn, ctime, cdtime, etime, message

  end subroutine

  subroutine fstr_TimeInc_PrintSTATUS_final(success_flag)
    logical, intent(in) :: success_flag
    write(ISTA,'(A10,"-+-",A60,"-+-",A40)') repeat("-",10),repeat("-",60),repeat("-",40)
    if(success_flag) then
      write(ISTA,'(A)') 'FSTR_SOLVE_NLGEOM HAS COMPLETED SUCCESSFULLY'
    else
      write(ISTA,'(A)') 'FSTR_SOLVE_NLGEOM HAS NOT COMPLETED SUCCESSFULLY'
    end if
  end subroutine

  !! true if step finished
  logical function fstr_TimeInc_isStepFinished( stepinfo )
    type(step_info), intent(in)  :: stepinfo

    real(kind=kreal) :: endtime, remain_time

    endtime = stepinfo%starttime + stepinfo%elapsetime
    remain_time = (endtime-current_time)/stepinfo%elapsetime
    fstr_TimeInc_isStepFinished = (remain_time < 1.d-12)
  end function

  !! true if current time is time point
  logical function fstr_TimeInc_isTimePoint( stepinfo, fstrPARAM )
    type(step_info), intent(in)  :: stepinfo
    type(fstr_param), intent(in) :: fstrPARAM

    fstr_TimeInc_isTimePoint = .false.
    if( stepinfo%inc_type == stepFixedInc ) return
    if( stepinfo%timepoint_id == 0 ) return

    fstr_TimeInc_isTimePoint = is_at_timepoints(current_time,stepinfo%starttime, &
      &  fstrPARAM%timepoints(stepinfo%timepoint_id))
  end function

  !! set time_inc from stepinfo and Newton Raphson iteration status
  subroutine fstr_TimeInc_SetTimeIncrement( stepinfo, fstrPARAM, substep, NRstatI, NRstatR, AutoINC_stat, Cutback_stat )
    type(step_info), intent(in)       ::  stepinfo
    type(fstr_param), intent(in)      ::  fstrPARAM
    integer(kind=kint), intent(in)    ::  substep
    integer(kind=kint), intent(in)    ::  NRstatI(:)  !previous Newton Raphson Iteration Log
    real(kind=kreal), intent(in)      ::  NRstatR(:)  !previous Newton Raphson Iteration Log
    integer(kind=kint), intent(inout) ::  AutoINC_stat
    integer(kind=kint), intent(in)    ::  Cutback_stat

    real(kind=kreal) :: timeinc0
    real(kind=kreal) :: endtime, remain
    logical          :: to_be_decreased, to_be_increased
    type(tParamAutoInc) :: pAinc

    pAinc = fstrPARAM%ainc(stepinfo%AincParam_id)

    if( stepinfo%inc_type == stepFixedInc ) then
      timeinc0 = stepinfo%initdt
    else ! INCTYPE==AUTO
      if( Cutback_stat > 0 ) then
        timeinc0 = pAinc%ainc_Rc*time_inc
        if(myrank == 0) write(*,'(2(A,E10.3))') 'time increment is decreased from ', time_inc_base, ' to ', timeinc0
        AutoINC_stat = -1
      else
        if( substep == 1 ) then
          timeinc0 = stepinfo%initdt
          AutoINC_stat = 0
        else
          timeinc0 = time_inc_base

          !decrease condition
          to_be_decreased = .false.
          if( NRstatI(knstMAXIT) > pAinc%NRbound_s(knstMAXIT) ) to_be_decreased = .true.
          if( NRstatI(knstSUMIT) > pAinc%NRbound_s(knstSUMIT) ) to_be_decreased = .true.
          if( NRstatI(knstCITER) > pAinc%NRbound_s(knstCITER) ) to_be_decreased = .true.

          !increase condition
          to_be_increased = .true.
          if( NRstatI(knstMAXIT) > pAinc%NRbound_l(knstMAXIT) ) to_be_increased = .false.
          if( NRstatI(knstSUMIT) > pAinc%NRbound_l(knstSUMIT) ) to_be_increased = .false.
          if( NRstatI(knstCITER) > pAinc%NRbound_l(knstCITER) ) to_be_increased = .false.

          ! count # of times that increase/decrease condition has been satisfied
          ! AutoINC_stat < 0 ... decrease condition has been satisfied -AutoINC_stat times
          ! AutoINC_stat > 0 ... increase condition has been satisfied AutoINC_stat times
          if(to_be_decreased) then
            AutoINC_stat = min(AutoINC_stat,0) - 1
          else if(to_be_increased) then
            AutoINC_stat = max(AutoINC_stat,0) + 1
          else
            AutoINC_stat = 0
          end if

          if( AutoINC_stat <= -pAinc%NRtimes_s ) then
            timeinc0 = pAinc%ainc_Rs*timeinc0
            if(myrank == 0) write(*,'(2(A,E10.3))') 'time increment is decreased from ', time_inc_base, ' to ', timeinc0
          else if( AutoINC_stat >= pAinc%NRtimes_l ) then
            timeinc0 = dmin1(pAinc%ainc_Rl*timeinc0,stepinfo%maxdt)
            if(myrank == 0) write(*,'(2(A,E10.3))') 'time increment is increased from ', time_inc_base, ' to ', timeinc0
          end if
        end if
      end if
    end if

    if( timeinc0 < stepinfo%mindt ) then
      if(myrank == 0) then
        write(*,'(2(A,E10.3))') 'Error: current time increment ',timeinc0,' is smaller than lower bound',stepinfo%mindt
        call fstr_TimeInc_PrintSTATUS_final(.false.)
      endif
      stop
    end if

    endtime = stepinfo%starttime + stepinfo%elapsetime
    time_inc_base = dmin1(timeinc0,endtime-current_time)
    if( stepinfo%timepoint_id > 0 ) then
      remain = get_remain_to_next_timepoints(current_time,stepinfo%starttime,fstrPARAM%timepoints(stepinfo%timepoint_id))
      time_inc = dmin1(time_inc_base,remain)
    else
      time_inc = time_inc_base
    end if

  end subroutine

end module m_fstr_TimeInc

max_stars_repo_name: masae-hayashi/FrontISTR

max_issues_repo_name: masae-hayashi/FrontISTR

max_forks_repo_name: masae-hayashi/FrontISTR


text:
module alt_micro

contains
    subroutine alt_Qimp(id,ierr)
        use constants_def
        use NScool_def
        integer, intent(in) :: id
        integer, intent(out) :: ierr
        type(NScool_info), pointer :: s
        real(dp) :: rhoC

        call get_NScool_info_ptr(id,s,ierr)
        if (ierr /= 0) return
        rhoC = s% extra_real_controls(1)

        where(s% rho_bar > rhoC)
            s% ionic_bar% Q = s% extra_real_controls(2)
        end where
    end subroutine alt_Qimp

    subroutine alt_sf(id,kp,kn,Tc)
        use constants_def
        use superfluid_def
        use superfluid_lib
        use NScool_def

        integer, intent(in) :: id   ! for accessing parameters
        real(dp), intent(in) :: kp, kn  ! proton, neutron wavevectors (fm**-1)
        real(dp), dimension(max_number_sf_types), intent(out) :: Tc ! (K)
        type(NScool_info), pointer :: s
        real(dp) :: rhoC, kC, n
        integer :: ierr

        Tc = 1.0e10
        call get_NScool_info_ptr(id,s,ierr)
        if (ierr /= 0) return

        rhoC = s% extra_real_controls(1)
        n = rhoC*avogadro
        kC = (0.5*threepisquare*n)**onethird / cm_to_fm
        if (kn > kC) then
            Tc(neutron_1S0) = 0.0
        end if
    end subroutine alt_sf

end module alt_micro

ext: f

max_stars_repo_name: LBJ-Wade/dStar

max_issues_repo_name: LBJ-Wade/dStar

max_forks_repo_name: LBJ-Wade/dStar


text:
      program equiv11

C     Check that an external cannot appear in an EQUIVALENCE: PIPS does not

      real y(100)
      equivalence (x, z), (u, v)
      real u(100)
      external v

      do i = 2, 10
         y(i) = 0.
      enddo

      print *, i, j

      end


text:
subroutine doMove(NLINE, lSize, nTube, i, j, move, nmove, matrix)
  implicit none
  integer :: NLINE, i, j, nmove, nTube, lSize, movingIndex
  integer, dimension(nTube, NLINE) :: matrix
  integer, dimension(lSize, 2) :: move

  integer :: moving, nempty, ntomove
  movingIndex = 1 ! matrix(i,movingIndex) is the moving color
  do while(matrix(i,movingIndex) == 0)
     movingIndex = movingIndex + 1
  end do
  moving = matrix(i,movingIndex)

  ! do the movement
  nmove = nmove + 1
  move(nmove, 1) = i
  move(nmove, 2) = j

  ! calculate new colors
  nempty = count(matrix(j,:) == 0)
  ntomove = 1
  if (movingIndex /= NLINE) then
     do while(moving == matrix(i, movingIndex + ntomove))
        ntomove = ntomove + 1
     end do
  end if

  if (ntomove > nempty) then
     ntomove = nempty
     matrix(i, movingIndex : (movingIndex + ntomove - 1) ) = 0 ! becomes empty
     matrix(j, 1:nempty) = moving ! gets `moving` color
  else
     matrix(i, movingIndex : (movingIndex + ntomove - 1) ) = 0 ! becomes empty
     matrix(j, (nempty - ntomove + 1):nempty) = moving
  end if
end subroutine doMove

max_stars_repo_name: masoudshao/watersortpuzzle

max_issues_repo_name: masoudshao/watersortpuzzle

max_forks_repo_name: masoudshao/watersortpuzzle


text:
      subroutine option (icas,int,ital,ior)
      call optn (4hcase,icas)
      if (int .eq. 0) call optn (4hintn,3hlow)
      if (int .eq. 1) call optn (4hintn,4hhigh)
      call optn (4hfont,ital)
      call optn (4horen,ior)
      return
      end

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
module find_fit_module

!! This module contains a general function find_fit() for a nonlinear least
!! squares fitting. The function can fit any nonlinear expression to any data.

use minpack_module, only: wp, lmdif1

implicit none

private

abstract interface
    function expr_f(x, pars) result(y)
        import :: wp
        implicit none
        real(wp), intent(in) :: x(:)
        real(wp), intent(in) :: pars(:)
        real(wp) :: y(size(x))
    end function expr_f
end interface

public :: find_fit

contains

subroutine find_fit(data_x, data_y, expr, pars)
! Fits the (data_x, data_y) arrays with the function expr(x, pars).
! The user can provide any nonlinear function 'expr' depending on any number of
! parameters 'pars' and it must return the evaluated expression on the
! array 'x'. The arrays 'data_x' and 'data_y' must have the same
! length.
real(wp), intent(in) :: data_x(:)
real(wp), intent(in) :: data_y(:)
procedure(expr_f) :: expr
real(wp), intent(inout) :: pars(:)

real(wp) :: tol, fvec(size(data_x))
integer :: iwa(size(pars)), info, m, n
real(wp), allocatable :: wa(:)

tol = sqrt(epsilon(1.0_wp))
m = size(fvec)
n = size(pars)
allocate(wa(m*n + 5*n + m))
call lmdif1(fcn, m, n, pars, fvec, tol, info, iwa, wa, size(wa))
if (info /= 1) error stop "failed to converge"

contains

    subroutine fcn(m, n, x, fvec, iflag)
    integer, intent(in) :: m
    integer, intent(in) :: n
    integer, intent(inout) :: iflag
    real(wp), intent(in) :: x(n)
    real(wp), intent(out) :: fvec(m)
    fvec(1) = iflag ! Suppress compiler warning
    fvec = data_y - expr(data_x, x)
    end subroutine fcn

end subroutine find_fit

end module find_fit_module

program example_primes

!! Find a nonlinear fit of the form `a*x*log(b + c*x)` to a list of primes.

use find_fit_module, only: find_fit
use minpack_module, only: wp
use iso_fortran_env, only: nwrite => output_unit

implicit none

real(wp) :: pars(3)
real(wp), parameter :: y(*) = real([2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, &
                                    37, 41, 43, 47, 53, 59, 61, 67, 71], wp)
integer :: i
pars = [1.0_wp, 1.0_wp, 1.0_wp]
call find_fit([(real(i, wp), i=1,size(y))], y, expression, pars)

write(nwrite, '(1p,5x,a/(5x,3e15.7))') 'pars: ',pars

contains

function expression(x, pars) result(y)
real(wp), intent(in) :: x(:)
real(wp), intent(in) :: pars(:)
real(wp) :: y(size(x))
real(wp) :: a, b, c
a = pars(1)
b = pars(2)
c = pars(3)
y = a*x*log(b + c*x)
end function expression

end program example_primes

max_stars_repo_name: awvwgk/minpack

max_issues_repo_name: fortran-lang/minpack

max_forks_repo_name: awvwgk/minpack


text:
module HP_LaxFriedrichs
#include <messenger.h>
    use mod_kinds,              only: rk,ik
    use mod_constants,          only: ZERO,ONE,TWO,HALF,THREE
    use type_operator,          only: operator_t
    use type_chidg_worker,      only: chidg_worker_t
    use type_properties,        only: properties_t
    use DNAD_D
    implicit none

    !>
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/26/2019
    !!
    !-----------------------------------------------------------------------------------------------
    type, extends(operator_t), public :: HP_LaxFriedrichs_t

    contains

        procedure   :: init
        procedure   :: compute

    end type HP_LaxFriedrichs_t
    !***********************************************************************************************

    real(rk) :: p_param = 2._rk

contains

    !>
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   8/29/2016
    !!
    !-----------------------------------------------------------------------------------------------
    subroutine init(self)
        class(HP_LaxFriedrichs_t),   intent(inout)  :: self

        ! Set operator name
        call self%set_name('Hyperbolized Poisson LaxFriedrichs Operator')

        ! Set operator type
        call self%set_operator_type('Boundary Advective Operator')

        ! Set operator equations
        call self%add_primary_field('u')
        call self%add_primary_field('p')
        call self%add_primary_field('q')
        call self%add_primary_field('r')

    end subroutine init
    !***********************************************************************************************

    !>
    !!
    !!  @author Nathan A. Wukie
    !!
    !!
    !!
    !!
    !-----------------------------------------------------------------------------------------------
    subroutine compute(self,worker,prop)
        class(HP_LaxFriedrichs_t), intent(inout)   :: self
        type(chidg_worker_t),               intent(inout)   :: worker
        class(properties_t),                intent(inout)   :: prop

        type(AD_D), dimension(:), allocatable   ::  &
            u_m, u_p,                               &
            p_m, p_p,                               &
            q_m, q_p,                               &
            r_m, r_p,                               &
            dissipation_u, dissipation_p, dissipation_q, dissipation_r, &
            sumsqr_m, sumsqr_p, mag_m, mag_p,   &
            d_norm_m, d_norm_p, lambda_m, lambda_p

        real(rk)    :: alpha, beta

        !
        ! Interpolate solution to quadrature nodes
        !
        u_m = worker%get_field('u', 'value', 'face interior')
        u_p = worker%get_field('u', 'value', 'face exterior')

        p_m = worker%get_field('p', 'value', 'face interior')
        p_p = worker%get_field('p', 'value', 'face exterior')

        q_m = worker%get_field('q', 'value', 'face interior')
        q_p = worker%get_field('q', 'value', 'face exterior')

        r_m = worker%get_field('r', 'value', 'face interior')
        r_p = worker%get_field('r', 'value', 'face exterior')

        ! Compute normal for 1D system
        d_norm_m = p_m*worker%unit_normal(1) + q_m*worker%unit_normal(2) + r_m*worker%unit_normal(3)
        d_norm_p = p_p*worker%unit_normal(1) + q_p*worker%unit_normal(2) + r_p*worker%unit_normal(3)

        ! Eigenvalues for 1D p-Poisson system
        alpha = THREE + (p_param-TWO)/TWO
        beta  = TWO   + (p_param-TWO)/TWO

        lambda_m = sqrt(abs(alpha*d_norm_m**beta))
        lambda_p = sqrt(abs(alpha*d_norm_p**beta))

        ! OVERRIDE
        lambda_m = ONE
        lambda_p = ONE

        !
        ! Compute dissipation
        !
        !dissipation_u = -HALF*max(mag_m,mag_p)*(u_p-u_m)
        !dissipation_p = -HALF*max(mag_m,mag_p)*(p_p-p_m)
        !dissipation_q = -HALF*max(mag_m,mag_p)*(q_p-q_m)
        !dissipation_r = -HALF*max(mag_m,mag_p)*(r_p-r_m)

        dissipation_u = -HALF*max(lambda_m,lambda_p)*(u_p-u_m)
        dissipation_p = -HALF*max(lambda_m,lambda_p)*(p_p-p_m)
        dissipation_q = -HALF*max(lambda_m,lambda_p)*(q_p-q_m)
        dissipation_r = -HALF*max(lambda_m,lambda_p)*(r_p-r_m)

        !
        ! Integrate flux
        !
        call worker%integrate_boundary_upwind('u',dissipation_u)
        call worker%integrate_boundary_upwind('p',dissipation_p)
        call worker%integrate_boundary_upwind('q',dissipation_q)
        call worker%integrate_boundary_upwind('r',dissipation_r)

    end subroutine compute
    !*************************************************************************************************

end module HP_LaxFriedrichs

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
C=======================================================================
C  ML_PHENOL, Subroutine
C
C  Determines phenological stage
C-----------------------------------------------------------------------
C  Revision history
C
C  1. Written
C  2. Header revision and minor changes           P.W.W.      2-7-93
C  3. Added switch block, code cleanup            P.W.W.      2-7-93
C  4. Modified TT calculations to reduce line #'s P.W.W.      2-7-93
C  5. Modified for MILLET model                   W.T.B.      MAY 94
C  6. Converted to modular format                 W.D.B.      7-31-02
C  7. Major revisions to millet model             K.J.B April-May 2015
C-----------------------------------------------------------------------
C  INPUT  : YRDOY
C
C  LOCAL  : I,NDAS,L,L0,XANC,TEMPCN,TEMPCX,XS,TDIF,TCOR,
C           TTMP,YDL,DEC,DLV,CHGDL,SWSD,PDTT,RATEIN,PFLOWR,PSKER,
C           YIELDB,TWILEN,BARFAC,ABSTRES
C
C  OUTPUT : STGDOY
C-----------------------------------------------------------------------
C  Called : ML_CERES
C
C  Calls  : ML_PHASEI ML_COLD
C-----------------------------------------------------------------------
C
C                         DEFINITIONS
C
C  STGDOY :
C  YRDOY  :
C  NOUTDO :
C=======================================================================

      SUBROUTINE ML_PHENOL (APTNUP,BIOMAS,BIOMS2,
     & CSD1,
     & CNSD1, CNSD2, CTYPE, CUMDEP, DLAYR, DTT, SETTC, SETTYP,
     & EMAT, GNUP, GPP, GPSM, GRAINN, GRNWT, IDETO,
     & IDUR1, IPRINT, ISDATE, ISM, ISTAGE, ISWNIT, ISWWAT,
     & LAI, LEAFNO, LL, MAXLAI, MDATE, NOUTDO, NLAYR,
     & P1, P2O, P2R, P3, P4, P5, P9, PANWT, PGRNWT, PHINT,
     & PLTPOP, PWA, ROPT, RTDEP, SDEPTH, SI1, SI2, SI3,
     & SI4, SIND, SKERWT, SNOW, SRAD, STGDOY,STMWT, STOVER, STOVN,
     & SUMDTT, SW, TANC, TBASE, TEMPCR, TMAX, TMIN,
     & TOPT,TOTNUP, TPSM, YIELD, YRDOY,XGNP, XSTAGE,
C      Variables passed through PHENOL to phasei but not used in phenol
     & AGEFAC, BIOMS1, CUMDTT, CUMPH, G4, G5, GROLF,
     & GRORT, GROSTM, LFWT, MGROLF, MGROPAN, MGROSTM,
     & MLAG1, MLFWT, MPANWT, MPLAG, MPLA, MSTMWT, NSTRES, PAF,
     & PGC, PLA, PLAN, PLAO, PLATO, PLAMX, PTF, RANC,
     & RESERVE, RLV, ROOTN, RTWT, RWU, SEEDRV, SENLA, SLAN,
     & STOVWT, TCARBO, TGROLF,
     & TGROPAN, TGROSTM, TILN, TILSW, TLAG1, TLFWT, TLNO, TMNC,
     & TPANWT, TPLA, TPLAG, TSIZE, TSTMWT, VANC, VMNC, WSTR1,
     & XNTI,SWFAC,TURFAC,DGET,SWCG,DJTI,
     &    DAYL, TWILEN, CANWAA, CANNAA)

      USE MODULEDEFS

      IMPLICIT  NONE
      SAVE
C -----------------------------------------------------------------------
C VARIABLES ONLY USED IN PHASEI. THEY ARE PASSED THROUGH PHENOL TO PHASEI
C------------------------------------------------------------------------
      REAL AGEFAC
      REAL BIOMS1
      REAL CUMDTT
      REAL CUMPH
      REAL DGET
      REAL DJTI
      REAL G4
      REAL G5
      REAL GROLF
      REAL GRORT
      REAL GROSTM
      REAL LFWT
      REAL MGROLF
      REAL MGROPAN
      REAL MGROSTM
      REAL MLAG1
      REAL MLFWT
      REAL MPANWT
      REAL MPLAG
      REAL MPLA
      REAL MSTMWT
      REAL NSTRES
      REAL PAF
      REAL PGC
      REAL PLA
      REAL PLAN
      REAL PLAO
      REAL PLATO
      REAL PLAY
      REAL PLAMX
      REAL PTF
      REAL RANC
      REAL RESERVE
      REAL RLV(NL)
      REAL ROOTN
      REAL RTWT
      REAL RWU(NL)
      REAL SEEDRV
      REAL SENLA
      REAL SLAN
      REAL STOVWT
      REAL SWCG
      REAL SWFAC
      REAL TCARBO
      REAL TGROLF
      REAL TGROPAN
      REAL TGROSTM
      REAL TILN
      REAL TILSW
      REAL TLAG1
      REAL TLFWT
      REAL TLNO
      REAL TMNC
      REAL TPANWT
      REAL TPLA
      REAL TPLAG
      REAL TSIZE
      REAL TSTMWT
      REAL TURFAC
      REAL VANC
      REAL VMNC
      REAL WSTR1
      REAL XNTI
      REAL RGSET
      REAL CRGSET
      REAL CURV
      REAL SETTC(4)
      CHARACTER*3  SETTYP
      REAL      TEMPST
      REAL      XRGSET
C------------------------------------------------------------------------
C   VARIABLES IN OLD COMMON BLOCKS
C------------------------------------------------------------------------

      REAL APTNUP
      REAL BIOMAS
      REAL BIOMS2
      REAL CANNAA
      REAL CANWAA
      REAL CSD1
      REAL CSD2
      REAL CNSD1
      REAL CNSD2
      INTEGER CTYPE
      REAL CUMDEP
      REAL DLAYR(NL)
      REAL DTT
      REAL EMAT
      REAL GNUP
      REAL GPP
      REAL GPSM
      REAL GRAINN
      REAL GRNWT
      INTEGER ICSDUR
      CHARACTER*1 IDETO
      INTEGER IDUR1
      INTEGER IPRINT
      INTEGER ISDATE
      INTEGER ISM
      INTEGER ISTAGE
      CHARACTER*1 ISWNIT
      CHARACTER*1 ISWWAT
      REAL LAI
      INTEGER LEAFNO
      REAL LL(NL)
      REAL MAXLAI
      INTEGER MDATE
      INTEGER NOUTDO
      INTEGER NLAYR
      REAL P1
      REAL P2O
      REAL P2R
      REAL P3
      REAL P4
      REAL P5
      REAL P9
      REAL PANWT
      REAL PGRNWT
      REAL PHINT
      REAL PLTPOP
      REAL PWA
      REAL ROPT
      REAL RTDEP
      REAL SDEPTH
      REAL SI1(6)
      REAL SI2(6)
      REAL SI3(6)
      REAL SI4(6)
      REAL SIND
      REAL SKERWT
      REAL SNOW
      REAL SRAD
      REAL STMWT
      REAL STOVER
      REAL STOVN
      REAL SUMDTT
      REAL SW(NL)
      REAL TANC
      REAL TBASE
      REAL TEMPCR
      REAL TMAX
      REAL TMIN
      REAL TOPT
      REAL TOTNUP
      REAL TPSM
      REAL YIELD
      REAL XGNP
      REAL XSTAGE

C------------------------------------------------------------------------
C     LOCAL VARIABLES
C------------------------------------------------------------------------
      INTEGER   STGDOY(20),YRDOY,I,NDAS,L,L0
      INTEGER   NCOUNT
      REAL      XANC,TEMPCN,TEMPCX,XS,
     +          SWSD,PDTT,RATEIN,PFLOWR,
     +          YIELDB,TWILEN
      REAL      MAXLAI2,MAXLAI3
      REAL      DAYL,ACOEF,TNSOIL,TDSOIL,TMSOIL
      REAL      TH,DOPT

!     Variables needed to send message to WARNING.OUT file
      CHARACTER*78 MESSAGE(10)

C--------------------------------------------------------------------
C                         MAIN CODE
C--------------------------------------------------------------------
      XANC   = TANC*100.0
      APTNUP = STOVN*10.0*PLTPOP
      TOTNUP = APTNUP

      !--------------------------------------------------------------
      !Compute thermal time based on new method developed by J.T.R
      !at CYMMIT, 5/5/98.  TBASE, TOPT, and ROPT are read in from
      ! the species file, as is CTYPE (Cereal type).
      !--------------------------------------------------------------

      !Initially, set TEMPCN and TEMPCS to TMAX and TMIN

      TEMPCN = TMIN
      TEMPCX = TMAX
      XS     = SNOW
      XS     = AMIN1 (XS,15.0)

      !--------------------------------------------------------------
      ! Calculate crown temperature based on temperature and snow cover
      !--------------------------------------------------------------
      IF (TMIN .LT. 0.0) THEN
         TEMPCN = 2.0 + TMIN*(0.4+0.0018*(XS-15.0)**2)
      ENDIF
      IF (TMAX .LT. 0.0) THEN
         TEMPCX = 2.0 + TMAX*(0.4+0.0018*(XS-15.0)**2)
      ENDIF
      TEMPCR = (TEMPCX + TEMPCN)/2.0
      !
      ! DOPT, Devlopment optimum temperature, is set to TOPT during
      ! vegetative growth and to ROPT after anthesis
      !
      DOPT = TOPT
      IF ((ISTAGE .GT. 3) .AND. (ISTAGE .LE. 6)) THEN
         DOPT = ROPT
      ENDIF
      !
      ! Check basic temperature ranges and calculate DTT for
      ! develpment based on PC with JTR
      !
      IF (TMAX .LT. TBASE) THEN
         DTT = 0.0
       ELSEIF (TMIN .GT. DOPT) THEN
         !
         ! This statement replaces DTT = TOPT
         ! ... GoL and LAH, CIMMYT, 1999
         DTT = DOPT - TBASE
         !
         ! Now, modify TEMPCN, TEMPCX based on soil conditions or snow
         ! If wheat and barley is before terminal spiklett stage
         ! Or if corn and sorghum are before 10 leaves
         !
       ELSEIF ((CTYPE .EQ. 1 .AND. LEAFNO .LE. 10) .OR.
     &         (CTYPE .EQ. 2 .AND. ISTAGE .LT. 2)) THEN
         !
         ! Check for snow  (should following be GT.0 or GT.15 ?).
         ! Based on snow cover, calculate DTT for the day
         !
         IF (XS .GT. 0.0) THEN
            !
            ! Snow on the ground
            !
            DTT    = (TEMPCN + TEMPCX)/2.0 - TBASE
          ELSE
            !
            ! No snow, compute soil temperature
            !
            ACOEF  = 0.01061 * SRAD + 0.5902
            TDSOIL = ACOEF * TMAX + (1.0 - ACOEF) * TMIN
            TNSOIL = 0.36354 * TMAX + 0.63646 * TMIN
            IF (TDSOIL .LT. TBASE) THEN
               DTT = 0.0
             ELSE
               IF (TNSOIL .LT. TBASE) THEN
                  TNSOIL = TBASE
               ENDIF
               IF (TDSOIL .GT. DOPT) THEN
                  TDSOIL = DOPT
               ENDIF
!     chp - import DAYL from Weather module instead
               !CALL DAYLEN (DOY,XLAT,DAYL,DEC,SNDN,SNUP)
               TMSOIL = TDSOIL * (DAYL/24.) + TNSOIL * ((24.-DAYL)/24.)
               IF (TMSOIL .LT. TBASE) THEN
                  DTT = (TBASE+TDSOIL)/2.0 - TBASE
                ELSE
                  DTT = (TNSOIL+TDSOIL)/2.0 - TBASE
               ENDIF
               !
               ! Statement added ... GoL and LAH, CIMMYT, 1999
               !
               DTT = AMIN1 (DTT,DOPT-TBASE)
            ENDIF
         ENDIF
         !
         ! Now, compute DTT for when Tmax or Tmin out of range
         !
       ELSEIF (TMIN .LT. TBASE .OR. TMAX .GT. DOPT) THEN
          DTT = 0.0
          DO I = 1, 24
             TH = (TMAX+TMIN)/2. + (TMAX-TMIN)/2. * SIN(3.14/12.*I)
             IF (TH .LT. TBASE) THEN
                TH = TBASE
             ENDIF
             IF (TH .GT. DOPT) THEN
                TH = DOPT
             ENDIF
             DTT = DTT + (TH-TBASE)/24.0
          END DO
       ELSE
          DTT = (TMAX+TMIN)/2.0 - TBASE
      ENDIF

      SUMDTT  = SUMDTT  + DTT
      CUMDTT = CUMDTT + DTT
      CSD1   = CSD1 + 1.0 - SWFAC
      CSD2   = CSD2 + 1.0 - TURFAC
      CNSD1 = CNSD1 + 1.0-NSTRES
      CNSD2 = CNSD2 + 1.0-AGEFAC
      ICSDUR = ICSDUR + 1

C     KJB added
      TEMPST = (TMAX + TMIN)*0.5   !Mean air temperature, C

C--------------------------------------------------------------------
C       Definition of Stages
C
C     7 - Sowing date
C     8 - Germination
C     9 - Emergence
C     1 - End juvenile
C     2 - Pannicle initiation
C     3 - End leaf growth
C     4 - End pannicle growth
C     5 - Grain fill
C     6 - Maturity
C--------------------------------------------------------------------

C--------------------------------------------------------------------
C         ISTAGE = 7: Determine when sowing date occurs
C--------------------------------------------------------------------
      IF (ISTAGE .EQ. 7) THEN
          STGDOY(ISTAGE) = YRDOY
          NDAS           = 0.0
          CRGSET          = 0.0
          XRGSET = 1.0

          CALL ML_PHASEI (
     &    AGEFAC, BIOMAS, BIOMS1, BIOMS2, CNSD1, CNSD2,
     &    CSD2, CSD1, CUMDEP, CUMDTT, CUMPH, DLAYR, NCOUNT,
     &    DTT, EMAT, P5, G4, G5, XRGSET, GPP, GRAINN, GRNWT,
     &    GROLF, GRORT, GROSTM, ICSDUR, IDUR1, ISM,
     &    ISTAGE, ISWWAT,ISWNIT,LAI, LEAFNO, LFWT, MGROLF, MGROPAN,
     &    MGROSTM,MLAG1, MLFWT, MPANWT, MPLAG, MPLA, MSTMWT, NLAYR,
     &    NSTRES, P3, P4, P9, PAF, PANWT, PGC, PHINT, PLA, PLAN, PLAO,
     &    PLATO, PLAY, PLAMX, PLTPOP, PTF, PWA, RANC, RESERVE,
     &    RLV, ROOTN, RTDEP, RTWT, RWU, SDEPTH, SEEDRV, SENLA,
     &    SIND, SLAN, STMWT, STOVWT,
     &    TANC, TCARBO, TGROLF, TGROPAN,
     &    TGROSTM, TILN, TILSW, TLAG1, TLFWT, TLNO, TMNC, TPANWT,
     &    TPLA, TPLAG, TPSM, TSIZE, TSTMWT, SUMDTT, VANC,
     &    VMNC, WSTR1, XNTI)

          IF (ISWWAT .EQ. 'N') RETURN

          CUMDEP         = 0.0
          DO L = 1, NLAYR
             CUMDEP = CUMDEP + DLAYR(L)
             IF (SDEPTH .LT. CUMDEP) GO TO 100       ! Was EXIT
          END DO

  100     CONTINUE                                   ! Sun Fix

          L0 = L
          RETURN

C--------------------------------------------------------------------
C      ISTAGE = 8: Determine Germination Date
C--------------------------------------------------------------------
       ELSEIF (ISTAGE .EQ. 8) THEN
          IF (ISWWAT .NE. 'N') THEN
             IF (SW(L0) .LE. LL(L0)) THEN
                 SWSD = (SW(L0)-LL(L0))*0.65 + (SW(L0+1)-LL(L0+1))*0.35
                 NDAS = NDAS + 1

                 IF (SWSD .LT. SWCG) RETURN
             ENDIF
          ENDIF

          STGDOY(ISTAGE) = YRDOY
          IPRINT = 0

          CALL ML_PHASEI (
     &    AGEFAC, BIOMAS, BIOMS1, BIOMS2, CNSD1, CNSD2,
     &    CSD2, CSD1, CUMDEP, CUMDTT, CUMPH, DLAYR, NCOUNT,
     &    DTT, EMAT, P5, G4, G5, XRGSET, GPP, GRAINN, GRNWT,
     &    GROLF, GRORT, GROSTM, ICSDUR, IDUR1, ISM,
     &    ISTAGE, ISWWAT,ISWNIT,LAI, LEAFNO, LFWT, MGROLF, MGROPAN,
     &    MGROSTM,MLAG1, MLFWT, MPANWT, MPLAG, MPLA, MSTMWT, NLAYR,
     &    NSTRES, P3, P4, P9, PAF, PANWT, PGC, PHINT, PLA, PLAN, PLAO,
     &    PLATO, PLAY, PLAMX, PLTPOP, PTF, PWA, RANC, RESERVE,
     &    RLV, ROOTN, RTDEP, RTWT, RWU, SDEPTH, SEEDRV, SENLA,
     &    SIND, SLAN, STMWT, STOVWT,
     &    TANC, TCARBO, TGROLF, TGROPAN,
     &    TGROSTM, TILN, TILSW, TLAG1, TLFWT, TLNO, TMNC, TPANWT,
     &    TPLA, TPLAG, TPSM, TSIZE, TSTMWT, SUMDTT, VANC,
     &    VMNC, WSTR1, XNTI)

          RETURN

C--------------------------------------------------------------------
C      ISTAGE = 9: Determine seedling emergence date
C--------------------------------------------------------------------
       ELSEIF (ISTAGE .EQ. 9) THEN
          NDAS = NDAS + 1
          RTDEP  = RTDEP + 0.15*DTT
          IDUR1  = IDUR1 + 1

          IF (SUMDTT .LT. P9) RETURN
          IF (P9 .GT. DGET) THEN
             ISTAGE = 6
             PLTPOP = 0.00
             GPP    = 1.0
             GRNWT  = 0.0

             WRITE(MESSAGE(1),1399)
             CALL WARNING(1, 'MLPHEN', MESSAGE)

             WRITE (     *,1399)
             IF (IDETO .EQ. 'Y') THEN
                WRITE (NOUTDO,1399)
             ENDIF
             MDATE = YRDOY
             RETURN
          ENDIF

          STGDOY(ISTAGE) = YRDOY

          CALL ML_PHASEI (
     &    AGEFAC, BIOMAS, BIOMS1, BIOMS2, CNSD1, CNSD2,
     &    CSD2, CSD1, CUMDEP, CUMDTT, CUMPH, DLAYR, NCOUNT,
     &    DTT, EMAT, P5, G4, G5, XRGSET, GPP, GRAINN, GRNWT,
     &    GROLF, GRORT, GROSTM, ICSDUR, IDUR1, ISM,
     &    ISTAGE, ISWWAT,ISWNIT,LAI, LEAFNO, LFWT, MGROLF, MGROPAN,
     &    MGROSTM,MLAG1, MLFWT, MPANWT, MPLAG, MPLA, MSTMWT, NLAYR,
     &    NSTRES, P3, P4, P9, PAF, PANWT, PGC, PHINT, PLA, PLAN, PLAO,
     &    PLATO, PLAY, PLAMX, PLTPOP, PTF, PWA, RANC, RESERVE,
     &    RLV, ROOTN, RTDEP, RTWT, RWU, SDEPTH, SEEDRV, SENLA,
     &    SIND, SLAN, STMWT, STOVWT,
     &    TANC, TCARBO, TGROLF, TGROPAN,
     &    TGROSTM, TILN, TILSW, TLAG1, TLFWT, TLNO, TMNC, TPANWT,
     &    TPLA, TPLAG, TPSM, TSIZE, TSTMWT, SUMDTT, VANC,
     &    VMNC, WSTR1, XNTI)
          RETURN

C--------------------------------------------------------------------
C        ISTAGE = 1: Determine end of juvenile stage
C--------------------------------------------------------------------
       ELSEIF (ISTAGE .EQ. 1) THEN
          NDAS   = NDAS + 1
          XSTAGE = 2. * SUMDTT / P1
          IDUR1  = IDUR1 + 1
          IF (SUMDTT .LT. P1) RETURN
          STGDOY(ISTAGE) = YRDOY

C -------------------------------------------------------------------
C         ISTAGE = 2: Determine date of panicle initiation
C--------------------------------------------------------------------
       ELSEIF (ISTAGE .EQ. 2) THEN
          NDAS   = NDAS + 1
          XSTAGE = 2.0 + SIND
          IDUR1  = IDUR1 + 1
          PDTT   = DTT

          IF (ISWWAT .EQ. 'N') THEN
             ICSDUR = ICSDUR + 1
             ENDIF

          IF (ICSDUR .EQ. 1) THEN
             PDTT = SUMDTT - P1
          ENDIF

          !TWILEN = AMAX1 (TWILEN,P2O)
          IF (TWILEN .GT. P2O) THEN
            !RATEIN = 1.0/(68.0+P2R*(TWILEN-P2O))
            RATEIN = 1.0/(DJTI+P2R*(TWILEN-P2O))
          ELSE
            !RATEIN = 1.0/68.0
            RATEIN = 1.0/DJTI
          ENDIF

          SIND = SIND + RATEIN*PDTT
          IF (SIND .LT. 1.0) RETURN
          STGDOY(ISTAGE) = YRDOY

C--------------------------------------------------------------------
C         ISTAGE = 3: Determine end of leaf growth
C--------------------------------------------------------------------
       ELSEIF (ISTAGE .EQ. 3) THEN
          NDAS   = NDAS + 1
          IDUR1 = IDUR1 + 1
          XSTAGE = 3.0 + 2.0*SUMDTT/P3
          IF (SUMDTT .LT. P3) RETURN
          STGDOY(ISTAGE) = YRDOY
          MAXLAI = LAI

          !chp added these - Bill - check!!!! 07/02/03
              CANNAA = STOVN*PLTPOP
              CANWAA = BIOMAS

! KJB - to add temp effect on seed set
              NCOUNT = 0
              CRGSET = 0.0
C--------------------------------------------------------------------
C         ISTAGE = 4: Determine end of panicle growth
C--------------------------------------------------------------------
       ELSEIF (ISTAGE .EQ. 4) THEN
          NDAS = NDAS + 1
          PFLOWR = 50.0
          XSTAGE = 5.0 + SUMDTT/PFLOWR
! KJB - to add temp effect on seed set
          NCOUNT = NCOUNT + 1

C  KJB:  ADD TEMP EFFECT ON GRAIN-SET, AVERAGE TEMP PRIOR days of ISTAGE4

      RGSET = CURV(SETTYP,SETTC(1),SETTC(2),SETTC(3),SETTC(4),TEMPST)
      RGSET  = AMAX1(RGSET,0.0)
      RGSET = AMIN1(RGSET,1.0)
          CRGSET = CRGSET + RGSET
          XRGSET = CRGSET/NCOUNT
C  KJB
C  Previous code tried, but failed to compute WSTR1 in PHASEI, always csd1=0 and icsdur=1
C  WSTR1 was not passed in initial code.  No longer use it, as PCARB during ISTAGE4
C  carries water stress effect now via CARBO
          IF (ISWWAT .EQ. 'N') THEN
              WSTR1 = 0.0
          ELSE
             WSTR1 = CSD1/ICSDUR
          ENDIF
C  KJB
          IF (SUMDTT .LE. PFLOWR) THEN
             IDUR1 = IDUR1 + 1
             ENDIF

          IF (SUMDTT .GE. PFLOWR .AND. IPRINT .EQ. 0) THEN
             STGDOY(ISTAGE) = YRDOY
             BIOMS2 = BIOMAS/PLTPOP
             ISDATE = YRDOY
             MAXLAI2 = LAI
             IPRINT = 1
             ENDIF

          IF (SUMDTT .LT. P4) RETURN
             STGDOY(16) = YRDOY
             XSTAGE     = 6.0 + 0.5 * (SUMDTT - 50.0) / (P4 - 50.0)
          IF (GPP .LE. 0.0) THEN
             GPP   = 1.0
             ENDIF

          MAXLAI3 = LAI
!chp          MAXLAI  = AMAX1(MAXLAI1,MAXLAI2,MAXLAI3)
          MAXLAI  = AMAX1(MAXLAI2,MAXLAI3)

C--------------------------------------------------------------------
C        ISTAGE = 5: Determine end of effective filling period
C--------------------------------------------------------------------
       ELSEIF (ISTAGE .EQ. 5) THEN
          NDAS = NDAS + 1

C
C            Determine end of grain filling for millet and sorghum
C
          XSTAGE = 6.5 + 2.5*SUMDTT/P5
          IF (SUMDTT .LT. P5) RETURN

          IF (ISM .GT. 0) THEN
             IF (SUMDTT .LT. P5 + 90.0) RETURN
                 EMAT = EMAT + 1
                 IF (EMAT .LT. 2) RETURN
                 STGDOY (ISTAGE) = YRDOY
             ELSE
                 ISM         = 1
                 STGDOY (17) = YRDOY
             RETURN
             ENDIF

C--------------------------------------------------------------------
C       ISTAGE = 6: Determine physiological maturity
C--------------------------------------------------------------------
       ELSEIF (ISTAGE .EQ. 6) THEN
          IF (DTT    .LT. 2.0) RETURN
          IF (SUMDTT .LT. 2.0) RETURN

          STGDOY(ISTAGE) = YRDOY
          IPRINT         = 0
          MDATE          = YRDOY
          GRNWT  = PANWT - 1.1 * PWA
          GRNWT  = AMAX1(GRNWT,0.0)
          YIELD  = GRNWT*10.0*PLTPOP

          IF (PLTPOP .NE. 0.0) THEN
              IF (GPP .LE. 0.0) THEN
                  GPP = 1.0
                  ENDIF
              SKERWT = GRNWT/GPP
              GPSM = GPP*PLTPOP
              STOVER  = BIOMAS*10.0-YIELD

              YIELDB  = YIELD/62.8
              PGRNWT = SKERWT*1000.0
              IF (ISWNIT .EQ. 'N') THEN
                  XGNP = 0.0
                  GNUP = 0.0
                ELSE
                  IF (GRNWT .GT. 0.0) THEN
                   XGNP = (GRAINN/GRNWT)*100.0
                   GNUP = GRAINN*PLTPOP*10.0
                   ENDIF
              ENDIF

             TOTNUP      = GNUP + APTNUP

             SI1(ISTAGE) = 0.0
             SI2(ISTAGE) = 0.0
             SI3(ISTAGE) = 0.0
             SI4(ISTAGE) = 0.0
          ENDIF
      ENDIF       !End of main ISTAGE Logic

C--------------------------------------------------------------------
C     This code is run at the beginning of ISTAGE 1, 2, 3, 4, 5, 6
C
C     This code is NOT run for ISTAGE 7, 8, or 9 due to the RETURN
C     statements in the above sections of code for ISTAGE 7, 8, and 9
C--------------------------------------------------------------------

      IF (ISWWAT .NE. 'N'.AND.ISTAGE.GT.0) THEN
         SI1(ISTAGE) = CSD1  / ICSDUR
         SI2(ISTAGE) = CSD2  / ICSDUR
         SI3(ISTAGE) = CNSD1 / ICSDUR
         SI4(ISTAGE) = CNSD2 / ICSDUR
      ENDIF

      CALL ML_PHASEI (
     &    AGEFAC, BIOMAS, BIOMS1, BIOMS2, CNSD1, CNSD2,
     &    CSD2, CSD1, CUMDEP, CUMDTT, CUMPH, DLAYR, NCOUNT,
     &    DTT, EMAT, P5, G4, G5, XRGSET, GPP, GRAINN, GRNWT,
     &    GROLF, GRORT, GROSTM, ICSDUR, IDUR1, ISM,
     &    ISTAGE, ISWWAT,ISWNIT,LAI, LEAFNO, LFWT, MGROLF, MGROPAN,
     &    MGROSTM,MLAG1, MLFWT, MPANWT, MPLAG, MPLA, MSTMWT, NLAYR,
     &    NSTRES, P3, P4, P9, PAF, PANWT, PGC, PHINT, PLA, PLAN, PLAO,
     &    PLATO, PLAY, PLAMX, PLTPOP, PTF, PWA, RANC, RESERVE,
     &    RLV, ROOTN, RTDEP, RTWT, RWU, SDEPTH, SEEDRV, SENLA,
     &    SIND, SLAN, STMWT, STOVWT,
     &    TANC, TCARBO, TGROLF, TGROPAN,
     &    TGROSTM, TILN, TILSW, TLAG1, TLFWT, TLNO, TMNC, TPANWT,
     &    TPLA, TPLAG, TPSM, TSIZE, TSTMWT, SUMDTT, VANC,
     &    VMNC, WSTR1, XNTI)

      RETURN

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

1399  FORMAT (10X,'Seed ran out of metabolite due to deep planting')

      END SUBROUTINE ML_PHENOL

ext: for


text:
      SUBROUTINE STABL3(ITERM)
*
*
*       Stability test of three-body system.
*       ------------------------------------
*
      IMPLICIT REAL*8  (A-H,M,O-Z)
      COMMON/AZREG/  TIME3,TMAX,Q(8),P(8),R1,R2,R3,ENERGY,M(3),X3(3,3),
     &               XDOT3(3,3),CM(10),C11,C12,C19,C20,C24,C25,
     &               NSTEP3,NAME3(3),KZ15,KZ27
      COMMON/AZOUT/  MB,RB,RI,SEMI,SEMI1,E,E1
*
*
*       Transform to physical variables.
      CALL TRANS3(3)
*
*       Identify index of second binary component and the third body.
      IM = 1
      IF (R2.LT.R1) IM = 2
      I = 3 - IM
*
*       Form scalar product of distance & velocity for body #I.
      RIDOT = X3(1,I)*XDOT3(1,I) + X3(2,I)*XDOT3(2,I) +
     &                             X3(3,I)*XDOT3(3,I)
*
*       Set distance & radial velocity of body #I with respect to binary.
      MB = M(IM) + M(3)
      FAC = CM(7)/MB
      RI = SQRT(X3(1,I)**2 + X3(2,I)**2 + X3(3,I)**2)
      RIDOT = FAC*RIDOT/RI
      RI = FAC*RI
*
*       Evaluate orbital elements.
      RDOT = 0.0D0
      VREL2 = 0.0D0
      VI2 = 0.0D0
      DO 10 K = 1,3
          RDOT = RDOT +
     &               (X3(K,3) - X3(K,IM))*(XDOT3(K,3) - XDOT3(K,IM))
          VREL2 = VREL2 + (XDOT3(K,3) - XDOT3(K,IM))**2
          VI2 = VI2 + XDOT3(K,I)**2
   10 CONTINUE
*
*       Determine semi-major axis & eccentricity of inner binary.
      RB = MIN(R1,R2)
      SEMI = 2.0D0/RB - VREL2/MB
      SEMI = 1.0/SEMI
      E = SQRT((1.0D0 - RB/SEMI)**2 + RDOT**2/(SEMI*MB))
*
*       Form outer semi-major axis & eccentricity.
      VI2 = VI2*FAC**2
      SEMI1 = 2.0D0/RI - VI2/CM(7)
      SEMI1 = 1.0/SEMI1
      E1 = SQRT((1.0D0 - RI/SEMI1)**2 + (RI*RIDOT)**2/(SEMI1*CM(7)))
*
*       Obtain standard stability ratio (outer pericentre / inner apocentre).
      RATIO = SEMI1*(1.0D0 - E1)/(SEMI*(1.0D0 + E))
*
*       Form coefficients for stability test (Valtonen, Vistas Ast 32, 1988).
*     AM = (2.65 + E)*(1.0 + M(I)/MB)**0.3333
*     FM = (2.0*M(I) - MB)/(3.0*MB)
*
*       Expand natural logarithm for small arguments.
*     IF (ABS(FM).LT.0.67) THEN
*         BM = FM*(1.0 - (0.5 - ONE3*FM)*FM)
*     ELSE
*         BM = LOG(1.0D0 + FM)
*     END IF
*       Define mass dependent criterion of Harrington (A.J. 80) & Bailyn.
*     PCRIT = AM*(1.0 + 0.7*BM)*SEMI
*
*       Form hierarchical stability ratio (Kiseleva & Eggleton 1995).
*     Q0 = MB/M(I)
*     Q1 = MAX(M(3)/M(IM),M(IM)/M(3))
*     Q3 = Q0**0.33333
*     Q13 = Q1**0.33333
*     AR = 1.0 + 3.7/Q3 - 2.2/(1.0 + Q3) + 1.4/Q13*(Q3 - 1.0)/(Q3 + 1.0)
*     PCRIT = AR*SEMI*(1.0D0 + E)
*
*       Adopt the semi-analytical stability criterion (MA 1997).
      Q1 = M(I)/MB
      IF (E1.LT.1.0) THEN
          XFAC = (1.0 + Q1)*(1.0 + E1)/SQRT(1.0 - E1)
      ELSE
          XFAC = 40.0*(1.0 + Q1)
      END IF
      PCRIT = 2.8*XFAC**0.4*SEMI
      PMIN = SEMI1*(1.0D0 - E1)
*
*       Set negative termination index if system is stable and RB > SEMI.
      IF (PCRIT.LT.PMIN.AND.E1.LT.1.0.AND.1./RB.LT.1./SEMI) THEN
          ITERM = -1
*       Obtain Zare's stability parameter (valid for small inclinations).
          M1 = M(IM)
          M2 = M(3)
          M3 = M(I)
          CALL STABLZ(M1,M2,M3,SP)
          if(rank.eq.0)
     &    WRITE (6,20)  SEMI, SEMI1, E, E1, RI, RATIO, SP, PCRIT, PMIN
   20     FORMAT ('  STABT:    A A1 E E1 RI RATIO SP PCR PM ',
     &                           1P,2E10.2,0P,2F6.2,F9.5,2F6.2,1P,2E9.1)
*       Terminate if escaper is outside 3*SEMI (> 0).
      ELSE IF (E1.GT.1.0.AND.1./RI.LT.1./(3.0*SEMI)) THEN
          ITERM = -1
      ELSE
          ITERM = 0
      END IF
*
      RETURN
*
      END

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine dlwqn7 ( a     , j     , c     , lun   , lchar  ,
     &                    action, dlwqd , gridps)

!       Deltares Software Centre

!>\file
!>                         horizontally upwind, vertically central, direct stationary method (7)
!>
!>                         Stationary solution. Upwind 1st order horizontally, central
!>                         vertically. Fully implicit with a direct method.\n
!>                         Matrices become very large in 3D and method unworkable. In 2D
!>                         the method can be used. In 1D the method outperforms the
!>                         iterative methods.

!     CREATED            : june 1988 by L. Postma
!
!     LOGICAL UNITS      : LUN(19) , output, monitoring file
!                          LUN(20) , output, formatted dump file
!                          LUN(21) , output, unformatted hist. file
!                          LUN(22) , output, unformatted dump file
!
!     SUBROUTINES CALLED : DLWQTR, user transport routine
!                          DLWQWQ, user waterquality routine
!                          DLWQPP, user postprocessing routine
!                          DLWQ10, system monitoring routine
!                          DLWQ11, system dump routine
!                          DLWQ13, system postpro-dump routine
!                          DLWQ15, wasteload routine
!                          DLWQ60, scales water quality
!                          DLWQ61, clears the matrix
!                          DLWQ63, stores the results
!                          DLWQ65, computes closure error
!                          DLWQ66, makes masses
!                          DLWQ67, zeros the matrix
!                          DLWQ70, fills the matrix
!                          DLWQ71, performs mass balance computation
!                          DELMAT, inverts the matrix
!                          DHOPNF, opens files
!
!     PARAMETERS    :
!
!     NAME    KIND     LENGTH   FUNC.  DESCRIPTION
!     ---------------------------------------------------------
!     A       REAL       *      LOCAL  real      workspace array
!     J       INTEGER    *      LOCAL  integer   workspace array
!     C       CHARACTER  *      LOCAL  character workspace array
!     LUN     INTEGER    *      INPUT  array with unit numbers
!     LCHAR   CHARACTER  *      INPUT  filenames
!
!     Declaration of arguments
!
      use grids
      use timers
      use delwaq2_data
      use m_openda_exchange_items, only : get_openda_buffer
      use report_progress
      use waqmem          ! module with the more recently added arrays

      implicit none

      include 'actions.inc'
!
!     Declaration of arguments
!
      REAL, DIMENSION(*)          :: A
      INTEGER, DIMENSION(*)       :: J
      INTEGER, DIMENSION(*)       :: LUN
      CHARACTER*(*), DIMENSION(*) :: C
      CHARACTER*(*), DIMENSION(*) :: LCHAR
      INTEGER                     :: ACTION
      TYPE(DELWAQ_DATA)           :: DLWQD
      type(GridPointerColl)       :: GridPs               ! collection off all grid definitions

!
!     COMMON  /  SYSN   /   System characteristics
!
      INCLUDE 'sysn.inc'
!
!     COMMON  /  SYSI  /    Timer characteristics
!
      INCLUDE 'sysi.inc'
!
!     COMMON  /  SYSA   /   Pointers in real array workspace
!
      INCLUDE 'sysa.inc'
!
!     COMMON  /  SYSJ   /   Pointers in integer array workspace
!
      INCLUDE 'sysj.inc'
!
!     COMMON  /  SYSC   /   Pointers in character array workspace
!
      INCLUDE 'sysc.inc'
!
!     Local declarations
!
      LOGICAL         IMFLAG , IDFLAG , IHFLAG
      LOGICAL         LDUMMY , LSTREC , LREWIN

      INTEGER         ITIME
      INTEGER         ITIMEL
      INTEGER         IAFLAG
      INTEGER         IBFLAG
      INTEGER         ISYS
      INTEGER         ICSYS
      INTEGER         NSYS
      INTEGER         INWTYP
      INTEGER         I
      INTEGER          :: NOSSS
      INTEGER          :: NOQTT
      INTEGER         sindex

      integer       :: ithandl

      !
      ! Distinguishing the actions is superfluous:
      ! there is only one step
      !
      IF ( ACTION == ACTION_INITIALISATION .OR.
     &     ACTION == ACTION_FINALISATION        ) THEN
          RETURN
      ENDIF
!
!          some initialisation
!
      ithandl = 0
      if ( timon ) call timstrt ( "dlwqn7", ithandl )

      ITIMEL = ITSTRT
      ITIME  = ITSTRT+IDT
      IBFLAG = 0
      IF ( MOD(INTOPT,16) .GE. 8 ) IBFLAG = 1
      CALL ZERO ( A(IMAS2) , NOTOT*5 )
      LDUMMY = .FALSE.
      LSTREC = .FALSE.
      nosss  = noseg + nseg2
      NOQTT  = NOQ + NOQ4
      inwtyp = intyp + nobnd

!        Determine the volumes and areas that ran dry,
!        They cannot have explicit processes during this time step

         call hsurf  ( noseg    , nopa     , c(ipnam) , a(iparm) , nosfun   ,
     &                 c(isfna) , a(isfun) , surface  , sindex   , lun(19)  )
         call dryfld ( noseg    , nosss    , nolay    , a(ivol)  , noq1+noq2,
     &                 a(iarea) , nocons   , c(icnam) , a(icons) , sindex   ,
     &                 surface  , j(iknmr) , iknmkv   )
!
!          makes closure error
!
      IF ( IDT.EQ.0 ) THEN
         CALL ZERO ( A(IVOL2), NOSEG )
      ELSE IF ( J(INRH2+1).GE.0 .AND. IVFLAG.EQ.0 ) THEN
         CALL DLWQ41 ( LUN     , ITIME   , ITIMEL  , A(IHARM), A(IFARR),
     *                 J(INRHA), J(INRH2), J(INRFT), NOSEG   , A(IVOL2),
     *                 J(IBULK), LCHAR   , ftype   , ISFLAG  , IVFLAG  ,
     *                 LDUMMY  , J(INISP), A(INRSP), J(INTYP), J(IWORK),
     *                 LSTREC  , LREWIN  , A(IVOLL), MYPART  , dlwqd   )
         CALL DLWQ65 ( A(IVOL2), A(IVOL) , IDT     , NOSEG   )
      ELSE
         CALL ZERO ( A(IVOL2) , NOSEG )
         WRITE ( LUN(19), 1000 )
      ENDIF
!
!          loop over the systems
!
      NSYS   = 1
      IAFLAG = 1
      DO 10 ISYS = 1 , NOSYS
         IF ( ISYS .EQ. NOSYS ) NSYS   = 1 + NOTOT - NOSYS
!
!             do the user transport processes
!
         ICSYS = ISYS
         CALL DLWQTR ( NOTOT   , NOSYS   , NOSEG   , NOQ     , NOQ1    ,
     *                 NOQ2    , NOQ3    , NOPA    , NOSFUN  , NODISP  ,
     *                 NOVELO  , J(IXPNT), A(IVOL) , A(IAREA), A(IFLOW),
     *                 A(ILENG), A(ICONC), A(IDISP), A(ICONS), A(IPARM),
     *                 A(IFUNC), A(ISFUN), A(IDIFF), A(IVELO), ICSYS   ,
     *                 IDT     , C(ISNAM), NOCONS  , NOFUN   , C(ICNAM),
     *                 C(IPNAM), C(IFNAM), C(ISFNA), LDUMMY  , ILFLAG  ,
     *                 NPARTp  )
!
!             do the user water quality processes
!
         CALL DLWQWQ ( NOTOT   , NOSYS   , NOSEG   , NOPA    , NOSFUN  ,
     *                 A(IVOL) , A(ICONC), A(ICONS), A(IPARM), A(IFUNC),
     *                 A(ISFUN), A(IDERV), ICSYS   , IDT     , A(ISMAS),
     *                 IBFLAG  , C(ISNAM), NOCONS  , NOFUN   , C(ICNAM),
     *                 C(IPNAM), C(IFNAM), C(ISFNA), NODUMP  , J(IDUMP))
         CALL DLWQ60 ( A(IDERV), A(ICONC), NOTOT   , NOSEG   , ITFACT  ,
     *                 A(IMAS2), ISYS    , NSYS    , A(IDMPS), INTOPT  ,
     *                 J(ISDMP))
!
!             add the waste loads
!
         call dlwq15 ( nosys     , notot    , noseg    , noq      , nowst    ,
     &                 nowtyp    , ndmps    , intopt   ,     1    , itime    ,
     &                 iaflag    , c(isnam) , a(iconc) , a(ivol)  , a(ivol2) ,
     &                 a(iflow ) , j(ixpnt) , c(iwsid) , c(iwnam) , c(iwtyp) ,
     &                 j(inwtyp) , j(iwast) , iwstkind , a(iwste) , a(iderv) ,
     &                 iknmkv    , nopa     , c(ipnam) , a(iparm) , nosfun   ,
     &                 c(isfna ) , a(isfun) , j(isdmp) , a(idmps) , a(imas2) ,
     &                 a(iwdmp)  , isys     , nsys     , j(iowns ), mypart   )
!
!          fill the matrix
!
         CALL DLWQ61 ( A(ICONC), A(IDERV), A(IVOL2), A(ITIMR), NOSEG   ,
     *                           NOTOT   , ISYS    , NSYS    , JTRACK  )
         CALL DLWQ70 ( A(IDISP), A(IDIFF), A(IAREA), A(IFLOW), A(ILENG),
     *                 A(IVELO), A(IBOUN), J(IXPNT), NOTOT   , ISYS    ,
     *                 NSYS    , NOQ1    , NOQ2    , NOQ     , NODISP  ,
     *                 NOVELO  , J(IDPNT), J(IVPNT), A(IDERV), A(ITIMR),
     *                                     JTRACK  , INTOPT  , ILFLAG  )
         CALL DLWQ67 ( A(ITIMR), NOSEG   , JTRACK  )
!
!             invert the matrix and store the results
!
         CALL DELMAT ( NOSEG   , JTRACK  , JTRACK  , NSYS    , A(ITIMR),
     *                                               A(IDERV),    0    )
         CALL DLWQ63 ( A(ICONC), A(IDERV), A(IMAS2), NOSEG   , NOTOT   ,
     *                 ISYS    , NSYS    , A(IDMPS), INTOPT  , J(ISDMP))
   10 CONTINUE
!
!          mass balance
!
      IAFLAG = 1
      CALL DLWQ71 ( A(IDISP), A(IDIFF), A(IAREA), A(IFLOW), A(ILENG),
     *              A(IVELO), A(ICONC), A(IBOUN), J(IXPNT), NOSYS   ,
     *              NOTOT   , NOQ1    , NOQ2    , NOQ     , NODISP  ,
     *              NOVELO  , J(IDPNT), J(IVPNT), INTOPT  , A(IMAS2),
     *              ILFLAG  , A(IDMPQ), NDMPQ   , J(IQDMP))
      CALL DLWQ66 ( A(IDERV), A(IVOL) , A(ICONC), NOTOT   , NOSEG   )
!
!     Call OUTPUT system
!
      CALL DLWQO2 ( NOTOT   , NOSEG   , NOPA    , NOSFUN  , ITSTRT  ,
     +              C(IMNAM), C(ISNAM), C(IDNAM), J(IDUMP), NODUMP  ,
     +              A(ICONC), A(ICONS), A(IPARM), A(IFUNC), A(ISFUN),
     +              A(IVOL) , NOCONS  , NOFUN   , 1       , NOUTP   ,
     +              LCHAR   , LUN     , J(IIOUT), J(IIOPO), A(IRIOB),
     +              C(IOSNM), C(IOUNI), C(IODSC), C(ISSNM), C(ISUNI), C(ISDSC),
     +              C(IONAM), NX      , NY      , J(IGRID), C(IEDIT),
     +              NOSYS   , A(IBOUN), J(ILP)  , A(IDERV), A(IMAS2),
     +              A(ISMAS), NFLUX   , A(IFLXI), ISFLAG  , IAFLAG  ,
     +              IBFLAG  , IMSTRT  , IMSTOP  , IMSTEP  , IDSTRT  ,
     +              IDSTOP  , IDSTEP  , IHSTRT  , IHSTOP  , IHSTEP  ,
     +              IMFLAG  , IDFLAG  , IHFLAG  , NOLOC   , A(IPLOC),
     +              NODEF   , A(IDEFA), ITSTRT  , ITSTOP  , NDMPAR  ,
     +              C(IDANA), NDMPQ   , NDMPS   , J(IQDMP), J(ISDMP),
     +              J(IPDMP), A(IDMPQ), A(IDMPS), A(IFLXD), NTDMPQ  ,
     +              C(ICBUF), NORAAI  , NTRAAQ  , J(IORAA), J(NQRAA),
     +              J(IQRAA), A(ITRRA), C(IRNAM), A(ISTOC), NOGRID  ,
     +              NOVAR   , J(IVARR), J(IVIDX), J(IVTDA), J(IVDAG),
     +              J(IAKND), J(IAPOI), J(IADM1), J(IADM2), J(IVSET),
     +              J(IGNOS), J(IGSEG), A       , NOBND   , NOBTYP  ,
     +              C(IBTYP), J(INTYP), C(ICNAM), NOQ     , J(IXPNT),
     +              INTOPT  , C(IPNAM), C(IFNAM), C(ISFNA), J(IDMPB),
     +              NOWST   , NOWTYP  , C(IWTYP), J(IWAST), J(INWTYP),
     +              A(IWDMP), iknmkv  , J(IOWNS), MYPART  , isegcol )
!
!          close files, except monitor file
!
      call CloseHydroFiles( dlwqd%collcoll )
      call close_files( lun )
!
!          write restart file
!
      CALL DLWQ13 ( LUN      , LCHAR , A(ICONC) , ITSTRT, C(IMNAM) ,
     *              C(ISNAM) , NOTOT , NOSEG    )
!
!          user output routine
!
      CALL DLWQPP ( NOTOT   , NOSYS   , NOSEG   , NOPA    , NOSFUN  ,
     *              ITSTRT  , IMFLAG  , IDFLAG  , IHFLAG  , C(IMNAM),
     *              C(ISNAM), C(IDNAM), C(IWSID), J(IDUMP), NODUMP  ,
     *              J(IWAST), NOWST   , A(ICONC), A(ICONS), A(IPARM),
     *              A(IFUNC), A(ISFUN), A(IVOL ), A(IWSTE), A(IBOUN),
     *              NOBND   , ITSTRT  , ITSTOP  , NX      , NY      ,
     *              J(IGRID), NODISP  , NOVELO  , NOQ     , NOQ1    ,
     *              NOQ2    , NOQ3    , A(IDISP), A(IVELO), A(ISMAS),
     *              IBFLAG  , NOCONS  , NOFUN   , C(ICNAM), C(IPNAM),
     *              C(IFNAM), C(ISFNA), C(IBNID))
!
!          output formats
!
 1000 FORMAT ( 'No closure error corrections !' )
!
      if ( timon ) call timstop ( ithandl )
      RETURN
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine calcguesstincf(nface,dmin,vfa,umfa,cvfa,hcfa,ithermal,
     &        tincfguess,compressible)
!
!     calculates a guess for tincf based on the minimum of:
!       dmin/(local velocity)
!       density*dmin**2/(2*dynamic_viscosity)
!
!       dmin is the smallest edge length of the mesh
!
      implicit none
!
      integer nface,i,ithermal,compressible
!
      real*8 vfa(0:5,*),umax,dmin,umfa(*),tincfguess,cvfa(*),hcfa(*)
!
      tincfguess=1.d30
      do i=1,nface
!
!        convection
!
         umax=dsqrt(vfa(1,i)*vfa(1,i)+
     &              vfa(2,i)*vfa(2,i)+
     &              vfa(3,i)*vfa(3,i))
         if(umax.gt.1.d-30) tincfguess=min(tincfguess,dmin/umax)
!
!        viscous diffusion
!
         if((umfa(i).gt.0.d0).and.(vfa(5,i).gt.0.d0)) then
            tincfguess=min(tincfguess,vfa(5,i)*dmin*dmin/
     &           (2.d0*umfa(i)))
         endif
!
!        thermal diffusion
!
         if(ithermal.gt.0) then
            if((hcfa(i).gt.0.d0).and.(cvfa(i).gt.0.d0).and.
     &         (vfa(5,i).gt.0.d0)) then
               tincfguess=min(tincfguess,vfa(5,i)*cvfa(i)*dmin*dmin/
     &              (2.d0*hcfa(i)))
            endif
         endif
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
!$acc parallel loop collapse (2)
!$acc atomic

max_stars_repo_name: passlab/accparser

max_issues_repo_name: passlab/accparser

max_forks_repo_name: passlab/accparser


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

      SUBROUTINE NODORD (NPNODE, XN, YN, LISTN, NUID, NNN)
C***********************************************************************

C  SUBROUTINE NODORD = ORDER THE NODE TABLE INTO INCREASING VALUES OF
C                      THE VARIABLE LISTN

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

      DIMENSION LISTN (NPNODE)
      DIMENSION XN (NPNODE), YN (NPNODE), NUID (NPNODE)

      NN = NNN
      M = NN
  100 CONTINUE
      M =  (9 * M) / 16
      IF (M .LE. 0) RETURN
      M1 = M + 1
      DO 120 J = M1, NN
         L = J
         I = J - M
  110    CONTINUE
         IF (LISTN (L) .LT. LISTN (I)) THEN
            KLISTN = LISTN (I)
            KNUID = NUID (I)
            TXN = XN (I)
            TYN = YN (I)
            LISTN (I) = LISTN (L)
            NUID (I) = NUID (L)
            XN (I) = XN (L)
            YN (I) = YN (L)
            LISTN (L) = KLISTN
            NUID (L) = KNUID
            XN (L) = TXN
            YN (L) = TYN
            L = I
            I = I - M
            IF (I .GE. 1)GOTO 110
         ENDIF
  120 CONTINUE
      GOTO 100

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
	SUBROUTINE FWINDW(IWINDW)
C
	COMMON /WINBLK/ WINDOW(1024)
	DATA TWOPI /6.2831853/
C
	GO TO (100,200,300) IWINDW+1
	PRINT*,'REQUESTED WINDOW DOES NOT EXIST'
 100	CONTINUE
	PRINT*,'IWINDW =',IWINDW,'  NO WINDOWING'
	DO N = 1,1024
	  WINDOW(N) = 1.
	ENDDO
	RETURN
 200	CONTINUE
C	PRINT*,'IWINDW =',IWINDW,'  HAMMING'
	ALPHA = .54
	GO TO 310
 300	CONTINUE
	ALPHA = .5
	PRINT*,'IWINDW =',IWINDW,'  HANNING'
 310	CONTINUE
	SUMSQ = 0.
	DO N = 1,1024
	  WINDOW(N) = ALPHA + (ALPHA-1.)*COS(TWOPI*(N-1)/1024.)
	  SUMSQ = SUMSQ + WINDOW(N)**2
	ENDDO
	RMS = SQRT(SUMSQ/1024.)
	DO N = 1,1024
	  WINDOW(N) = WINDOW(N)/RMS
	ENDDO
	RETURN
	END


text:
!==========================================================================
elemental subroutine gsw_rho_second_derivatives_wrt_enthalpy (sa, ct, p, &
                                              rho_sa_sa, rho_sa_h, rho_h_h)
! =========================================================================
!
!  Calculates three second-order derivatives of rho with respect to enthalpy.
!  Note that this function uses the using the computationally-efficient
!  expression for specific volume (Roquet et al., 2014).
!
!  SA  =  Absolute Salinity                                        [ g/kg ]
!  CT  =  Conservative Temperature (ITS-90)                       [ deg C ]
!  p   =  sea pressure                                             [ dbar ]
!         ( i.e. absolute pressure - 10.1325 dbar )
!
!  rho_SA_SA = The second-order derivative of rho with respect to
!              Absolute Salinity at constant h & p.     [ J/(kg (g/kg)^2) ]
!  rho_SA_h  = The second-order derivative of rho with respect to
!              SA and h at constant p.                   [ J/(kg K(g/kg)) ]
!  rho_h_h   = The second-order derivative of rho with respect to h at
!              constant SA & p
!--------------------------------------------------------------------------

use gsw_mod_toolbox, only : gsw_specvol
use gsw_mod_toolbox, only : gsw_specvol_first_derivatives_wrt_enthalpy
use gsw_mod_toolbox, only : gsw_specvol_second_derivatives_wrt_enthalpy

use gsw_mod_kinds

implicit none

real (r8), intent(in) :: sa, ct, p
real (r8), intent(out), optional :: rho_sa_sa, rho_sa_h, rho_h_h

integer :: iflag1, iflag2
real (r8) :: rec_v, rec_v2, rec_v3, v_h, v_h_h, v_sa, v_sa_h, v_sa_sa

iflag1 = 0
if (present(rho_sa_sa) .or. present(rho_sa_h)) iflag1 = ibset(iflag1,1)
if (present(rho_sa_h) .or. present(rho_h_h)) iflag1 = ibset(iflag1,2)

call gsw_specvol_first_derivatives_wrt_enthalpy(sa,ct,p,v_sa,v_h,iflag=iflag1)

iflag2 = 0
if (present(rho_sa_sa)) iflag2 = ibset(iflag2,1)
if (present(rho_sa_h)) iflag2 = ibset(iflag2,2)
if (present(rho_h_h)) iflag2 = ibset(iflag2,3)

call gsw_specvol_second_derivatives_wrt_enthalpy(sa,ct,p,v_sa_sa,v_sa_h,v_h_h, &
                                                 iflag=iflag2)

rec_v = 1.0_r8/gsw_specvol(sa,ct,p)
rec_v2 = rec_v**2
rec_v3 = rec_v2*rec_v

if (present(rho_sa_sa)) rho_sa_sa = -v_sa_sa*rec_v2 + 2.0_r8*v_sa*v_sa*rec_v3

if (present(rho_sa_h)) rho_sa_h = -v_sa_h*rec_v2 + 2.0_r8*v_sa*v_h*rec_v3

if (present(rho_h_h)) rho_h_h = -v_h_h*rec_v2 + 2.0_r8*v_h*v_h*rec_v3

return
end subroutine

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

max_stars_repo_name: minsukji/ci-debug

max_issues_repo_name: minsukji/ci-debug

max_forks_repo_name: minsukji/ci-debug


          program demo_unexpandviewport
          use M_draw
          use M_draw, only : page
          use M_draw,    only  : D_BLACK,   D_WHITE
          use M_draw,    only  : D_RED,     D_GREEN,    D_BLUE
          use M_draw,    only  : D_YELLOW,  D_MAGENTA,  D_CYAN
          implicit none
          integer :: ipaws

          !! set up graphics area
          call prefsize(1000,200)
          call vinit(' ') ! start graphics using device $M_draw_DEVICE

          !! draw circle that fills default world coordinate window
          call polyfill(.true.)
          call color(D_GREEN)
          call circle(0.0,0.0,1.0)
          ipaws=getkey() !! pause

          !! set new scales for viewport so <-1,-1> and <1,1> are at
          !! corners of display instead of corners of largest square
          !! that can fit on display
          call expandviewport()
          call viewport(-1.0,1.0,-1.0,1.0)

          !! draw circle that fills default world coordinate window again
          !! instead of a circle, the mapping now produces an ellipse unless
          !! this device has a square display
          call polyhatch(.true.)
          call hatchpitch(0.1)
          call hatchang(30.0)
          call linewidth(40)
          call color(D_CYAN)
          call circle(0.0,0.0,1.0)
          ipaws=getkey() !! pause

          !! set new scales for viewport so <-1,-1> and <1,1> are at
          !! corners of largest square that fits on display
          call unexpandviewport()
          !! actually change to the new viewport
          call viewport(-1.0,1.0,-1.0,1.0)

          !! now the same circle should draw where the original one did
          call color(D_BLACK)
          call hatchang(120.0)
          call linewidth(40)
          call circle(0.0,0.0,1.0)
          ipaws=getkey() !! pause

          call vexit()                  !! wrap up graphics

          end program demo_unexpandviewport