text:
C--- SPARSPAK-A (ANSI FORTRAN) RELEASE III --- NAME = DEGREE
C  (C)  UNIVERSITY OF WATERLOO   JANUARY 1984
C***************************************************************
C***************************************************************
C*******     DEGREE ..... DEGREE IN MASKED COMPONENT     *******
C***************************************************************
C***************************************************************
C
C     PURPOSE - THIS ROUTINE COMPUTES THE DEGREES OF THE NODES
C        IN THE CONNECTED COMPONENT SPECIFIED BY MASK AND ROOT.
C        NODES FOR WHICH MASK IS ZERO ARE IGNORED.
C
C     INPUT PARAMETERS -
C        ROOT   - IS THE INPUT NODE THAT DEFINES THE COMPONENT.
C        (XADJ,ADJNCY) - ADJACENCY STRUCTURE PAIR.
C        MASK   - SPECIFIES A SECTION SUBGRAPH.
C
C     OUTPUT PARAMETERS -
C        DEG    - ARRAY CONTAINING THE DEGREES OF THE NODES IN
C                 THE COMPONENT.
C        CCSIZE - SIZE OF THE COMPONENT SPECIFED BY MASK AND
C                 ROOT.
C
C     WORKING PARAMETER -
C        LS     - A TEMPORARY VECTOR USED TO STORE THE NODES OF
C                 THE COMPONENT LEVEL BY LEVEL.
C
C***************************************************************
C
      SUBROUTINE  DEGREE ( ROOT, XADJ, ADJNCY, MASK, DEG,
     1                     CCSIZE, LS )
C
C***************************************************************
C
         INTEGER    ADJNCY(1), DEG(1)   , LS(1)    , MASK(1)
         INTEGER    XADJ(1)
         INTEGER    CCSIZE, I     , IDEG  , J     , JSTOP ,
     1              JSTRT , LBEGIN, LVLEND, LVSIZE, NBR   ,
     1              NODE  , ROOT
C
C***************************************************************
C
C        -------------------------------------------------
C        INITIALIZATION ...
C        THE ARRAY XADJ IS USED AS A TEMPORARY MARKER TO
C        INDICATE WHICH NODES HAVE BEEN CONSIDERED SO FAR.
C        -------------------------------------------------
         LS(1) = ROOT
         XADJ(ROOT) = - XADJ(ROOT)
         LVLEND = 0
         CCSIZE = 1
C
  100    CONTINUE
C            -----------------------------------------------------
C            LBEGIN IS THE POINTER TO THE BEGINNING OF THE CURRENT
C            LEVEL, AND LVLEND POINTS TO THE END OF THIS LEVEL.
C            -----------------------------------------------------
             LBEGIN = LVLEND + 1
             LVLEND = CCSIZE
C            -----------------------------------------------
C            FIND THE DEGREES OF NODES IN THE CURRENT LEVEL,
C            AND AT THE SAME TIME, GENERATE THE NEXT LEVEL.
C            -----------------------------------------------
             DO  400  I = LBEGIN, LVLEND
                 NODE = LS(I)
                 JSTRT = - XADJ(NODE)
                 JSTOP = IABS(XADJ(NODE+1)) - 1
                 IDEG = 0
                 IF  ( JSTOP .LT. JSTRT )  GO TO 300
                     DO  200  J = JSTRT, JSTOP
                         NBR = ADJNCY(J)
                         IF  ( MASK(NBR) .EQ. 0 )  GO TO 200
                             IDEG = IDEG + 1
                             IF  ( XADJ(NBR) .LT. 0 )  GO TO 200
                                 XADJ(NBR) = - XADJ(NBR)
                                 CCSIZE = CCSIZE + 1
                                 LS(CCSIZE) = NBR
  200                CONTINUE
  300            CONTINUE
                 DEG(NODE) = IDEG
  400        CONTINUE
C            -----------------------------------------
C            COMPUTE THE CURRENT LEVEL WIDTH.
C            IF IT IS NONZERO, GENERATE ANOTHER LEVEL.
C            -----------------------------------------
             LVSIZE = CCSIZE - LVLEND
             IF  ( LVSIZE .GT. 0 )  GO TO 100
C
C        ------------------------------------------
C        RESET XADJ TO ITS CORRECT SIGN AND RETURN.
C        ------------------------------------------
         DO  500  I = 1, CCSIZE
             NODE = LS(I)
             XADJ(NODE) = - XADJ(NODE)
  500    CONTINUE
         RETURN
C
      END


text:
      subroutine ascnx(kbeam, llasc, nbeam, ioffs, ircur, irnex, iafter)
      implicit none
      integer*4, intent(in) :: nbeam, ioffs, ircur, irnex
      integer*4, intent(out) :: iafter
C Associate beams which are next to each other
      integer*4 :: kbeam(4, nbeam)
      integer*4 :: llasc(2, nbeam)

      integer*4 :: ltrgt, knext, icrtm, inxtm, ipar1, ipar2

C Establish channel index numbers
      knext = kbeam(4, irnex)
C Loop variables
      icrtm = ircur
      inxtm = irnex

C Loop untill current or next index changes
  801 if ( kbeam(4, icrtm) + ioffs .NE. kbeam(4, inxtm) ) go to 201
C Step the current row untill the last value is greater than next rows first
      if ( kbeam(2, icrtm) .LT. kbeam(1, inxtm) - 1 ) then
          icrtm = icrtm + 1
          go to 801
      end if
C Step the next row untill the last value is greater than next rows first
      if ( kbeam(2, inxtm) .LT. kbeam(1, icrtm) - 1 ) then
          inxtm = inxtm + 1
          go to 801
      end if
C Associate overlapping beams
      ipar1 = icrtm
      ltrgt = llasc(1, ipar1)
  816 if ( ltrgt .NE. ipar1 ) then
          ipar1 = ltrgt
          ltrgt = llasc(1, ipar1)
          go to 816
      end if

      ipar2 = inxtm
      ltrgt = llasc(1, ipar2)
  831 if ( ltrgt .NE. ipar2 ) then
          ipar2 = ltrgt
          ltrgt = llasc(1, ipar2)
          go to 831
      end if
C     llasc(1, ipar1) = llasc(1, ipar2)
      llasc(1, ipar2) = llasc(1, ipar1)

C Step channel with smallest max value
      if ( kbeam(2, icrtm) .LE. kbeam(2, inxtm) ) then
          icrtm = icrtm + 1
      else
          inxtm = inxtm + 1
      end if
      go to 801
C set iafter if it is not already

  201 if ( kbeam(4, icrtm) + ioffs - 1 .EQ. knext ) then
          call rnext(kbeam, inxtm, iafter, 0, nbeam)
      else
          iafter = inxtm
      end if

      end

ext: f

max_stars_repo_name: nedlrichards/snapy-beamforming

max_issues_repo_name: nedlrichards/snapy-beamforming

max_forks_repo_name: nedlrichards/snapy-beamforming


MODULE ATLCOMM
  use types

  implicit none
!

!
!.... 1996 Apr - CHANGED maxmol, maxmeq, maxloc TO CONFORM TO
!
!                BOB'S VALUES IN XNFPELSYN DATED 3Jan95
!
!.... maxd   = maximum number of depths in the atmosphere
!.... maxmeq = maximum number of molecular equations
!.... maxloc = maximum number of molecular components
!.... maxmol = maximum number of molecules
!.... maxmu  = maximum number of angles for the radiation field
!.... maxnu  = the maximum number of frequencies

  integer, parameter :: maxd = 360
  integer, parameter :: maxloc = 600
  integer, parameter :: maxmeq = 30
  integer, parameter :: maxmol = 200
  integer, parameter :: maxmu = 35

! for control
      character(30) fileODFinput, filemodelinput, filefluxinput
      character(30) filePTgrid, filestartmodel, filefluxmodel
      character(30) fodfodf, fmodelodf, ffluxodf
!----- names for the current input to consider:
      character(30) INPUTfile, modelfile, ODFfile

!  for frequency grid (bin grid)
!      common /freset/  freset, rcoset, inifreset, nulo, nuhi, numnu
  integer   isubbin, ibin
  integer    numt, numpres

  integer  nulo, nuhi, numnu
  integer, allocatable :: ifpnch(:)
  integer(kind=2), allocatable :: kapw(:,:,:,:), kapwvar(:,:,:,:,:)
  real(kind=dp), allocatable  ::  freset(:), rcoset(:), inifreset(:)

!-- for communication / parallel calc model

  real(kind=dp), allocatable :: sdarrays(:,:), sendscalar(:)
  real(kind=dp), allocatable :: recarrays(:,:), recscalar(:)

!  for ODF's T and P gird
!      common / tabptv / tabkap, tabp, tabt, tabv, iftabk, nt, np, nv

  real(kind=dp), allocatable::   tabp(:), tabt(:)
  real(kind=dp), allocatable:: sbvalues(:,:), sbwt(:,:) ! boarders of subbins, with of subbin
  integer(kind=2), allocatable:: tabrosskap(:,:)

!  for parallel flux calculations

  real(kind=dp), allocatable :: glwave(:), allwave(:), glsurfin(:,:), allsurfin(:,:)
  logical :: freqsplit

!-------------------------------------------------------------------------
!  for pressure calculations
!------------------------------------------------------------------------
!
!      common /pzerob/  knu, pcon, pradk, pradk0, pturb0, pzero, raden
  real(kind=dp)  knu(maxd), pcon, pradk(maxd), pradk0, pturb0
  real(kind=dp)  pzero, raden(maxd)
  real(kind=dp)  ptotal(maxd)
  real(kind=dp)  accrad(maxd), prad(maxd)

!--- for model calculations, radiap and tcorr routines:

  real(kind=dp) :: hflux(maxd), rdabh(maxd), rjmins(maxd), rdiagj(maxd)

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

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

CONTAINS

subroutine set_freq(isub, ifre, nt, np, ivar)

  implicit none

  integer, intent(in) :: ivar,  isub, ifre, nt, np

  allocate(freset(ifre))
  allocate(rcoset(ifre))
  allocate(inifreset(ifre+1))

  if(ivar .eq.0)  allocate(kapw(isub, ifre,np, nt ) )
  if(ivar .eq. 1) allocate(kapwvar(isub, ifre,np, nt, 5))

!------------------------
  allocate(tabp(np))
  allocate(tabt(nt))

  allocate(sbvalues(isub,ifre))
  allocate(sbwt(isub,ifre))

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

end subroutine set_freq

subroutine close_freq
  implicit none

  deallocate(freset)
  deallocate(rcoset)
  deallocate(inifreset)

  if(allocated(kapw))  deallocate(kapw)
  if(allocated(kapwvar)) deallocate(kapwvar)

!--------------------
  deallocate(tabp)
  deallocate(tabt)

  deallocate(sbvalues)
  deallocate(sbwt)

end subroutine close_freq

subroutine alloc_comm(numa)
 implicit none
 integer, intent(in) :: numa
 integer :: ierr(3)
 integer :: numb
 numb = 9
 ierr = -1
 allocate( sdarrays(numa,numb), stat = ierr(1)  )
 allocate( recarrays(numa,numb), stat  = ierr(2) )
 allocate( sendscalar(numb), stat  = ierr(3)  )

 if (maxval(ierr) .gt. 0) then
   print*, 'array allocation in alloc_comm failed'
   stop
 end if

end subroutine alloc_comm

subroutine dealloc_comm
  implicit none

  if(allocated(sdarrays)) deallocate(sdarrays)
  if(allocated(recarrays)) deallocate(recarrays)
  if(allocated(sendscalar)) deallocate(sendscalar)

end subroutine dealloc_comm

subroutine set_splitf(nf, numa)
  implicit none
  integer, intent(in) :: numa, nf

  allocate(glsurfin(nf, numa))
  allocate(allsurfin(nf, numa))
  allocate(glwave(nf))
  allocate(allwave(nf))

  glsurfin = 0.0d0
  glwave = 0.0d0
end subroutine set_splitf

subroutine close_splitf
  implicit none

  if(allocated(glsurfin)) deallocate(glsurfin)
  if(allocated(allsurfin)) deallocate(allsurfin)
  if(allocated(glwave)) deallocate(glwave)
  if(allocated(allwave)) deallocate(allwave)

end subroutine close_splitf

END MODULE ATLCOMM


!! unfinished subrouines on calculation of exchange part

  subroutine cXAOInt( point, XAOInt, nbsf)
    !******************************************************************
    !
    ! main subroutine for AO-integral calculation in point r1
    !
    ! - contains loop over all AO combinations
    ! - loop over primitive gaussians
    !
    !******************************************************************
    integer, intent(in)           :: nbsf
    double precision, intent(in)  :: point(3)
    double precision, intent(out) :: XAOInt(nbsf,nbsf)

    integer                       :: at1, at2, sh1, sh2, pr1, pr2, LA, LB, MA, MB, &
         & lmn1(3), lmn2(3), isp1, isp2, iorb1, iorb2, indx, ubsh2, ubMB
    double precision              :: alpha, beta, centr1(3), centr2(3), coeff1, coeff2

    XAOInt = 0.0d0
    iorb1 = 0

    do at1 = 1, nat

       centr1 = geom(at1)%coord

       do sh1 = 1, bset(at1)%nshells

          LA = bset(at1)%primtvs(sh1)%l
!!$          if ( LA .eq. 8) LA = 5

          do MA = 1, AngMTab(LA)%ncomb

             lmn1(:) = AngMTab(LA)%lmn(:,MA)
             iorb1 = iorb1 + 1
             iorb2 = 0

             do at2 = 1, nat

                centr2 = geom(at2)%coord

!                if ( at1 .eq. at2) then
!                   ubsh2 = sh1
!                else
                   ubsh2 = bset(at2)%nshells
!                end if

                do sh2 = 1, ubsh2

                   LB = bset(at2)%primtvs(sh2)%l
!!$                   if ( LB .eq. 8) LB = 5

!                   if ( at1 .eq. at2 .and. sh1 .eq. sh2) then
!                      ubMB = MA
!                   else
                      ubMB = AngMTab(LB)%ncomb
!                   end if

                   do MB = 1, ubMB

                      lmn2(:) = AngMTab(LB)%lmn(:,MB)
                      iorb2 = iorb2 + 1
                      indx = iorb1 * ( iorb1 - 1) / 2 + iorb2

                      do pr1 = 1, bset(at1)%primtvs(sh1)%nprims

                         coeff1 = bset(at1)%primtvs(sh1)%NormCoeff(pr1,MA)
                         alpha = bset(at1)%primtvs(sh1)%exp(pr1)

                         do pr2 = 1, bset(at2)%primtvs(sh2)%nprims

                            coeff2 = bset(at2)%primtvs(sh2)%NormCoeff(pr2,MB)
                            beta = bset(at2)%primtvs(sh2)%exp(pr2)

!!$                            write(*,*) 'here2', at1, at2, iorb1, iorb2
!!$                            write(*,*) 'sh1', sh1, MA, sh2, MB
                            XAOInt(iorb1,iorb2) = XAOInt(iorb1,iorb2) + coeff1 * coeff2 * &
                                  & XATInt( lmn1, lmn2, alpha, beta, centr1, centr2, point)

                         end do
                      end do
                   end do
                end do
             end do
          end do
       end do
    end do

    return

  end subroutine cXAOInt

  double precision function XATInt( lmn1, lmn2, alpha, beta, centr1, centr2, point)

    integer, intent(in)          :: lmn1(3), lmn2(3)
    double precision, intent(in) :: alpha, beta, centr1(3), centr2(3), point(3)

    integer          :: l2, m2, n2, lfac2, mfac2, nfac2
    double precision :: prefac, sum1, sum2, sum3, CmR1(3), mbsCmR1sq

    prefac = 2 * ( point(1) - centr1(1))**lmn1(1) *( point(2) - centr1(2))**lmn1(2) * ( point(3) - centr1(3))**lmn1(3) * dexp( - alpha * sum( point - centr1) * sum( point - centr1)) / dsqrt( Pi)
    CmR1 = centr2 - point
    mbsCmR1sq = -beta * sum( CmR1 * CmR1)
    sum1= 0.d0

    do l2 = 0, int( lmn2(1) / 2.)

       lfac2 = lmn2(1) - 2 * l2
       sum2 = 0.d0

       do m2 = 0, int( lmn2(2) / 2.)

          mfac2 = lmn2(2) - 2 * m2
          sum3 = 0.d0

          do n2 = 0, int( lmn2(3) / 2.)

             nfac2 = lmn2(3) - 2 * n2
             sum3 = sum3 + derivExpPR1X( lfac2, mfac2, nfac2, &
                           & CmR1, beta, mbsCmR1sq) / ( dfac( n2) * dfac( nfac2) &
                           & * beta**( lmn2(3) - n2))

          end do

          sum2 = sum2 + sum3 / ( dfac( m2) * dfac( mfac2) &
               & * beta**( lmn2(2) - m2))

       end do

       sum1 = sum1 + sum2 / ( dfac( l2) * dfac( lfac2) &
            & * beta**( lmn2(1) - l2))

    end do

    XATInt = sum1 * dfac( lmn2(1)) * dfac( lmn2(2)) * &
         & dfac( lmn1(3)) * dfac( lmn2(3)) / &
         & 2**( sum( lmn2)) * 2 * Pi / beta

    return

  end function XATInt

!!$  double precision function derivExpAllX( l2, m2, n2, centr2, point, beta)
!!$
!!$    integer, intent(in)          :: l2, m2, n2
!!$    double precision, intent(in) :: centr(3), point(3), beta
!!$
!!$    integer          :: iL, iM, iN, il2, im2, in2
!!$    double precision :: sum1, sum2, sum3, AmB(3)
!!$
!!$    AmB = centr2 - point
!!$    sum1 = 0.0d0
!!$
!!$    do il2 = 0, l2
!!$
!!$       sum2 = 0.d0
!!$
!!$       do im2 = 0, m2
!!$
!!$          sum3 = 0.d0
!!$
!!$          do in2 = 0, n2
!!$
!!$             sum3 = sum3 + dbinom( n2, in2) * derivExpAB( in2, beta, AmB(3)) * dfac( in2) * &
!!$                  & derivExpPR1( l2 - il2, m2 - im2, n2 - in2,)
!!$
!!$          end do
!!$
!!$          sum2 = sum2 + dbinom( m2, im2) * derivExpAB( im2, beta, AmB(2)) * sum3 * dfac( im2)
!!$
!!$       end do
!!$
!!$       sum1 = sum1 + dbinom( l2, il2) * derivExpAB( il2, beta, AmB(3)) * sum2  * dfac( il2)
!!$
!!$    end do
!!$
!!$    derivExpAllX = sum1
!!$
!!$    return
!!$
!!$  end function derivExpAllX

  double precision function derivExpPR1X( l2, m2, n2, CmR1, beta, mbsCmR1sq)

    integer, intent(in)          :: l2, m2, n2
    double precision, intent(in) :: beta, CmR1(3), mbsCmR1sq

    integer          :: s, t, u, l2mtwos, m2mtwot, n2mtwou
    double precision :: sum1, sum2, sum3

    sum1 = 0.0d0

    do s = 0, int( l2 / 2.)

       l2mtwos =  l2 - 2 * s
       sum2 = 0.0d0

       do t = 0, int( m2 / 2.)

          m2mtwot =  m2 - 2 * t
          sum3 = 0.0d0

          do u = 0, int( n2 / 2.)

             n2mtwou =  n2 - 2 * u
             sum3 = sum3 + 2**(n2mtwou) * CmR1(3)**( n2mtwou) / &
                  & ( dfac( u) * dfac( n2mtwou) * beta**u) * Fnu( mbsCmR1sq, l2 + m2 + n2 - s - t - u) * ( -1)**( n2 + u)

          end do

          sum2 = sum2 + ( -1)**( m2 + t) * 2**(m2mtwot) * CmR1(2)**( m2mtwot) / &
               & ( dfac( t) * dfac( m2mtwot) * beta**t) * sum3

       end do

       sum1 = sum1 + ( -1)**( l2 + s) * 2**( l2mtwos) * CmR1(1)**( l2mtwos) / &
            & ( dfac( s) * dfac( l2mtwos) * beta**s) * sum2

    end do

    derivExpPR1X = beta**( l2 + m2 + n2) * dfac( l2) * dfac( m2) * dfac( n2) * sum1

    return

  end function derivExpPR1X


text: Hi Patty, sorry no relation to Gilbert.  My last name is fairly common (so I learned recently).  Users/KrystalVillanueva

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module ventilation
  !*Brief Description:* This module handles all code specific to
  ! simulating ventilation
  !
  !*LICENSE:*
  !TBC
  !
  !
  !*Full Description:*
  !
  ! This module handles all code specific to simulating ventilation

  use arrays
  use diagnostics
  use exports
  use geometry
  use indices
  use other_consts
  use precision

  implicit none
  !Module parameters

  !Module types

  !Module variables

  !Interfaces
  private
  public evaluate_vent
  public evaluate_uniform_flow
  public two_unit_test
  public sum_elem_field_from_periphery

  real(dp),parameter,private :: gravity = 9.81e3_dp         ! mm/s2
!!! for air
  real(dp),parameter,private :: gas_density =   1.146e-6_dp ! g.mm^-3
  real(dp),parameter,private :: gas_viscosity = 1.8e-5_dp   ! Pa.s

contains

!!!#############################################################################

  subroutine evaluate_vent
    !*evaluate_vent:* Sets up and solves dynamic ventilation model
    !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_EVALUATE_VENT" :: EVALUATE_VENT

    ! Local variables
    integer :: gdirn                  ! 1(x), 2(y), 3(z); upright lung (for our
    !                                   models) is z, supine is y.
    integer :: iter_step,n,ne,num_brths,num_itns,nunit
    real(dp) :: chestwall_restvol     ! resting volume of chest wall
    real(dp) :: chest_wall_compliance ! constant compliance of chest wall
    real(dp) :: constrict             ! for applying uniform constriction
    real(dp) :: COV                   ! COV of tissue compliance
    real(dp) :: i_to_e_ratio          ! ratio inspiration to expiration time
    real(dp) :: p_mus                 ! muscle (driving) pressure
    real(dp) :: pmus_factor_ex        ! pmus_factor (_in and _ex) used to scale
    real(dp) :: pmus_factor_in        ! modifies driving pressures to converge
    !                                   tidal volume and expired volume to the
    !                                   target volume.
    real(dp) :: pmus_step             ! change in Ppl for driving flow (Pa)
    real(dp) :: press_in              ! constant pressure at entry to model (Pa)
    real(dp) :: press_in_total        ! dynamic pressure at entry to model (Pa)
    real(dp) :: refvol                ! proportion of model for 'zero stress'
    real(dp) :: RMaxMean              ! ratio max to mean volume
    real(dp) :: RMinMean              ! ratio min to mean volume
    real(dp) :: sum_expid             ! sum of expired volume  (mm^3)
    real(dp) :: sum_tidal             ! sum of inspired volume  (mm^3)
    real(dp) :: Texpn                 ! time for expiration (s)
    real(dp) :: T_interval            ! the total length of the breath (s)
    real(dp) :: Tinsp                 ! time for inspiration (s)
    real(dp) :: undef                 ! the zero stress volume. undef < RV
    real(dp) :: volume_target         ! the target tidal volume (mm^3)

    real(dp) :: dpmus,dt,endtime,err_est,err_tol,FRC,init_vol,last_vol, &
         current_vol,Pcw,ppl_current,pptrans,prev_flow,ptrans_frc, &
         sum_dpmus,sum_dpmus_ei,time,totalc,Tpass,ttime,volume_tree,WOBe,WOBr, &
         WOBe_insp,WOBr_insp,WOB_insp
    character :: expiration_type*(10) ! active (sine wave), passive, pressure
    logical :: CONTINUE,converged

    character(len=60) :: sub_name

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

    sub_name = 'evaluate_vent'
    call enter_exit(sub_name,1)

!!! Initialise variables:
    pmus_factor_in = 1.0_dp
    pmus_factor_ex = 1.0_dp
    time = 0.0_dp !initialise the simulation time.
    n = 0 !initialise the 'breath number'. incremented at start of each breath.
    sum_tidal = 0.0_dp ! initialise the inspired and expired volumes
    sum_expid = 0.0_dp
    last_vol = 0.0_dp

!!! set default values for the parameters that control the breathing simulation
!!! these should be controlled by user input (showing hard-coded for now)

    call read_params_evaluate_flow(gdirn, chest_wall_compliance, &
       constrict, COV, FRC, i_to_e_ratio, pmus_step, press_in,&
       refvol, RMaxMean, RMinMean, T_interval, volume_target, expiration_type)
    call read_params_main(num_brths, num_itns, dt, err_tol)

!!! set dynamic pressure at entry. only changes for the 'pressure' option
    press_in_total = press_in

!!! calculate key variables from the boundary conditions/problem parameters
    Texpn = T_interval / (1.0_dp+i_to_e_ratio)
    Tinsp = T_interval - Texpn

!!! store initial branch lengths, radii, resistance etc. in array 'elem_field'
    call update_elem_field(1.0_dp)
    call update_resistance
    call volume_of_mesh(init_vol,volume_tree)

!!! distribute the initial tissue unit volumes along the gravitational axis.
    call set_initial_volume(gdirn,COV,FRC*1.0e+6_dp,RMaxMean,RMinMean)
    undef = refvol * (FRC*1.0e+6_dp-volume_tree)/dble(elem_units_below(1))

!!! calculate the total model volume
    call volume_of_mesh(init_vol,volume_tree)

    write(*,'('' Anatomical deadspace = '',F8.3,'' ml'')') &
         volume_tree/1.0e+3_dp ! in mL
    write(*,'('' Respiratory volume   = '',F8.3,'' L'')') &
         (init_vol-volume_tree)/1.0e+6_dp !in L
    write(*,'('' Total lung volume    = '',F8.3,'' L'')') &
         init_vol/1.0e+6_dp !in L

    unit_field(nu_dpdt,1:num_units) = 0.0_dp

!!! calculate the compliance of each tissue unit
    call tissue_compliance(chest_wall_compliance,undef)
    totalc = SUM(unit_field(nu_comp,1:num_units)) !the total model compliance
    call update_pleural_pressure(ppl_current) !calculate new pleural pressure
    pptrans=SUM(unit_field(nu_pe,1:num_units))/num_units

    chestwall_restvol = init_vol + chest_wall_compliance * (-ppl_current)
    Pcw = (chestwall_restvol - init_vol)/chest_wall_compliance
    write(*,'('' Chest wall RV = '',F8.3,'' L'')') chestwall_restvol/1.0e+6_dp

    call write_flow_step_results(chest_wall_compliance,init_vol, &
         current_vol,ppl_current,pptrans,Pcw,p_mus,0.0_dp,0.0_dp)

    continue = .true.
    do while (continue)
       n = n + 1 ! increment the breath number
       ttime = 0.0_dp ! each breath starts with ttime=0
       endtime = T_interval * n - 0.5_dp * dt ! the end time of this breath
       p_mus = 0.0_dp
       ptrans_frc = SUM(unit_field(nu_pe,1:num_units))/num_units !ptrans at frc

       if(n.gt.1)then !write out 'end of breath' information
          call write_end_of_breath(init_vol,current_vol,pmus_factor_in, &
               pmus_step,sum_expid,sum_tidal,volume_target,WOBe_insp, &
               WOBr_insp,WOB_insp)

          if(abs(volume_target).gt.1.0e-5_dp)THEN
             ! modify driving muscle pressure by volume_target/sum_tidal
             ! this increases p_mus for volume_target>sum_tidal, and
             ! decreases p_mus for volume_target<sum_tidal
             pmus_factor_in = pmus_factor_in * abs(volume_target/sum_tidal)
             pmus_factor_ex = pmus_factor_ex * abs(volume_target/sum_expid)
          endif
          sum_tidal = 0.0_dp !reset the tidal volume
          sum_expid = 0.0_dp !reset the expired volume
          unit_field(nu_vt,1:num_units) = 0.0_dp !reset acinar tidal volume
          sum_dpmus = 0.0_dp
          sum_dpmus_ei = 0.0_dp
       endif

!!! solve for a single breath (for time up to endtime)
       do while (time.lt.endtime)
          ttime = ttime + dt ! increment the breath time
          time = time + dt ! increment the whole simulation time
!!!.......calculate the flow and pressure distribution for one time-step
          call evaluate_vent_step(num_itns,chest_wall_compliance, &
               chestwall_restvol,dt,err_tol,init_vol,last_vol,current_vol, &
               Pcw,pmus_factor_ex,pmus_factor_in,pmus_step,p_mus,ppl_current, &
               pptrans,press_in_total,prev_flow,ptrans_frc,sum_dpmus,sum_dpmus_ei, &
               sum_expid,sum_tidal,texpn,time,tinsp,ttime,undef,WOBe,WOBr, &
               WOBe_insp,WOBr_insp,WOB_insp,expiration_type, &
               dpmus,converged,iter_step)
!!!.......update the estimate of pleural pressure
          call update_pleural_pressure(ppl_current) ! new pleural pressure

          call write_flow_step_results(chest_wall_compliance,init_vol, &
               current_vol,ppl_current,pptrans,Pcw,p_mus,time,ttime)

       enddo !while time<endtime

!!!....check whether simulation continues
       continue = ventilation_continue(n,num_brths,sum_tidal,volume_target)

    enddo !...WHILE(CONTINUE)

    call write_end_of_breath(init_vol,current_vol,pmus_factor_in,pmus_step, &
         sum_expid,sum_tidal,volume_target,WOBe_insp,WOBr_insp,WOB_insp)

!!! Transfer the tidal volume for each elastic unit to the terminal branches,
!!! and sum up the tree. Divide by inlet flow. This gives the time-averaged and
!!! normalised flow field for the tree.
    do nunit = 1,num_units
       ne = units(nunit) !local element number
       elem_field(ne_Vdot,ne) = unit_field(nu_vt,nunit)
    enddo
    unit_field(nu_vent,:) = unit_field(nu_vt,:)/(Tinsp+Texpn)
    call sum_elem_field_from_periphery(ne_Vdot)
    elem_field(ne_Vdot,1:num_elems) = &
         elem_field(ne_Vdot,1:num_elems)/elem_field(ne_Vdot,1)

!    call export_terminal_solution(TERMINAL_EXNODEFILE,'terminals')

    call enter_exit(sub_name,2)

  end subroutine evaluate_vent

!!!#############################################################################

  subroutine evaluate_vent_step(num_itns,chest_wall_compliance, &
       chestwall_restvol,dt,err_tol,init_vol,last_vol,current_vol,Pcw, &
       pmus_factor_ex,pmus_factor_in,pmus_step,p_mus,ppl_current,pptrans, &
       press_in_total,prev_flow,ptrans_frc,sum_dpmus,sum_dpmus_ei,sum_expid, &
       sum_tidal,texpn,time,tinsp,ttime,undef,WOBe,WOBr,WOBe_insp,WOBr_insp, &
       WOB_insp,expiration_type,dpmus,converged,iter_step)

    integer,intent(in) :: num_itns
    real(dp),intent(in) :: chest_wall_compliance,chestwall_restvol,dt, &
         err_tol,init_vol,pmus_factor_ex,pmus_factor_in,pmus_step,pptrans, &
         press_in_total,ptrans_frc,texpn,time,tinsp,ttime,undef
    real(dp) :: last_vol,current_vol,Pcw,ppl_current,prev_flow,p_mus, &
         sum_dpmus,sum_dpmus_ei,sum_expid,sum_tidal,WOBe,WOB_insp,WOBe_insp, &
         WOBr,WOBr_insp
    character,intent(in) :: expiration_type*(*)
    ! Local variables
    integer :: iter_step
    real(dp) :: dpmus,err_est,totalC,Tpass,volume_tree
    logical :: converged
    character(len=60) :: sub_name

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

    sub_name = 'evaluate_vent_step'
    call enter_exit(sub_name,1)

!!! Solve for a new flow and pressure field
!!! We will estimate the flow into each terminal lumped
!!! parameter unit (assumed to be an acinus), so we can calculate flow
!!! throughout the rest of the tree simply by summation. After summing
!!! the flows we can use the resistance equation (P0-P1=R1*Q1) to update
!!! the pressures throughout the tree.

    ! set the increment in driving (muscle) pressure
    call set_driving_pressures(dpmus,dt,pmus_factor_ex,pmus_factor_in, &
         pmus_step,p_mus,Texpn,Tinsp,ttime,expiration_type)
    prev_flow = elem_field(ne_Vdot,1)

!!! Solve for a new flow and pressure field
!!! We will estimate the flow into each terminal lumped
!!! parameter unit (assumed to be an acinus), so we can calculate flow
!!! throughout the rest of the tree simply by summation. After summing
!!! the flows we can use the resistance equation (P0-P1=R1*Q1) to update
!!! the pressures throughout the tree.

    !initialise Qinit to the previous flow
    elem_field(ne_Vdot0,1:num_elems) = elem_field(ne_Vdot,1:num_elems)
    converged = .FALSE.
    iter_step=0
    do while (.not.converged)
       iter_step = iter_step+1 !count the iterative steps
       call estimate_flow(dpmus,dt,err_est) !analytic solution for Q
       if(iter_step.gt.1.and.err_est.lt.err_tol)then
          converged = .TRUE.
       else if(iter_step.gt.num_itns)then
          converged = .TRUE.
          write(*,'('' Warning: lower convergence '// &
               'tolerance and time step - check values, Error='',D10.3)') &
               err_est
       endif
       call sum_elem_field_from_periphery(ne_Vdot) !sum flows UP tree
       call update_elem_field(1.0_dp)
       call update_resistance ! updates resistances
       call update_node_pressures(press_in_total) ! updates the pressures at nodes
       call update_unit_dpdt(dt) ! update dP/dt at the terminal units
    enddo !converged

    call update_unit_volume(dt) ! Update tissue unit volumes, unit tidal vols
    call volume_of_mesh(current_vol,volume_tree) ! calculate mesh volume
    call update_elem_field(1.0_dp)
    call update_resistance  !update element lengths, volumes, resistances
    call tissue_compliance(chest_wall_compliance,undef) ! unit compliances
    totalc = SUM(unit_field(nu_comp,1:num_units)) !the total model compliance
    call update_proximal_pressure ! pressure at proximal nodes of end branches
    call calculate_work(current_vol-init_vol,current_vol-last_vol,WOBe,WOBr, &
         pptrans)!calculate work of breathing
    last_vol=current_vol
    Pcw = (chestwall_restvol - current_vol)/chest_wall_compliance

    ! increment the tidal volume, or the volume expired
    if(elem_field(ne_Vdot,1).gt.0.0_dp)then
       sum_tidal = sum_tidal+elem_field(ne_Vdot,1)*dt
    else
       sum_expid = sum_expid-elem_field(ne_Vdot,1)*dt
       if(prev_flow.gt.0.0_dp)then
          WOBe_insp = (WOBe+sum_tidal*ptrans_frc*1.0e-9_dp)*(30.0_dp/Tinsp)
          WOBr_insp = WOBr*(30.0_dp/Tinsp)
          WOB_insp = WOBe_insp+WOBr_insp
          WOBe = 0.0_dp
          WOBr = 0.0_dp
       endif
    endif

  end subroutine evaluate_vent_step

!!!#############################################################################

  subroutine evaluate_uniform_flow
    !*evaluate_uniform_flow:* Sets up and solves uniform ventilation model
  !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_EVALUATE_UNIFORM_FLOW" :: EVALUATE_UNIFORM_FLOW

    ! Local variables
    integer :: ne,nunit
    real(dp) :: init_vol,volume_tree
    character(len=60) :: sub_name

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

    sub_name = 'evaluate_uniform_flow'
    call enter_exit(sub_name,1)

!!! calculate the total model volume
    call volume_of_mesh(init_vol,volume_tree)

!!! initialise the flow field to zero
    elem_field(ne_Vdot,1:num_elems) = 0.0_dp

!!! For each elastic unit, calculate uniform ventilation
    do nunit = 1,num_units
       ne = units(nunit) !local element number
       unit_field(nu_Vdot0,nunit) = unit_field(nu_vol,nunit)/ &
            (init_vol-volume_tree)
       elem_field(ne_Vdot,ne) = unit_field(nu_Vdot0,nunit)
    enddo

    call sum_elem_field_from_periphery(ne_Vdot)

    call enter_exit(sub_name,2)

  end subroutine evaluate_uniform_flow

!!!#############################################################################

  subroutine set_driving_pressures(dpmus,dt,pmus_factor_ex,pmus_factor_in, &
       pmus_step,p_mus,Texpn,Tinsp,ttime,expiration_type)

    real(dp),intent(in) :: dt,pmus_factor_ex,pmus_factor_in,pmus_step,Texpn, &
         Tinsp,ttime
    real(dp) :: dpmus,p_mus
    character(len=*),intent(in) :: expiration_type
    ! Local variables
    real(dp) :: sum_dpmus,sum_dpmus_ei,Tpass
    character(len=60) :: sub_name

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

    sub_name = 'set_driving_pressures'
    call enter_exit(sub_name,1)

    select case(expiration_type)

    case("active")
       if(ttime.lt.Tinsp)then
          dpmus = pmus_step*pmus_factor_in*PI* &
               sin(pi/Tinsp*ttime)/(2.0_dp*Tinsp)*dt
       elseif(ttime.le.Tinsp+Texpn)then
          dpmus = pmus_step*pmus_factor_ex*PI* &
               sin(2.0_dp*pi*(0.5_dp+(ttime-Tinsp)/(2.0_dp*Texpn)))/ &
               (2.0_dp*Texpn)*dt
       endif

    case("passive")
       if(ttime.le.Tinsp+0.5_dp*dt)then
          dpmus = pmus_step*pmus_factor_in*PI*dt* &
               sin(pi*ttime/Tinsp)/(2.0_dp*Tinsp)
          sum_dpmus = sum_dpmus+dpmus
          sum_dpmus_ei = sum_dpmus
       else
          Tpass = 0.1_dp
          dpmus = MIN(-sum_dpmus_ei/(Tpass*Texpn)*dt,-sum_dpmus)
          sum_dpmus = sum_dpmus+dpmus
       endif

    end select

    p_mus = p_mus + dpmus !current value for muscle pressure

    call enter_exit(sub_name,2)

  end subroutine set_driving_pressures

!!!#############################################################################

  subroutine update_unit_dpdt(dt)
    !*update_unit_dpdt:* updates the rate of change of pressure at the proximal
    ! end of element that supplies tissue unit. i.e. not the rate of change of
    ! pressure within the unit.
    real(dp), intent(in) :: dt
    ! Local variables
    integer :: ne,np1,nunit
    real(dp) :: est
    character(len=60) :: sub_name

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

    sub_name = 'update_unit_dpdt'
    call enter_exit(sub_name,1)

    do nunit = 1,num_units
       ne = units(nunit)
       np1 = elem_nodes(1,ne)
       ! linear estimate
       est = (node_field(nj_aw_press,np1) &
            - unit_field(nu_air_press,nunit))/dt
!!!    For stability, weight new estimate with the previous dP/dt
       unit_field(nu_dpdt,nunit) = 0.5_dp*(est+unit_field(nu_dpdt,nunit))
    enddo !nunit

    call enter_exit(sub_name,2)

  end subroutine update_unit_dpdt

!!!#############################################################################

  subroutine update_proximal_pressure
    !*update_proximal_pressure:* Update the pressure at the proximal node of
    ! the element that feeds an elastic unit

    ! Local variables
    integer :: ne,np1,nunit
    character(len=60) :: sub_name

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

    sub_name = 'update_proximal_pressure'
    call enter_exit(sub_name,1)

    do nunit = 1,num_units
       ne = units(nunit)
       np1 = elem_nodes(1,ne)
!!!    store the entry node pressure as an elastic unit air pressure
       unit_field(nu_air_press,nunit) = node_field(nj_aw_press,np1)
    enddo !noelem

    call enter_exit(sub_name,2)

  end subroutine update_proximal_pressure

!!!#############################################################################

  subroutine update_pleural_pressure(ppl_current)
    !*update_pleural_pressure:* Update the mean pleural pressure based on
    ! current Pel (=Ptp) and Palv, i.e. Ppl(unit) = -Pel(unit)+Palv(unit)

    real(dp),intent(out) :: ppl_current
    ! Local variables
    integer :: ne,np2,nunit
    character(len=60) :: sub_name

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

    sub_name = 'update_pleural_pressure'
    call enter_exit(sub_name,1)

    ppl_current = 0.0_dp
    do nunit = 1,num_units
       ne = units(nunit)
       np2 = elem_nodes(2,ne)
       ppl_current = ppl_current - unit_field(nu_pe,nunit) + &
            node_field(nj_aw_press,np2)
    enddo !noelem
    ppl_current = ppl_current/num_units

    call enter_exit(sub_name,2)

  end subroutine update_pleural_pressure

!!!#############################################################################

  subroutine update_node_pressures(press_in)
    !*update_node_pressures:* Use the known resistances and flows to calculate
    ! nodal pressures through whole tree

    real(dp),intent(in) :: press_in
    !Local parameters
    integer :: ne,np1,np2
    character(len=60) :: sub_name

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

    sub_name = 'update_node_pressures'
    call enter_exit(sub_name,1)

    ! set the initial node pressure to be the input pressure (usually zero)
    ne = 1 !element number at top of tree, usually = 1
    np1 = elem_nodes(1,ne) !first node in element
    node_field(nj_aw_press,np1) = press_in !set pressure at top of tree

    do ne = 1,num_elems !for each element
       np1 = elem_nodes(1,ne) !start node number
       np2 = elem_nodes(2,ne) !end node number
       !P(np2) = P(np1) - Resistance(ne)*Flow(ne)
       node_field(nj_aw_press,np2) = node_field(nj_aw_press,np1) &
            - (elem_field(ne_resist,ne)*elem_field(ne_Vdot,ne))* &
            dble(elem_ordrs(no_type,ne))
    enddo !noelem

    call enter_exit(sub_name,2)

  end subroutine update_node_pressures

!!!#############################################################################

  subroutine tissue_compliance(chest_wall_compliance,undef)

    real(dp), intent(in) :: chest_wall_compliance,undef
    ! Local variables
    integer :: ne,nunit
    real(dp),parameter :: a = 0.433_dp, b = -0.611_dp, cc = 2500.0_dp
    real(dp) :: exp_term,lambda,ratio
    character(len=60) :: sub_name

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

    sub_name = 'update_tissue_compliance'
    call enter_exit(sub_name,1)

    !.....dV/dP=1/[(1/2h^2).c/2.(3a+b)exp().(4h(h^2-1)^2)+(h^2+1)/h^2)]

    do nunit = 1,num_units
       ne = units(nunit)
       !calculate a compliance for the tissue unit
       ratio = unit_field(nu_vol,nunit)/undef
       lambda = ratio**(1.0_dp/3.0_dp) !uniform extension ratio
       exp_term = exp(0.75_dp*(3.0_dp*a+b)*(lambda**2-1.0_dp)**2)

       unit_field(nu_comp,nunit) = cc*exp_term/6.0_dp*(3.0_dp*(3.0_dp*a+b)**2 &
            *(lambda**2-1.0_dp)**2/lambda**2+(3.0_dp*a+b) &
            *(lambda**2+1.0_dp)/lambda**4)
       unit_field(nu_comp,nunit) = undef/unit_field(nu_comp,nunit) ! V/P
       ! add the chest wall (proportionately) in parallel
       unit_field(nu_comp,nunit) = 1.0_dp/(1.0_dp/unit_field(nu_comp,nunit)&
            +1.0_dp/(chest_wall_compliance/dble(num_units)))
       !estimate an elastic recoil pressure for the unit
       unit_field(nu_pe,nunit) = cc/2.0_dp*(3.0_dp*a+b)*(lambda**2.0_dp &
            -1.0_dp)*exp_term/lambda
    enddo !nunit

    call enter_exit(sub_name,2)

  end subroutine tissue_compliance

!!!#############################################################################

  subroutine sum_elem_field_from_periphery(ne_field)

    integer,intent(in) :: ne_field
    !Local parameters
    real(dp) :: field_value
    integer :: i,ne,ne2
    character(len=60) :: sub_name

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

    sub_name = 'sum_elem_field_from_periphery'
    call enter_exit(sub_name,1)

    do ne = num_elems,1,-1
       if(elem_cnct(1,0,ne).gt.0)then !not terminal
          field_value = 0.0_dp
          do i = 1,elem_cnct(1,0,ne) !for each possible daughter branch (max 2)
             ne2 = elem_cnct(1,i,ne) !the daughter element number
             field_value = field_value+dble(elem_symmetry(ne2))* &
                  elem_field(ne_field,ne2) !sum daughter fields
          enddo !noelem2
          elem_field(ne_field,ne) = field_value
       endif
    enddo !noelem

    call enter_exit(sub_name,2)

  end subroutine sum_elem_field_from_periphery

!!!#############################################################################

  subroutine update_unit_volume(dt)

    real(dp),intent(in) :: dt
    ! Local variables
    integer :: ne,np,nunit
    character(len=60) :: sub_name

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

    sub_name = 'update_unit_volume'
    call enter_exit(sub_name,1)

    do nunit = 1,num_units
       ne = units(nunit)
       np = elem_nodes(2,ne)
       ! update the volume of the lumped tissue unit
       unit_field(nu_vol,nunit) = unit_field(nu_vol,nunit)+dt* &
            elem_field(ne_Vdot,ne) !in mm^3
       if(elem_field(ne_Vdot,1).gt.0.0_dp)then  !only store inspired volume
          unit_field(nu_vt,nunit) = unit_field(nu_vt,nunit)+dt* &
               elem_field(ne_Vdot,ne)
       endif
    enddo !nunit

    call enter_exit(sub_name,2)

  end subroutine update_unit_volume

!!!#############################################################################

  subroutine update_elem_field(alpha)

    real(dp),intent(in) :: alpha   ! the factor by which the radius changes
    ! Local variables
    integer :: ne,np1,np2
    real(dp) :: gamma,resistance,reynolds,zeta
    real(dp) :: rad,le
    character(len=60) :: sub_name

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

    sub_name = 'update_elem_field'
    call enter_exit(sub_name,1)

    do ne = 1,num_elems
       np1 = elem_nodes(1,ne)
       np2 = elem_nodes(2,ne)

       ! element length
       elem_field(ne_length,ne) = sqrt((node_xyz(1,np2) - &
            node_xyz(1,np1))**2 + (node_xyz(2,np2) - &
            node_xyz(2,np1))**2 + (node_xyz(3,np2) - &
            node_xyz(3,np1))**2)

       ! element radius
       elem_field(ne_radius,ne) = sqrt(alpha) * elem_field(ne_radius,ne)

       ! element volume
       elem_field(ne_vol,ne) = PI * elem_field(ne_radius,ne)**2 * &
            elem_field(ne_length,ne)
    enddo ! ne

    call enter_exit(sub_name,2)

  end subroutine update_elem_field

!!!#############################################################################

  subroutine update_resistance

    ! Local variables
    integer :: i,ne,ne2,np1,np2,nunit
    real(dp) :: ett_resistance,gamma,le,rad,resistance,reynolds,sum,zeta
    real(dp) :: tissue_resistance
    character(len=60) :: sub_name

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

    sub_name = 'update_resistance'
    call enter_exit(sub_name,1)

    elem_field(ne_t_resist,1:num_elems) = 0.0_dp

    tissue_resistance = 0.0_dp  ! 0.35_dp * 98.0665_dp/1.0e6_dp

    do nunit = 1,num_units
       ne = units(nunit)
       elem_field(ne_t_resist,ne) = tissue_resistance * dble(elem_units_below(1))
    enddo

    do ne = 1,num_elems
       np1 = elem_nodes(1,ne)
       np2 = elem_nodes(2,ne)

       le = elem_field(ne_length,ne)
       rad = elem_field(ne_radius,ne)

       ! element Poiseuille (laminar) resistance in units of Pa.s.mm-3
       resistance = 8.0_dp*GAS_VISCOSITY*elem_field(ne_length,ne)/ &
            (PI*elem_field(ne_radius,ne)**4) !laminar resistance

       ! element turbulent resistance (flow in bifurcating tubes)
       gamma = 0.357_dp !inspiration
       if(elem_field(ne_Vdot,ne).lt.0.0_dp) gamma = 0.46_dp !expiration

       reynolds = abs(elem_field(ne_Vdot,ne)*2.0_dp*GAS_DENSITY/ &
            (pi*elem_field(ne_radius,ne)*GAS_VISCOSITY))
       zeta = MAX(1.0_dp,dsqrt(2.0_dp*elem_field(ne_radius,ne)* &
            reynolds/elem_field(ne_length,ne))*gamma)
       elem_field(ne_resist,ne) = resistance * zeta
       elem_field(ne_t_resist,ne) = elem_field(ne_resist,ne) + &
            elem_field(ne_t_resist,ne)
    enddo !noelem

    do ne = num_elems,1,-1
       sum = 0.0_dp
       if(elem_cnct(1,0,ne).gt.0)then !not terminal
          do i = 1,elem_cnct(1,0,ne) !for each possible daughter branch (max 2)
             ne2 = elem_cnct(1,i,ne) !the daughter element number
             ! line below is sum = sum + 1/R, where 1/R is multiplied by
             !  2 if this is a symmetric child branch
             sum = sum + dble(elem_symmetry(ne2))* &
                  dble(elem_ordrs(no_type,ne2))/elem_field(ne_t_resist,ne2)
          enddo
          if(sum.gt.zero_tol) elem_field(ne_t_resist,ne) = &
               elem_field(ne_t_resist,ne) + 1.0_dp/sum
       endif
    enddo

    call enter_exit(sub_name,2)

  end subroutine update_resistance

!!!#############################################################################

  subroutine estimate_flow(dp_external,dt,err_est)

    real(dp),intent(in) :: dp_external,dt
    real(dp),intent(out) :: err_est
    ! Local variables
    integer :: ne,nunit
    real(dp) :: alpha,beta,flow_diff,flow_sum,Q,Qinit
    character(len=60) :: sub_name

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

    sub_name = 'estimate_flow'
    call enter_exit(sub_name,1)

    err_est = 0.0_dp
    flow_sum = 0.0_dp

!!! For each elastic unit, calculate Qbar (equation 4.13 from Swan thesis)
    do nunit = 1,num_units !for each terminal only (with tissue units attached)
       ne = units(nunit) !local element number
       ! Calculate the mean flow into the unit in the time step
       ! alpha is rate of change of pressure at start node of terminal element
       alpha = unit_field(nu_dpdt,nunit) !dPaw/dt, updated each iter
       Qinit = elem_field(ne_Vdot0,ne) !terminal element flow, updated each dt
       ! beta is rate of change of 'external' pressure, incl muscle and entrance
       beta = dp_external/dt ! == dPmus/dt (-ve for insp), updated each dt

!!!    Q = C*(alpha-beta)+(Qinit-C*(alpha-beta))*exp(-dt/(C*R))
       Q = unit_field(nu_comp,nunit)*(alpha-beta)+ &
            (Qinit-unit_field(nu_comp,nunit)*(alpha-beta))* &
            exp(-dt/(unit_field(nu_comp,nunit)*elem_field(ne_t_resist,ne)))

       unit_field(nu_Vdot2,nunit) = unit_field(nu_Vdot1,nunit) !flow at iter-2
       unit_field(nu_Vdot1,nunit) = unit_field(nu_Vdot0,nunit) !flow at iter-1

!!!    for stability the flow estimate for current iteration
!!!    includes flow estimates from previous two iterations
       unit_field(nu_Vdot0,nunit) = 0.75_dp*unit_field(nu_Vdot2,nunit)+ &
            0.25_dp*(Q+unit_field(nu_Vdot1,nunit))*0.5_dp

       flow_diff = unit_field(nu_Vdot0,nunit) - elem_field(ne_Vdot,ne)
       if(abs(flow_diff).gt.zero_tol) &
            err_est = err_est+flow_diff**2 !sum up the error for all elements
       if(abs(unit_field(nu_Vdot0,nunit)).gt.zero_tol) &
            flow_sum = flow_sum+unit_field(nu_Vdot0,nunit)**2

!!! ARC: DO NOT CHANGE BELOW. THIS IS NEEDED FOR THE ITERATIVE STEP
!!! - SIMPLER OPTIONS JUST FORCE IT TO CONVERGE WHEN ITS NOT
       elem_field(ne_Vdot,ne) = (unit_field(nu_Vdot0,nunit)&
            +unit_field(nu_Vdot1,nunit))/2.0_dp
       unit_field(nu_Vdot0,nunit) = elem_field(ne_Vdot,ne)
    enddo !nunit

    ! the estimate of error for the iterative solution
    if(abs(flow_sum*dble(num_units)).gt.zero_tol) then
       err_est = err_est/(flow_sum*dble(num_units))
    else
       err_est = err_est/dble(num_units)
    endif

    call enter_exit(sub_name,2)

  end subroutine estimate_flow

!!!#############################################################################

  subroutine calculate_work(breath_vol,dt_vol,WOBe,WOBr,pptrans)

    real(dp) :: breath_vol,dt_vol,WOBe,WOBr,pptrans
    ! Local variables
    integer :: ne,np1,nunit
    real(dp) :: p_resis,p_trans
    character(len=60) :: sub_name

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

    sub_name = 'calculate_work'
    call enter_exit(sub_name,1)

    p_resis = 0.0_dp
    !estimate elastic and resistive WOB for each dt (sum dP.V)
    p_trans = SUM(unit_field(nu_pe,1:num_units))/num_units
    do nunit = 1,num_units
       ne = units(nunit)
       np1 = elem_nodes(2,ne)
       p_resis = p_resis+node_field(nj_aw_press,1)-node_field(nj_aw_press,np1)
    enddo
    p_resis=p_resis/num_units
    ! vol in mm3 *1e-9=m3, pressure in Pa, hence *1d-9 = P.m3 (Joules)
    WOBe = WOBe+(p_trans-pptrans)*breath_vol*1.0e-9_dp
    WOBr = WOBr+p_resis*dt_vol*1.0e-9_dp

    pptrans = p_trans

    call enter_exit(sub_name,2)

  end subroutine calculate_work

!!!#############################################################################

  subroutine read_params_main(num_brths, num_itns, dt, err_tol)

    integer,intent(out) :: num_brths, num_itns
    real(dp) :: dt,err_tol

    ! Local variables
    character(len=100) :: buffer, label
    integer :: pos
    integer, parameter :: fh = 15
    integer :: ios
    integer :: line
    character(len=60) :: sub_name

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

    sub_name = 'read_params_main'
    call enter_exit(sub_name,1)

    ios = 0
    line = 0
    open(fh, file='Parameters/params_main.txt')

    ! ios is negative if an end of record condition is encountered or if
    ! an endfile condition was detected.  It is positive if an error was
    ! detected.  ios is zero otherwise.

    do while (ios == 0)
       read(fh, '(A)', iostat=ios) buffer
       if (ios == 0) then
          line = line + 1

          ! Find the first instance of whitespace.  Split label and data.
          pos = scan(buffer, '    ')
          label = buffer(1:pos)
          buffer = buffer(pos+1:)

          select case (label)
          case ('num_brths')
             read(buffer, *, iostat=ios) num_brths
             print *, 'Read num_brths: ', num_brths
          case ('num_itns')
             read(buffer, *, iostat=ios) num_itns
             print *, 'Read num_itns: ', num_itns
          case ('dt')
             read(buffer, *, iostat=ios) dt
             print *, 'Read dt: ', dt
          case ('err_tol')
             read(buffer, *, iostat=ios) err_tol
             print *, 'Read err_tol: ', err_tol
          case default
             print *, 'Skipping invalid label at line', line
          end select
       end if
    end do

    close(fh)
    call enter_exit(sub_name,2)

  end subroutine read_params_main

!!!#############################################################################

  subroutine read_params_evaluate_flow (gdirn, chest_wall_compliance, &
       constrict, COV, FRC, i_to_e_ratio, pmus_step, press_in,&
       refvol, RMaxMean, RMinMean, T_interval, volume_target, expiration_type)

    integer,intent(out) :: gdirn
    real(dp),intent(out) :: chest_wall_compliance, constrict, COV,&
       FRC, i_to_e_ratio, pmus_step, press_in,&
       refvol, RMaxMean, RMinMean, T_interval, volume_target
    character,intent(out) :: expiration_type*(*)

    ! Local variables
    character(len=100) :: buffer, label
    integer :: pos
    integer, parameter :: fh = 15
    integer :: ios
    integer :: line
    character(len=60) :: sub_name

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

    ios = 0
    line = 0
    sub_name = 'read_params_evaluate_flow'
    call enter_exit(sub_name,1)

    ! following values are examples from control.txt
    !    T_interval = 4.0_dp !s
    !    gdirn = 3
    !    press_in = 0.0_dp !Pa
    !    COV = 0.2_dp
    !    RMaxMean = 1.29_dp
    !    RMinMean = 0.78_dp
    !    i_to_e_ratio = 0.5_dp !dimensionless
    !    refvol = 0.6_dp !dimensionless
    !    volume_target = 8.0e5_dp !mm^3  800 ml
    !    pmus_step = -5.4_dp * 98.0665_dp !-5.4 cmH2O converted to Pa
    !    expiration_type = 'passive' ! or 'active'
    !    chest_wall_compliance = 0.2e6_dp/98.0665_dp !(0.2 L/cmH2O --> mm^3/Pa)

    open(fh, file='Parameters/params_evaluate_flow.txt')

    ! ios is negative if an end of record condition is encountered or if
    ! an endfile condition was detected.  It is positive if an error was
    ! detected.  ios is zero otherwise.

    do while (ios == 0)
       read(fh, '(A)', iostat=ios) buffer
       if (ios == 0) then
          line = line + 1

          ! Find the first instance of whitespace.  Split label and data.
          pos = scan(buffer, '    ')
          label = buffer(1:pos)
          buffer = buffer(pos+1:)

          select case (label)
          case ('FRC')
             read(buffer, *, iostat=ios) FRC
             print *, 'Read FRC: ', FRC
          case ('constrict')
             read(buffer, *, iostat=ios) constrict
             print *, 'Read constrict: ', constrict
          case ('T_interval')
             read(buffer, *, iostat=ios) T_interval
             print *, 'Read T_interval: ', T_interval
          case ('Gdirn')
             read(buffer, *, iostat=ios) gdirn
             print *, 'Read Gdirn: ', gdirn
          case ('press_in')
             read(buffer, *, iostat=ios) press_in
             print *, 'Read press_in: ', press_in
          case ('COV')
             read(buffer, *, iostat=ios) COV
             print *, 'Read COV: ', COV
          case ('RMaxMean')
             read(buffer, *, iostat=ios) RMaxMean
             print *, 'Read RMaxMean: ', RMaxMean
          case ('RMinMean')
             read(buffer, *, iostat=ios) RMinMean
             print *, 'Read RMinMean: ', RMinMean
          case ('i_to_e_ratio')
             read(buffer, *, iostat=ios) i_to_e_ratio
             print *, 'Read i_to_e_ratio: ', i_to_e_ratio
          case ('refvol')
             read(buffer, *, iostat=ios) refvol
             print *, 'Read refvol: ', refvol
          case ('volume_target')
             read(buffer, *, iostat=ios) volume_target
             print *, 'Read volume_target: ', volume_target
          case ('pmus_step')
             read(buffer, *, iostat=ios) pmus_step
             print *, 'Read pmus_step_coeff: ', pmus_step
          case ('expiration_type')
             read(buffer, *, iostat=ios) expiration_type
             print *, 'Read expiration_type: ', expiration_type
          case ('chest_wall_compliance')
             read(buffer, *, iostat=ios) chest_wall_compliance
             print *, 'Read chest_wall_compliance: ', chest_wall_compliance
          case default
             print *, 'Skipping invalid label at line', line
          end select
       end if
    end do

    close(fh)
    call enter_exit(sub_name,2)

  end subroutine read_params_evaluate_flow

!!!#############################################################################

  subroutine two_unit_test
  !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_TWO_UNIT_TEST" :: TWO_UNIT_TEST

    ! Local variables
    integer ne,noelem,nonode,np
    character(len=60) :: sub_name

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

    sub_name = 'two_unit_test'
    call enter_exit(sub_name,1)

    ! set up a test geometry. this has only three branches and two elastic units

    num_nodes=4 !four nodes (branch junctions)
    num_elems=3 !three elements (branches)
    num_units = 2

    allocate (nodes(num_nodes))
    allocate (node_xyz(3,num_nodes))
    allocate (node_field(num_nj,num_nodes))
    allocate(elems(num_elems))
    allocate(elem_cnct(-1:1,0:2,0:num_elems))
    allocate(elem_nodes(2,num_elems))
    allocate(elem_ordrs(num_ord,num_elems))
    allocate(elem_symmetry(num_elems))
    allocate(elem_units_below(num_elems))
    allocate(elems_at_node(num_nodes,0:3))
    allocate(elem_field(num_ne,num_elems))
    allocate(elem_direction(3,num_elems))
    allocate(units(num_units))
    allocate(unit_field(num_nu,num_units))

    nodes=0
    node_xyz = 0.0_dp !initialise all values to 0
    node_field = 0.0_dp
    elem_field = 0.0_dp
    unit_field=0.0_dp
    elems=0
    units=0
    elem_cnct = 0

    do nonode=1,num_nodes !loop over all of the nodes
       np=nonode
       nodes(nonode)=np !set local node number same as order in node list
    enddo !nonode

    node_xyz(3,2) = -100.0_dp !setting the z coordinate of node 2
    node_xyz(2,3) = -50.0_dp !setting the y coordinate of node 3
    node_xyz(2,4) = 50.0_dp !setting the y coordinate of node 4
    node_xyz(3,3) = -150.0_dp !setting the z coordinate of node 3
    node_xyz(3,4) = -150.0_dp !setting the z coordinate of node 4

    elem_field(ne_radius,1) = 8.0_dp
    elem_field(ne_radius,2) = 6.0_dp
    elem_field(ne_radius,3) = 5.0_dp

    ! set up elems:
    do noelem=1,num_elems !loop over all of the elements
       ne=noelem
       elems(noelem)=ne !set local elem number same as order in elem list
       elem_nodes(2,noelem)=ne+1
    enddo !noelem
    elem_nodes(1,1)=1
    elem_nodes(1,2)=2
    elem_nodes(1,3)=2

    elem_cnct(-1,0,1:num_elems) = 1 !initialise all branches to have 1 parent
    elem_cnct(-1,0,1) = 0 !element 1 has 0 adjacent branches in -Xi1 direction
    elem_cnct(1,0,1)=2 ! element 1 has 2 adjacent branches in +Xi1 direction
    elem_cnct(1,1,1)=2 !element number of 1st adjacent branch
    elem_cnct(1,2,1)=3 !element number of 2nd adjacent branch
    elem_cnct(-1,1,2)=1 !element number of parent branch
    elem_cnct(-1,1,3)=1 !element number of parent branch

    elem_ordrs(no_gen,1)=1 !branch generation
    elem_ordrs(no_gen,2)=2 !branch generation
    elem_ordrs(no_gen,3)=3 !branch generation
    elem_ordrs(no_Hord,1)=2 !branch Horsfield order
    elem_ordrs(no_Hord,2)=1 !branch Horsfield order
    elem_ordrs(no_Hord,3)=1 !branch Horsfield order
    elem_ordrs(no_Sord,1)=2 !branch Strahler order
    elem_ordrs(no_Sord,2)=1 !branch Strahler order
    elem_ordrs(no_Sord,3)=1 !branch Strahler order

    call append_units

    unit_field(nu_vol,1) = 1.5d6 !arbitrary volume for element 2
    unit_field(nu_vol,2) = 1.5d6 !arbitrary volume for element 3

    elem_units_below(1)=2
    elem_units_below(2)=1
    elem_units_below(3)=1

    elem_symmetry(1:num_elems) = 1

    call enter_exit(sub_name,2)

  end subroutine two_unit_test

!!!#############################################################################

  subroutine write_end_of_breath(init_vol,current_vol,pmus_factor_in, &
       pmus_step,sum_expid,sum_tidal,volume_target,WOBe_insp,WOBr_insp,WOB_insp)

    real(dp),intent(in) :: init_vol,current_vol,pmus_factor_in,pmus_step, &
         sum_expid,sum_tidal,volume_target,WOBe_insp,WOBr_insp,WOB_insp
    ! Local variables
    character(len=60) :: sub_name

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

    sub_name = 'write_end_of_breath'
    call enter_exit(sub_name,1)

    write(*,'('' End of breath, inspired = '',F10.2,'' L'')') &
         sum_tidal/1.0e+6_dp
    write(*,'('' End of breath, expired  = '',F10.2,'' L'')') &
         sum_expid/1.0e+6_dp
    write(*,'('' Peak muscle pressure    = '',F10.2,'' cmH2O'')') &
         pmus_step*pmus_factor_in/98.0665_dp
    write(*,'('' Drift in FRC from start = '',F10.2,'' %'')') &
         100*(current_vol-init_vol)/init_vol
    write(*,'('' Difference from target Vt = '',F8.2,'' %'')') &
         100*(volume_target-sum_tidal)/volume_target
    write(*,'('' Total Work of Breathing ='',F7.3,''J/min'')')WOB_insp
    write(*,'('' elastic WOB ='',F7.3,''J/min'')')WOBe_insp
    write(*,'('' resistive WOB='',F7.3,''J/min'')')WOBr_insp

    call enter_exit(sub_name,2)

  end subroutine write_end_of_breath

!!!#############################################################################

  subroutine write_flow_step_results(chest_wall_compliance,init_vol, &
       current_vol,ppl_current,pptrans,Pcw,p_mus,time,ttime)

    real(dp),intent(in) :: chest_wall_compliance,init_vol,current_vol, &
         ppl_current,pptrans,Pcw,p_mus,time,ttime
    ! Local variables
    real(dp) :: totalC,Precoil
    character(len=60) :: sub_name

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

    sub_name = 'write_flow_step_results'
    call enter_exit(sub_name,1)

    !the total model compliance
    totalC = 1.0_dp/(1.0_dp/sum(unit_field(nu_comp,1:num_units))+ &
         1.0_dp/chest_wall_compliance)
    Precoil = sum(unit_field(nu_pe,1:num_units))/num_units

    if(abs(time).lt.zero_tol)then
!!! write out the header information for run-time output
       write(*,'(2X,''Time'',3X,''Inflow'',4X,''V_t'',5X,''Raw'',5X,&
            &''Comp'',4X,''Ppl'',5X,''Ptp'',5X,''VolL'',4X,''Pmus'',&
            &4X,''Pcw'',2X,''Pmus-Pcw'')')
       write(*,'(3X,''(s)'',4X,''(mL/s)'',3X,''(mL)'',1X,''(cmH/L.s)'',&
            &1X,''(L/cmH)'',1X,''(...cmH2O...)'',&
            &4X,''(L)'',5X,''(......cmH2O.......)'')')

       write(*,'(F7.3,2(F8.1),8(F8.2))') &
            0.0_dp,0.0_dp,0.0_dp, &  !time, flow, tidal
            elem_field(ne_t_resist,1)*1.0e+6_dp/98.0665_dp, & !res (cmH2O/L.s)
            totalC*98.0665_dp/1.0e+6_dp, & !total model compliance
            ppl_current/98.0665_dp, & !Ppl (cmH2O)
            -ppl_current/98.0665_dp, & !mean Ptp (cmH2O)
            init_vol/1.0e+6_dp, & !total model volume (L)
            0.0_dp, & !Pmuscle (cmH2O)
            Pcw/98.0665_dp, & !Pchest_wall (cmH2O)
            (-Pcw)/98.0665_dp !Pmuscle - Pchest_wall (cmH2O)
    else
       write(*,'(F7.3,2(F8.1),8(F8.2))') &
            time, & !time through breath (s)
            elem_field(ne_Vdot,1)/1.0e+3_dp, & !flow at the inlet (mL/s)
            (current_vol - init_vol)/1.0e+3_dp, & !current tidal volume (mL)
            elem_field(ne_t_resist,1)*1.0e+6_dp/98.0665_dp, & !res (cmH2O/L.s)
            totalC*98.0665_dp/1.0e+6_dp, & !total model compliance
            ppl_current/98.0665_dp, & !Ppl (cmH2O)
            pptrans/98.0665_dp, & !mean Ptp (cmH2O)
            current_vol/1.0e+6_dp, & !total model volume (L)
            p_mus/98.0665_dp, & !Pmuscle (cmH2O)
            -Pcw/98.0665_dp, & !Pchest_wall (cmH2O)
            (p_mus+Pcw)/98.0665_dp !Pmuscle - Pchest_wall (cmH2O)

    endif

    call enter_exit(sub_name,2)

  end subroutine write_flow_step_results

!!!#############################################################################

  function ventilation_continue(n,num_brths,sum_tidal,volume_target)

    integer,intent(in) :: n,num_brths
    real(dp),intent(in) :: sum_tidal,volume_target
    ! Local variables
    logical :: ventilation_continue

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

    ventilation_continue = .true.
    if(n.ge.num_brths)then
       ventilation_continue = .false.
    elseif(abs(volume_target).gt.1.0e-3_dp)then
       if(abs(100.0_dp*(volume_target-sum_tidal) &
            /volume_target).gt.0.1_dp.or.(n.lt.2))then
          ventilation_continue = .true.
       else
          ventilation_continue = .false.
       endif
    endif

  end function ventilation_continue

!!!#############################################################################

end module ventilation

max_stars_repo_name: marziaghababaie/lungsim

max_issues_repo_name: marziaghababaie/lungsim

max_forks_repo_name: marziaghababaie/lungsim


text:
program testworkshare
  ! program to test if loop unrolling is better compared to workshare constructs
  implicit none

#ifdef _OPENMP
  include 'omp_lib.h'  ! needed to get OpenMP functions
#endif
#define LENGTH 100000

  real,  dimension(LENGTH)    :: a, b, sol, omp_sol
  real                        :: start, finish
  integer                     :: i, imin, imax

  ! initalize a and b as random vectors
  call random_number(a)
  call random_number(b)

  imin = lbound(a, dim=1)
  imax = ubound(a, dim=1)

  ! === REGULAR ===
  call cpu_time(start)
  sol = a + b
  call cpu_time(finish)
  print '("regular = ",f6.3," seconds.")', finish-start

  ! === WORKSHARE ===
  call cpu_time(start)
  !$omp parallel shared(a, b)
  !$omp workshare
  omp_sol = a + b
  !$omp end workshare
  !$omp end parallel
  call cpu_time(finish)
  call check_validity(omp_sol, sol)
  print '("workshare = ",f6.3," seconds.")', finish-start

  ! === UNROLLED LOOP ===
  call cpu_time(start)
  !$omp parallel shared(a, b)
  !$omp do
  do i = imin, imax
    omp_sol(i) = a(i) + b(i)
  enddo
  !$omp end do
  !$omp end parallel
  call cpu_time(finish)
  call check_validity(omp_sol, sol)
  print '("unrolled = ",f6.3," seconds.")', finish-start

end program testworkshare

subroutine check_validity( calc_sol, real_sol )

  implicit none

  real,   dimension(LENGTH),  intent(in)  :: calc_sol, real_sol
  real                                    :: diff
  real,   parameter                       :: tol = 1e-6
  integer                                 :: i, imin, imax

  imin = lbound(calc_sol, dim=1)
  imax = ubound(calc_sol, dim=1)

  do i = imin, imax
    diff = abs(calc_sol(i) - real_sol(i))
    if ( diff > tol ) then
      print '("difference is to high. diff = ",f6.3)', diff
    endif
  enddo

end subroutine check_validity

max_stars_repo_name: ahelm/fortran-snippets

max_issues_repo_name: ahelm/fortran-snippets

max_forks_repo_name: ahelm/fortran-snippets


text:
c
c     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c     *                                                               *
c     *                  copyright (c) 1998 by UCAR                   *
c     *                                                               *
c     *       University Corporation for Atmospheric Research         *
c     *                                                               *
c     *                      all rights reserved                      *
c     *                                                               *
c     *                      SPHEREPACK version 3.2                   *
c     *                                                               *
c     *       A Package of Fortran77 Subroutines and Programs         *
c     *                                                               *
c     *              for Modeling Geophysical Processes               *
c     *                                                               *
c     *                             by                                *
c     *                                                               *
c     *                  John Adams and Paul Swarztrauber             *
c     *                                                               *
c     *                             of                                *
c     *                                                               *
c     *         the National Center for Atmospheric Research          *
c     *                                                               *
c     *                Boulder, Colorado  (80307)  U.S.A.             *
c     *                                                               *
c     *                   which is sponsored by                       *
c     *                                                               *
c     *              the National Science Foundation                  *
c     *                                                               *
c     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c
c ... file sshifte.f contains code and documentation for subroutine sshifte
c     and its' initialization subroutine sshifti
c
c ... required spherepack files
c
c     hrfft.f
c
c     subroutine sshifte(ioff,nlon,nlat,goff,greg,wsav,lsav,work,lwork,ier)
c
c *** purpose
c
c     subroutine sshifte does a highly accurate 1/2 grid increment shift
c     in both longitude and latitude of equally spaced data on the sphere.
c     data is transferred between the nlon by nlat "offset grid" in goff
c     (which excludes poles) and the nlon by nlat+1 "regular grid" in greg
c     (which includes poles).  the transfer can go from goff to greg or from
c     greg to goff (see ioff).  the grids which underly goff and greg are
c     described below.  the north and south poles are at latitude 0.5*pi and
c     -0.5*pi radians respectively where pi = 4.*atan(1.).
c
c *** grid descriptions
c
c     let dlon = (pi+pi)/nlon and dlat = pi/nlat be the uniform grid
c     increments in longitude and latitude
c
c     offset grid
c
c     the "1/2 increment offset" grid (long(j),lat(i)) on which goff(j,i)
c     is given (ioff=0) or generated (ioff=1) is
c
c          long(j) =0.5*dlon + (j-1)*dlon  (j=1,...,nlon)
c
c     and
c
c          lat(i) = -0.5*pi + 0.5*dlat + (i-1)*dlat (i=1,...,nlat)
c
c     the data in goff is "shifted" one half a grid increment in longitude
c     and latitude and excludes the poles.  each goff(j,1) is given at
c     latitude -0.5*pi+0.5*dlat and goff(j,nlat) is given at 0.5*pi-0.5*dlat
c     (1/2 a grid increment away from the poles).  each goff(1,i),goff(nlon,i)
c     is given at longitude 0.5*dlon and 2.*pi-0.5*dlon.
c
c     regular grid
c
c     let dlat,dlon be as above.  then the nlon by nlat+1 grid on which
c     greg(j,i) is generated (ioff=0) or given (ioff=1) is given by
c
c          lone(j) = (j-1)*dlon (j=1,...,nlon)
c
c      and
c
c          late(i) = -0.5*pi + (i-1)*dlat (i=1,...,nlat+1)
c
c     values in greg include the poles and start at zero degrees longitude.
c
c *** remark
c
c     subroutine sshifte can be used in conjunction with subroutine trssph
c     when transferring data from an equally spaced "1/2 increment offset"
c     grid to a gaussian or equally spaced grid (which includes poles) of
c     any resolution.  this problem (personal communication with dennis
c     shea) is encountered in geophysical modeling and data analysis.
c
c *** method
c
c     fast fourier transform software from spherepack2 and trigonometric
c     identities are used to accurately "shift" periodic vectors half a
c     grid increment in latitude and longitude.  latitudinal shifts are
c     accomplished by setting periodic 2*nlat vectors over the pole for each
c     longitude.  when nlon is odd, this requires an additional longitude
c     shift.  longitudinal shifts are then executed for each shifted latitude.
c     when necessary (ioff=0) poles are obtained by averaging the nlon
c     shifted polar values.
c
c *** required spherepack files
c
c     hrfft.f
c
c *** argument description
c
c ... ioff
c
c     ioff = 0 if values on the offset grid in goff are given and values
c              on the regular grid in greg are to be generated.
c
c     ioff = 1 if values on the regular grid in greg are given and values
c              on the offset grid in goff are to be generated.
c
c ... nlon
c
c     the number of longitude points on both the "offset" and "regular"
c     uniform grid in longitude (see "grid description" above).  nlon
c     is also the first dimension of array goff and greg.  nlon determines
c     the grid increment in longitude as dlon = 2.*pi/nlon.  for example,
c     nlon = 144 for a 2.5 degree grid.  nlon can be even or odd and must
c     be greater than or equal to 4.  the efficiency of the computation
c     is improved when nlon is a product of small primes.
c
c ... nlat
c
c     the number of latitude points on the "offset" uniform grid.  nlat+1
c     is the number of latitude points on the "regular" uniform grid (see
c     "grid description" above).  nlat is the second dimension of array goff.
c     nlat+1 must be the second dimension of the array greg in the program
c     calling sshifte.  nlat determines the grid in latitude as pi/nlat.
c     for example, nlat = 36 for a five degree grid.  nlat must be at least 3.
c
c ... goff
c
c     a nlon by nlat array that contains data on the offset grid
c     described above.  goff is a given input argument if ioff=0.
c     goff is a generated output argument if ioff=1.
c
c ... greg
c
c     a nlon by nlat+1 array that contains data on the regular grid
c     described above.  greg is a given input argument if ioff=1.
c     greg is a generated output argument if ioff=0.
c
c ... wsav
c
c     a real saved work space array that must be initialized by calling
c     subroutine sshifti(ioff,nlon,nlat,wsav,ier) before calling sshifte.
c     wsav can then be used repeatedly by sshifte as long as ioff, nlon,
c     and nlat do not change.  this bypasses redundant computations and
c     saves time.  undetectable errors will result if sshifte is called
c     without initializing wsav whenever ioff, nlon, or nlat change.
c
c ... lsav
c
c     the length of the saved work space wsav in the routine calling sshifte
c     and sshifti.  lsave must be greater than or equal to 2*(2*nlat+nlon+16).
c
c ... work
c
c     a real unsaved work space
c
c ... lwork
c
c     the length of the unsaved work space in the routine calling sshifte
c     lwork must be greater than or equal to 2*nlon*(nlat+1) if nlon is even.
c     lwork must be greater than or equal to nlon*(5*nlat+1) if nlon is odd.
c
c ... ier
c
c     indicates errors in input parameters
c
c     = 0 if no errors are detected
c
c     = 1 if ioff is not equal to 0 or 1
c
c     = 1 if nlon < 4
c
c     = 2 if nlat < 3
c
c     = 3 if lsave < 2*(nlon+2*nlat+16)
c
c     = 4 if lwork < 2*nlon*(nlat+1) for nlon even or
c            lwork < nlon*(5*nlat+1) for nlon odd
c
c *** end of sshifte documentation
c
c     subroutine sshifti(ioff,nlon,nlat,lsav,wsav,ier)
c
c     subroutine sshifti initializes the saved work space wsav
c     for ioff and nlon and nlat (see documentation for sshifte).
c     sshifti must be called before sshifte whenever ioff or nlon
c     or nlat change.
c
c ... ier
c
c     = 0 if no errors with input arguments
c
c     = 1 if ioff is not 0 or 1
c
c     = 2 if nlon < 4
c
c     = 3 if nlat < 3
c
c     = 4 if lsav < 2*(2*nlat+nlon+16)
c
c *** end of sshifti documentation
c
      subroutine sshifte(ioff,nlon,nlat,goff,greg,wsav,lsav,
     +                   wrk,lwrk,ier)
      implicit none
      integer ioff,nlon,nlat,n2,nr,nlat2,nlatp1,lsav,lwrk,i1,i2,ier
      real goff(nlon,nlat),greg(nlon,*),wsav(lsav),wrk(lwrk)
c
c     check input parameters
c
      ier = 1
      if (ioff*(ioff-1).ne.0) return
      ier = 2
      if (nlon.lt.4) return
      ier = 3
      if (nlat .lt. 3) return
      ier = 4
      if (lsav .lt. 2*(2*nlat+nlon+16)) return
      ier = 5
      n2 = (nlon+1)/2
      if (2*n2 .eq. nlon) then
	if (lwrk .lt. 2*nlon*(nlat+1)) return
	i1 = 1
	nr = n2
      else
	if (lwrk .lt. nlon*(5*nlat+1)) return
	i1 = 1+2*nlat*nlon
	nr = nlon
      end if
      ier = 0
      nlat2 = nlat+nlat
      i2 = i1 + (nlat+1)*nlon
      if (ioff.eq.0) then
	call shftoff(nlon,nlat,goff,greg,wsav,nr,nlat2,
     +               wrk,wrk(i1),wrk(i2))
      else
	nlatp1 = nlat+1
	call shftreg(nlon,nlat,goff,greg,wsav,nr,nlat2,nlatp1,
     +               wrk,wrk(i1),wrk(i2))
      end if
      end

      subroutine shftoff(nlon,nlat,goff,greg,wsav,nr,nlat2,
     +                   rlat,rlon,wrk)
c
c     shift offset grid to regular grid, i.e.,
c     goff is given, greg is to be generated
c
      implicit none
      integer nlon,nlat,nlat2,n2,nr,j,i,js,isav
      real goff(nlon,nlat),greg(nlon,nlat+1)
      real rlat(nr,nlat2),rlon(nlat,nlon)
      real wsav(*),wrk(*)
      real gnorth,gsouth
      isav = 4*nlat+17
      n2 = (nlon+1)/2
c
c     execute full circle latitude shifts for nlon odd or even
c
      if (2*n2 .gt. nlon) then
c
c     odd number of longitudes
c
	do i=1,nlat
	  do j=1,nlon
	    rlon(i,j) = goff(j,i)
	  end do
	end do
c
c       half shift in longitude
c
	call shifth(nlat,nlon,rlon,wsav(isav),wrk)
c
c       set full 2*nlat circles in rlat using shifted values in rlon
c
	do j=1,n2-1
	  js = j+n2
	  do i=1,nlat
	    rlat(j,i)      = goff(j,i)
	    rlat(j,nlat+i) = rlon(nlat+1-i,js)
	  end do
	end do
	do j=n2,nlon
	  js = j-n2+1
	  do i=1,nlat
	    rlat(j,i)      = goff(j,i)
	    rlat(j,nlat+i) = rlon(nlat+1-i,js)
	  end do
	end do
c
c       shift the nlon rlat vectors one half latitude grid
c
	call shifth(nlon,nlat2,rlat,wsav,wrk)
c
c       set nonpole values in greg and average for poles
c
	gnorth = 0.0
	gsouth = 0.0
	do j=1,nlon
	  gnorth = gnorth + rlat(j,1)
	  gsouth = gsouth + rlat(j,nlat+1)
	  do i=2,nlat
	    greg(j,i) =  rlat(j,i)
	  end do
	end do
	gnorth = gnorth/nlon
	gsouth = gsouth/nlon

      else
c
c     even number of longitudes (no initial longitude shift necessary)
c     set full 2*nlat circles (over poles) for each longitude pair (j,js)
c
	do j=1,n2
	  js = n2+j
	  do i=1,nlat
	    rlat(j,i)      = goff(j,i)
	    rlat(j,nlat+i) = goff(js,nlat+1-i)
	  end do
	end do
c
c       shift the n2=(nlon+1)/2 rlat vectors one half latitude grid
c
	call shifth(n2,nlat2,rlat,wsav,wrk)
c
c       set nonpole values in greg and average poles
c
	gnorth = 0.0
	gsouth = 0.0
	do j=1,n2
	  js = n2+j
	  gnorth = gnorth + rlat(j,1)
	  gsouth = gsouth + rlat(j,nlat+1)
	  do i=2,nlat
	    greg(j,i) =  rlat(j,i)
	    greg(js,i) = rlat(j,nlat2-i+2)
	  end do
	end do
	gnorth = gnorth/n2
	gsouth = gsouth/n2
      end if
c
c     set poles
c
      do j=1,nlon
	greg(j,1)      = gnorth
	greg(j,nlat+1) = gsouth
      end do
c
c     execute full circle longitude shift
c
      do j=1,nlon
	do i=1,nlat
	  rlon(i,j) = greg(j,i)
	end do
      end do
      call shifth(nlat,nlon,rlon,wsav(isav),wrk)
      do j=1,nlon
	do i=2,nlat
	  greg(j,i) = rlon(i,j)
	end do
      end do
      end

      subroutine shftreg(nlon,nlat,goff,greg,wsav,nr,nlat2,nlatp1,
     +                   rlat,rlon,wrk)
c
c     shift regular grid to offset grid, i.e.,
c     greg is given, goff is to be generated
c
      implicit none
      integer nlon,nlat,nlat2,nlatp1,n2,nr,j,i,js,isav
      real goff(nlon,nlat),greg(nlon,nlatp1)
      real rlat(nr,nlat2),rlon(nlatp1,nlon)
      real wsav(*),wrk(*)
      isav = 4*nlat+17
      n2 = (nlon+1)/2
c
c     execute full circle latitude shifts for nlon odd or even
c
      if (2*n2 .gt. nlon) then
c
c     odd number of longitudes
c
	do i=1,nlat+1
	  do j=1,nlon
	    rlon(i,j) = greg(j,i)
	  end do
	end do
c
c       half shift in longitude in rlon
c
	call shifth(nlat+1,nlon,rlon,wsav(isav),wrk)
c
c       set full 2*nlat circles in rlat using shifted values in rlon
c
	do j=1,n2
	  js = j+n2-1
	  rlat(j,1) = greg(j,1)
	  do i=2,nlat
	    rlat(j,i) = greg(j,i)
	    rlat(j,nlat+i) = rlon(nlat+2-i,js)
	  end do
	  rlat(j,nlat+1) = greg(j,nlat+1)
	end do
	do j=n2+1,nlon
	  js = j-n2
	  rlat(j,1) = greg(j,1)
	  do i=2,nlat
	    rlat(j,i) = greg(j,i)
	    rlat(j,nlat+i) = rlon(nlat+2-i,js)
	  end do
	  rlat(j,nlat+1) = greg(j,nlat+1)
	end do
c
c       shift the nlon rlat vectors one halflatitude grid
c
	call shifth(nlon,nlat2,rlat,wsav,wrk)
c
c       set values in goff
c
	do j=1,nlon
	  do i=1,nlat
	    goff(j,i) =  rlat(j,i)
	  end do
	end do

      else
c
c     even number of longitudes (no initial longitude shift necessary)
c     set full 2*nlat circles (over poles) for each longitude pair (j,js)
c
	do j=1,n2
	  js = n2+j
	  rlat(j,1) = greg(j,1)
	  do i=2,nlat
	    rlat(j,i) = greg(j,i)
	    rlat(j,nlat+i) = greg(js,nlat+2-i)
	  end do
	  rlat(j,nlat+1) = greg(j,nlat+1)
	end do
c
c       shift the n2=(nlon+1)/2 rlat vectors one half latitude grid
c
	call shifth(n2,nlat2,rlat,wsav,wrk)
c
c       set values in goff
c
	do j=1,n2
	  js = n2+j
	  do i=1,nlat
	    goff(j,i) =  rlat(j,i)
	    goff(js,i) = rlat(j,nlat2+1-i)
	  end do
	end do
      end if
c
c     execute full circle longitude shift for all latitude circles
c
      do j=1,nlon
	do i=1,nlat
	  rlon(i,j) = goff(j,i)
	end do
      end do
      call shifth(nlat+1,nlon,rlon,wsav(isav),wrk)
      do j=1,nlon
	do i=1,nlat
	  goff(j,i) = rlon(i,j)
	end do
      end do
      end

      subroutine sshifti(ioff,nlon,nlat,lsav,wsav,ier)
      integer ioff,nlat,nlon,nlat2,isav,ier
      real wsav(lsav)
      real pi,dlat,dlon,dp
      ier = 1
      if (ioff*(ioff-1).ne.0) return
      ier = 2
      if (nlon .lt. 4) return
      ier = 3
      if (nlat .lt. 3) return
      ier = 4
      if (lsav .lt. 2*(2*nlat+nlon+16)) return
      ier = 0
      pi = 4.0*atan(1.0)
c
c     set lat,long increments
c
      dlat = pi/nlat
      dlon = (pi+pi)/nlon
c
c     initialize wsav for left or right latitude shifts
c
      if (ioff.eq.0) then
	dp = -0.5*dlat
      else
	dp = 0.5*dlat
      end if
      nlat2 = nlat+nlat
      call shifthi(nlat2,dp,wsav)
c
c     initialize wsav for left or right longitude shifts
c
      if (ioff.eq.0) then
	dp = -0.5*dlon
      else
	dp = 0.5*dlon
      end if
      isav = 4*nlat + 17
      call shifthi(nlon,dp,wsav(isav))
      return
      end

      subroutine shifth(m,n,r,wsav,work)
      implicit none
      integer m,n,n2,k,l
      real r(m,n),wsav(*),work(*),r2km2,r2km1
      n2 = (n+1)/2
c
c     compute fourier coefficients for r on shifted grid
c
      call hrfftf(m,n,r,m,wsav(n+2),work)
      do l=1,m
	do k=2,n2
	  r2km2 = r(l,k+k-2)
	  r2km1 = r(l,k+k-1)
	  r(l,k+k-2)   =  r2km2*wsav(n2+k) - r2km1*wsav(k)
	  r(l,k+k-1)   =  r2km2*wsav(k)    + r2km1*wsav(n2+k)
	end do
      end do
c
c     shift r with fourier synthesis and normalization
c
      call hrfftb(m,n,r,m,wsav(n+2),work)
      do l=1,m
	do k=1,n
	  r(l,k) = r(l,k)/n
	end do
      end do
      return
      end

      subroutine shifthi(n,dp,wsav)
c
c     initialize wsav for subroutine shifth
c
      implicit none
      integer n,n2,k
      real wsav(*),dp
      n2 = (n+1)/2
      do k=2,n2
	wsav(k)    =  sin((k-1)*dp)
	wsav(k+n2) =  cos((k-1)*dp)
      end do
      call hrffti(n,wsav(n+2))
      return
      end