ccc	

ccc	Title:	A validated model for the dispersion of exhaust emissions from aircraft in their take-off roll

ccc	Authors:	A. Graham1, M. Bennett*2, S. Christie2 and D.W. Raper2

ccc

ccc	* To whom correspondence should be addressed: m.bennett@mmu.ac.uk

ccc

ccc	Affliation 1:	Norwegian Centre of Offshore Wind Energy, Postboks 6031, N-5892 Bergen, Norway angus.graham@uni.no

ccc	Affiliation 2	Centre for Aviation, Transport and the Environment, Manchester Metropolitan University, Manchester M1 5GD, UK  
ccc

ccc	Copyright. Manchester Metropolitan University 2018.
ccc
ccc
ccc	The following is a FORTRAN subroutine which codes the mathematical model described in the above paper 

ccc	Specimen output files can be found in: TOdispResScenA.txt; and TOdispResScenB.txt.

ccc
ccc
ccc

ccc	The source code is F77 including the use of mathematical intrinsic functions such as hyperbolic, which most code libraries should support.

ccc	It is fully commented so as to be reasonably easily adaptable by future users. 
ccc

ccc


        real*8 fricv,wdir,invMO,F0max,MTOM,wingar,dum

        print *,'friction velocity (m s-1) ?'
        read *,fricv
        print *,'wind dir rel to aircraft hding (+- 180 degs) ?'
        read *,wdir
        print *,'inverse Monin Obukhov length (m-1) ?'
        read *,invMO
        print *,'boundary-layer height (m) ?'
        read *,dum
        iblht = idint(dum)
           if (iblht.lt.1) then
            print *,'value too low: prog ends'
            stop
           endif
        print *,'top static thrust from pooled engines (kN) ?'
        read *,dum
        F0max = 1000*dum
        print *,'maximum takeoff mass (tonnes) ?'
        read *,dum
        MTOM = 1000*dum
        print *,'basic wing planform area (m2) ?'
        read *,wingar

        call getTOd(fricv,wdir,invMO,iblht,F0max,MTOM,wingar)

        end

************************************************************************

        subroutine getTOd(fricv,wdir,invMO,iblht,F0max,MTOM,wingar)

        real*8 zero,sig,sig1,dtor,sixth,quart,third,twothd,quarpi,one,
     +         cthrhf
        parameter (zero=0,sig=1e-9,sig1=1e-3,quart=.25,one=1,
     +             sixth=one/6,third=2*sixth,twothd=2*third,cthrhf=1.5)

        real*8 kmF,kmM,rollco,gacc,corot,CLmax,rho,therc1,therc2,therc3,
     +         therc4,tau,xres,coke,cokF,corat,therc5,therc6,therc7,
     +         therc8
        parameter (kmF=.85,kmM=.90,rollco=.025,gacc=9.8,corot=1.2,
     +             CLmax=2.25,rho=1.25,therc1=11.5,therc2=59,therc3=75,
     +             therc4=5,tau=80,xres=.5,coke=.3,cokF=.7,therc5=34,
     +             therc6=third,therc7=1.5,therc8=1.4)
ccc     kmF     ..      fraction of max thrust characterisrtically
ccc                     used during takeoff
ccc     kmM     ..      fraction of max takeoff weight
ccc                     characteristically holding for departing
ccc                     aircraft
ccc     rollco  ..      coefficient of friction of aircraft tyres on
ccc                     runway
ccc     gacc    ..      free acceleration due to gravity (m s-2)
ccc     corot   ..      airspeed of aircraft when upward rotation of
ccc                     aircraft is begun, expressed as a factor of
ccc                     stall speed
ccc     CLmax   ..      maximum attainable lift coefficient of aircraft
ccc                     at takeoff flap setting
ccc     rho     ..      density of ambient air at ground level (kg m-3)
ccc     therc1  ..      log_e (z_ref / z_r) / k_vK, where z_ref is
ccc                     reference height of 5 m, z_r is representative
ccc                     ambient roughness length of 5 cm and k_vK is
ccc                     von Karman's constant, 0.4
ccc     therc2  ..      alpha_s * z_ref / k_vK (m), where alpha_s is
ccc                     4.7 (in stable ambient stratification, shear at
ccc                     height z_ref, as normalised on
ccc                     fricv / (k_vK * z_ref), is equal to
ccc                     1 + alpha_s * z_ref / L_M, L_M being Monin-
ccc                     Obukhov length)
ccc     therc3  ..      alpha_c * z_ref (m), where alpha_c is 15
ccc                     (in ambient convection, shear at height z_ref,
ccc                     as normalised on fricv / (k_vK * zref), is equal
ccc                     to 1 / (1 - alpha_c * z_ref / L_M) ** [1 / 4])
ccc     therc4  ..      2 / k_vK
ccc     tau     ..      F0 / B_F (s), where B_F is buoyant force
ccc                     imparted to exhaust pool in unit release time
ccc     xres    ..      resolution of grid in direction of release over
ccc                     which exhaust ages are solved for during ground
ccc                     run, prior to upward rotation of aircraft (m)
ccc     coke    ..      rate of change (as total derivative) of square
ccc                     root of cross-sectional area of exhausts, as
ccc                     normalised on their excess speed in direction
ccc                     of release
ccc     cokF    ..      factor by which excess momentum of exhausts in
ccc                     in direction of release is reduced as result of
ccc                     drag losses at ground
ccc     therc5  ..      value in stable ambient stratification of k_rho,
ccc                     equal to fractional rate of change with height
ccc                     of potential density at height, L_M, as
ccc                     normalised on -(fricv / L_M) ** 2 / gacc (a
ccc                     thermal layer is identified between ground and
ccc                     height L_M in which density gradient is treated
ccc                     as constant, and assigned a value as in reality
ccc                     holding at height L_M)
ccc     therc6  ..      in stable ambient stratification, potential-
ccc                     density gradient in upper thermal layer holding
ccc                     between heights L_M and iblht (where L_M<iblht),
ccc                     expressed as a fraction of gradient in lower
ccc                     layer
ccc     therc7  ..      value in ambient convection of k_rho, equal to
ccc                     fractional rate of change with height of
ccc                     potential density at height, -L_M, as
ccc                     normalised on (fricv / L_M) ** 2 / gacc (a
ccc                     thermal layer is identified between ground and
ccc                     height -L_M in which density gradient is treated
ccc                     as constant, and assigned a value as in reality
ccc                     holding at height -L_M)
ccc     therc8  ..      in ambient convection,
ccc                     -(w_c / fricv) / (iblht / L_M) ** (1 / 3), where
ccc                     w_c is velocity scale of free convection

        real*8 fricv,wdir,invMO,F0max,MTOM,wingar,F0,TOM,acc,ts,xs,zs,
     +         zsndMO,xresnd,concs,QBs,dQBTs,Tresnd,karho,karho2,karhou,
     +         karhu2,Tnd,delXnd,Zhnd,xnd,concf,agennd,dagxnd,dagTnd,cX,
     +         c1,c2,c3,dum,dum1,del,del1
ccc     F0      ..      kmF * F0max (N)
ccc     TOM     ..      kmM * MTOM (kg)
ccc     acc     ..      aircraft acceleration during ground run (m s-2)
ccc     ts      ..      [F0 / (rho * acc ** 4)] ** (1 / 6) (s), a
ccc                     nondimensionalising time scale in solution of
ccc                     dispersion before aircraft is rotated upward
ccc     xs      ..      [F0 / (rho * acc)] ** (1 / 3) (m), a
ccc                     nondimensionalising length scale in direction of
ccc                     release in solution of dispersion before
ccc                     aircraft is rotated upward
ccc     zs      ..      xs * ts / tau (m), a nondimensionalising length
ccc                     scale in vertical in solution of dispersion
ccc                     before aircraft is rotated upward
ccc     zsndMO  ..      zs as nondimensionalised on L_M
ccc     xresnd  ..      xres as nondimensionalised on xs
ccc     concs   ..      xs * ts / (F0 / rho) (s m-3), a
ccc                     nondimensionalising scale of concentration of
ccc                     conservative tracer in exhaust per unit rate
ccc                     tracer-mass release; used in solution of
ccc                     dispersion before aircraft is rotated upward
ccc     QBs     ..      F0 / (acc * tau) (kg s-1), a nondimensionalising
ccc                     scale of QB; used in converting solution for
ccc                     dispersion when upward rotation of aircraft
ccc                     begins to dimensional form appropriate for
ccc                     modelling dispersion of these exhausts
ccc                     subsequently
ccc     dQBTs   ..      QBs / ts (kg s-2), a nondimensionalising scale
ccc                     of rate of change of QB, and used as per QBs
ccc     Tresnd  ..      Tres nondimensionalised on ts, where Tres is
ccc                     temporal resolution of grid over which exhaust
ccc                     ages are solved for before aircraft is rotated
ccc                     upward, equal to 1 / ngsps (in s)
ccc     karho   ..      k_rho ** (1 / 2) * fricv * ts / abs(L_M)
ccc     karho2  ..      k_rho ** (1 / 2) * fricv * Tres2 / abs(L_M),
ccc                     where Tres2 is time step over which dispersion
ccc                     of exhausts released before aircraft is rotated
ccc                     upward is evaluated after rotation begins
ccc     karhou  ..      in stable ambient stratification,
ccc                     therc6 ** (1 / 2) * karho, ie. analogue of karho
ccc                     in upper of two thermal layers identified
ccc     karhu2  ..      in stable ambient stratification,
ccc                     therc6 ** (1 / 2) * karho2, ie analogue of
ccc                     karho2 in upper of two thermal layers identified
ccc     Tnd     ..      time since start of ground run, as
ccc                     nondimensionalised on ts
ccc     delXnd  ..      before aircraft is rotated upward, distance of
ccc                     exhaust-plume head in direction of release (in
ccc                     absence of wind), as measured from grid location
ccc                     furthest from source at which exhausts were
ccc                     present at last grid time processed, and as
ccc                     nondimensionalised on xs
ccc     Zhnd    ..      before aircraft is rotated upward, mean height
ccc                     exhausts rise buoyantly to in plume head, as
ccc                     nondimensionalised on zs
ccc     xnd     ..      before aircraft is rotated upward, distance
ccc                     from source in direction of release, as
ccc                     nondimensionalised on xs
ccc     concf   ..      before aircraft is rotated upward, factor used
ccc                     to calculate concentrations of conservative
ccc                     tracer per unit rate tracer-mass release
ccc     agennd  ..      before aircraft is rotated upward, solution for
ccc                     exhaust age (at given grid point), as
ccc                     nondimensionalised on ts
ccc     dagxnd  ..      before aircraft is rotated upward, partial
ccc                     derivative of exhaust age with respect to
ccc                     distance in direction of release, as
ccc                     nondimensionalised on ts / xs
ccc     dagTnd  ..      before aircraft is rotated upward, partial
ccc                     derivative of exhaust age with respect to time
ccc     cX      ..      before aircraft is rotated upward, distance in
ccc                     direction of release of grid location furthest
ccc                     from source at which exhausts were present at
ccc                     last grid time processed, measured from location
ccc                     of exhaust-plume head at this time (so cX<=0),
ccc                     and as nondimensionalised on xs

        parameter (nx=8000,ngsps=200,ng2sps=100,nitmax=50,iTrotm=60,
     +             nt=ngsps*iTrotm)
ccc     nx      ..      max number grid locations permitted in solution
ccc                     for dispersion before aircraft is rotated upward
ccc     ngsps   ..      number grid times over which solution for
ccc                     dispersion is obtained in each 1 s of ground run
ccc                     before aircraft is rotated upward
ccc     ng2sps  ..      number grid times over which solution for
ccc                     dispersion is obtained in each 1 s after upward
ccc                     rotation of aircraft begins
ccc     nitmax  ..      max number iterations of Newton's method
ccc                     permitted
ccc     iTrotm  ..      max duration ground run before aircraft is
ccc                     rotated upward (integral s)
ccc     nt      ..      ngsps * iTrotm, max number grid times stepped
ccc                     through before aircraft is rotated upward

        integer itrac(iTrotm)
ccc     itrac   ..      status of dispersion as a function of exhaust
ccc                     age, set to 0 when model ceases to apply

        real*8 agend(0:nx),znd(0:nx),area(iTrotm),z(iTrotm),
     +         xLol(iTrotm),conc(iTrotm),MBnd(nt),dMBTnd(nt),xL(iTrotm),
     +         QB(iTrotm),dQBT(iTrotm),QLF(iTrotm)
ccc     agend   ..      before aircraft is rotated upward, exhaust age
ccc                     as function distance from source in direction of
ccc                     release, and as nondimensionalised on ts
ccc     znd     ..      before aircraft is rotated upward, mean height
ccc                     of exhausts following from buoyant rise as
ccc                     function distance from source in direction of
ccc                     release, and as nondimensionalised on zs
ccc     area    ..      cross-sectional area of exhaust plume as
ccc                     function exhaust age (m2)
ccc     z       ..      mean height of exhausts following from buoyant
ccc                     rise as function exhaust age (m)
ccc     xLol    ..      value of xL at previous output time (m), as
ccc                     function previous exhaust age
ccc     conc    ..      concentration of conservative tracer in exhaust
ccc                     per unit rate tracer-mass release, as function
ccc                     exhaust age (s m-3)
ccc     MBnd    ..      rho * area * v / dagx, before aircraft is
ccc                     rotated upward and as function exhaust age
ccc                     (vals for oldest exhausts stored first), where v
ccc                     is (total) rate of change of z and dagx is
ccc                     partial differential of exhaust age with respect
ccc                     to distance in direction of release; ie. buoyant
ccc                     (vertical) momentum per unit difference in age
ccc                     of exhausts; and as nondimensionalised on
ccc                     B_F * ts
ccc     dMBTnd  ..      (total) rate of change of buoyant (vertical)
ccc                     momentum per unit difference in age of exhausts,
ccc                     before aircract is rotated upward and as
ccc                     function exhaust age (vals for oldest exhausts
ccc                     stored first); and as nondimensionalised on B_F
ccc     xL      ..      distance exhausts travel in direction of release
ccc                     from reference point collocated with them at
ccc                     release but static with respect to undisturbed
ccc                     air
ccc     QB      ..      rho * area * v (kg s-1), after upward rotation
ccc                     of aircraft begins and as function exhaust
ccc                     age; ie. buoyant (vertical) momentum per unit
ccc                     length in direction of release of exhausts
ccc                     released before rotation, as evaluated after
ccc                     rotation begins
ccc     dQBT    ..      rate of change of QB (kg s-2)
ccc     QLF     ..      rho * area * delta_u (kg s-1), after upward
ccc                     rotation of aircraft begins and as function
ccc                     exhaust age, where delta_u is speed of exhausts
ccc                     in direction of release with respect to
ccc                     undisturbed air; ie. excess release momentum
ccc                     per unit length in direction of release of
ccc                     exhausts released before rotation, as evaluated
ccc                     after rotation begins

        quarpi = datan(one)
        dtor = quarpi/45

        F0 = kmF*F0max
        TOM = kmM*MTOM
        acc = F0/TOM - rollco*gacc
           if (acc.lt.sig) then
            print *,'algorithm fails: prog ends'
            stop
           endif

ccc     calc airspeed of aircraft when rotated upward from runway
ccc     (m s-1)
        dum = corot*dsqrt(2*TOM*gacc/(CLmax*rho*wingar))

ccc        when ambient convection is absent
           if (invMO.ge.zero) then
            dum1 = (therc1 + therc2*invMO)*fricv

           else
            dum1 = one - therc3*invMO
            del = dum1**quart
            dum1 = (therc1 +
     +              therc4*(-quarpi + datan(del) + dlog(2/(one + del)) +
     +                      dlog(2/(one + dsqrt(dum1)))))*fricv

           endif

ccc     calc duration of ground run before rotation (s)
        iTrot = idint((dum - dum1*dcos(dtor*wdir))/acc)
           if (iTrot.lt.1) then
            print *,'algorithm fails: prog ends'
            stop
           endif

ccc     initialise last grid time since start of ground run at which
ccc     outputs were generated (integral s)
        iTlas = 0

ccc     initialise furthest grid location in direction of release from
ccc     source to have been processed previously
        mx = 0

ccc     initialise value of mx to be passed to next grid time to be
ccc     processed
        mxnew = 0

        dum = F0/rho
        xs = (dum/acc)**third
        ts = (dum/(acc**4))**sixth
        zs = xs*ts/tau
        zsndMO = zs*invMO
        concs = xs/(dum*ts)
        QBs = F0/(acc*tau)
        dQBTs = QBs/ts
        xresnd = xres/xs
        Tresnd = 1/(ngsps*ts)
        corat = (coke**2)/cokF
        dum = fricv*dabs(invMO)
        dum1 = dum/ng2sps
        dum = dum*ts

ccc        when ambient convection is absent
           if (invMO.ge.zero) then
            del = dsqrt(therc5)
            karho = del*dum
            karho2 = del*dum1
            del = dsqrt(therc6*therc5)
            karhou = del*dum
            karhu2 = del*dum1

           else
            del = dsqrt(therc7)
            karho = del*dum
            karho2 = del*dum1

           endif
        Tnd = zero
        delXnd = zero
        Zhnd = zero
        agend(0) = zero
        znd(0) = zero

        open (unit=50,file='TOdispRes.txt')

        write (50,*) 'nomenclature:'
        write (50,*) 't         .. exhaust age (s)'
        write (50,*)
     +  'Xs        .. distance source has moved through at time of'
        write (50,*) '             release (m)'
        write (50,*)
     +  'delta_xL  .. distance of exhausts in release direction from'
        write (50,*)
     +  '             reference point collocated with them at old time'
        write (50,*)
     +  '             but static with respect to undisturbed air (m)'
        write (50,*)
     +  'xL        .. distance of exhausts in release direction from'
        write (50,*)
     +  '             reference point collocated with them at releasce'
        write (50,*)
     +  '             but static with respect to undisturbed air (m)'
        write (50,*)
     +  'zb        .. lowest height exhausts extend to (m)'
        write (50,*)
     +  'zt        .. greatest height exhausts extend to (m)'
        write (50,*) 'delta_y   .. lateral spread of exhausts (m)'
        write (50,*)
     +  'conc      .. concentration of passive conservative tracer per'
        write (50,*)
     +  '             unit rate tracer-mass release (s m-3)'
        write (50,*) 'duration of ground run before rotation (s)'
        write (50,*) iTrot

100         iTfir = iTlas

            print *,'time since beg of grd run (integl secs, >',
     +      iTfir,' to continue) ?'
            read *,iTlas

               if (iTlas.gt.iTfir) then

                iTold = iTfir

ccc                    when some of ground run prior to rotation
ccc                    remains to be processed
120                    if (iTold.lt.iTrot) then

ccc                     set start grid time to be processed
                        mTolgs = iTfir*ngsps + 1

ccc                     set end grid time to be processed
                        mTgs = min(iTrot,iTlas)*ngsps

ccc                         step through grid times
                            do 500 iTgs=mTolgs,mTgs

ccc                         update cX
                            cX = (mxnew - mx)*xresnd - delXnd
                               if (cX.lt.-xresnd-sig .or. cX.gt.sig)
     +                         then
                                print *,'algorithm fails: prog ends'
                                stop
                               endif

ccc                         update furthest grid location in direction
ccc                         of release from source to have been
ccc                         processed
                            mx = mxnew

                            m = mx + 1

ccc                         initialise youngest appropriate gridded age
                            iolagt = 1

                            Tnd = Tnd + Tresnd

ccc                             update MBnd and dMBTnd in case of
ccc                             exhausts released before last grid time
ccc                             processed, according to exhaust
ccc                             characteristics at each grid location,
ccc                             ix, at last grid time
                                do 140 ix=1,m

ccc                                calc oldest appropriate gridded age
                                   if (ix.le.mx) then
                                    n = idint(agend(ix)/Tresnd)
                                   else
                                    n = iTgs - 1
                                   endif
                                   if (n.ge.iTgs) then
                                    print *,'algorithm fails: prog ends'
                                    stop
                                   endif

                                    do 130 i=iolagt,n
                                    j = iTgs - i
                                       if (j.lt.1 .or. j.ge.iTgs) then
                                        print *,
     +                                  'algorithm fails: prog ends'
                                        stop
                                       endif
                                       if (ix.le.mx) then
                                        dum = znd(ix)
                                       else
                                        dum = Zhnd
                                       endif

ccc                                 obtain mean height exhausts have
ccc                                 risen buoyantly to, as
ccc                                 nondimensionalised on Monin-Obukhov
ccc                                 length
                                    dum = dum*zsndMO

ccc                                    when centroid in lower ambient
ccc                                    thermal layer
                                       if (dabs(dum).lt.one) then
                                        del = karho

ccc                                    when centroid in upper ambient
ccc                                    thermal layer (in stable
ccc                                    stratification)
                                       else if (dum.ge.one) then
                                        del = karhou

                                       else
                                        del = zero
                                       endif
                                    dum1 = del*Tresnd

ccc                                    when stable ambient
ccc                                    stratification is absent
                                       if (invMO.le.zero) then
                                           if (dum1.lt.sig) then
                                            dum = Tresnd
                                           else
                                            dum = dsinh(dum1)/del
                                           endif
                                        dum1 = dcosh(dum1)
                                        k = 1

                                       else
                                           if (dum1.lt.sig) then
                                            dum = Tresnd
                                           else
                                            dum = dsin(dum1)/del
                                           endif
                                        dum1 = dcos(dum1)
                                        k = -1
                                       endif

                                    del1 = dum*dMBTnd(j) + dum1*MBnd(j)
                                    dMBTnd(j) = dum1*dMBTnd(j) +
     +                                          k*(del**2)*dum*MBnd(j)
                                    MBnd(j) = del1
130                                 continue

ccc                             update iolagt as necessary
                                if (ix.le.mx) iolagt = n + 1

140                             continue

                            dum = karho*Tresnd

ccc                             update MBnd and dMBTnd in case of
ccc                             exhausts released at last grid time
ccc                             processed
                                if (dum.lt.sig) then
                                 MBnd(iTgs) = Tresnd
                                else if (invMO.le.zero) then
                                 MBnd(iTgs) = dsinh(dum)/karho
                                else
                                 MBnd(iTgs) = dsin(dum)/karho
                                endif
                                if (invMO.le.zero) then
                                 dMBTnd(iTgs) = dcosh(dum)
                                else
                                 dMBTnd(iTgs) = dcos(dum)
                                endif

ccc                         initialise value of mx to be inherited at
ccc                         next grid time
                            mxnew = mx + 1

ccc                         initialise grid location to be processed
                            ix = 0

                            xnd = zero
ccc                             step through grid locations, from
ccc                             nearest to furthest from source in
ccc                             direction of release
200                             ix = ix + 1

                                xnd = xnd + xresnd

ccc                                when first processing a grid location
ccc                                at which exhausts were not present
ccc                                at last grid time processed
                                   if (ix.eq.mx+1) then
                                    i = 0

ccc                                 set iteration invariants
                                    c1 = mx*xresnd - (Tnd**2)/2
                                    c2 = Tnd - agend(mx)
                                    c3 = cX/Tresnd - Tnd

ccc                                 set initial estimate of delXnd for
ccc                                 improvement by Newton's method
                                    delXnd = Tnd*Tresnd - cX

ccc                                     step through Newton's method
ccc                                     iterations
250                                     i = i + 1
                                        dum = c1 + delXnd
                                        dum1 = c2/delXnd
                                        del1 = delXnd/Tresnd + c3

ccc                                     calc function whose root with
ccc                                     respect to dimensionless
ccc                                     location of exhaust-plume head
ccc                                     is to be obtained
                                        del = corat*(dum**2)*del1 - dum1

ccc                                     calc derivative of this function
ccc                                     with respect to dimensionless
ccc                                     location of exhaust-plume head
                                        del1 =
     +                                  corat*dum*(dum/Tresnd +
     +                                             2*del1) + dum1/delXnd
                                           if (del1.lt.sig) then
                                            print *,
     +                                      'algorithm fails: prog ends'
                                            stop
                                           endif

ccc                                        when function is
ccc                                        significantly valued
                                           if (dabs(del).ge.sig) then

ccc                                         calc improvement on root
                                            delXnd = delXnd - del/del1
                                               if (delXnd.lt.
     +                                             Tnd*Tresnd-cX) then
                                                print *,
     +                                          'alg fails: prog ends'
                                                stop
                                               endif

                                           endif

                                        if (i.lt.nitmax .and.
     +                                      dabs(del).ge.sig) goto 250
                                       if (i.gt.nitmax .or.
     +                                     (i.eq.nitmax .and.
     +                                      dabs(del).ge.sig)) then
                                        print *,
     +                                  'algorithm fails: prog ends'
                                        stop
                                       endif

ccc                                calc value of mx to be inherited at
ccc                                next grid time
                                    mxnew = mx + idint(delXnd/xresnd)
                                       if (mxnew.gt.nx) then
                                        print *,
     +                                  'algorithm fails: prog ends'
                                        stop
                                       endif

                                    dum = (delXnd + cX -
     +                                     Tnd*Tresnd)/cokF

ccc                                 update mean height of exhaust-plume
ccc                                 head resulting from buoyant rise
                                    Zhnd = Zhnd + dum*MBnd(1)

                                    concf = dum/Tresnd
                                    if (iTgs.eq.mTgs .and.
     +                                  iTlas.ge.iTrot) dagxnd = dum1
                                   endif

ccc                                when exhausts are present at least
ccc                                one grid location further from source
ccc                                than at last grid time processed,
ccc                                solve at associated grid points by
ccc                                linear extrapolation
                                   if (ix.gt.mx .and. mxnew.gt.mx) then
                                    dum = (ix - mx)*xresnd/delXnd
                                    agend(ix) = agend(mx) + dum*c2
                                    znd(ix) = znd(mx) +
     +                                        dum*(Zhnd - znd(mx))

ccc                                when processing grid locations at
ccc                                which exhausts were present at last
ccc                                grid time processed
                                   else if (ix.le.mx) then

                                    i = 0

ccc                                 set initial estimate of
ccc                                 dimensionless exhaust age for
ccc                                 improvement by Newton's method
                                    agennd = agend(ix-1)

ccc                                     step through Newton's method
ccc                                     iterations
300                                     i = i + 1
                                        dum = (agennd -
     +                                         agend(ix-1))/xresnd
                                        dum1 = one -
     +                                         (agennd - agend(ix))/
     +                                          Tresnd - dum*Tnd
                                        del = xnd -
     +                                        (Tnd - agennd/2)*agennd

ccc                                     calc derivative of function
ccc                                     whose root with respect to
ccc                                     dimensionless exhaust age is
ccc                                     sought
                                        del1 = -corat*del*
     +                                          (del*(Tnd/xresnd +
     +                                                1/Tresnd) +
     +                                           2*dum1*
     +                                            (Tnd - agennd)) -
     +                                         2*dum/xresnd

ccc                                     calc function
                                        dum = corat*dum1*del**2 - dum**2

                                           if (dum1.lt.sig .or.
     +                                         del.lt.sig .or.
     +                                         del1.gt.-sig) then
                                            print *,
     +                                      'algorithm fails: prog ends'
                                            stop
                                           endif

ccc                                     when function is significantly
ccc                                     valued, calc improvement on root
                                        if (dabs(dum).ge.sig)
     +                                  agennd = agennd - dum/del1

                                        if (i.lt.nitmax .and.
     +                                      dabs(dum).ge.sig) goto 300

                                       if (i.gt.nitmax .or.
     +                                     (i.eq.nitmax .and.
     +                                      dabs(dum).ge.sig) .or.
     +                                     agennd.lt.
     +                                     agend(ix-1)-sig1 .or.
     +                                     agennd.gt.agend(ix)+sig1)
     +                                 then
                                        print *,
     +                                  'algorithm fails: prog ends'
                                        print *,'try increasing resn'
                                        stop
                                       endif

                                    dagTnd = (agennd - agend(ix))/Tresnd
                                    agend(ix) = agennd
                                    dagxnd = (agend(ix) -
     +                                        agend(ix-1))/xresnd

ccc                                 establish appropriate value of MBnd
ccc                                 by linear extrapolation, for
ccc                                 calculation of mean height of
ccc                                 exhausts resulting from buoyant rise
                                    del = agend(ix)/Tresnd
                                    i = idint(del)
                                    j = iTgs - i + 1
                                       if (j.lt.1 .or. j.gt.iTgs+1) then
                                        print *,'alg fails: prog ends'
                                        stop
                                       endif
                                       if (j.gt.iTgs) then
                                        dum = zero
                                       else
                                        dum = MBnd(j)
                                       endif
                                       if (j.eq.1) then
                                        dum1 = MBnd(1)
                                       else
                                        dum1 = MBnd(j-1)
                                       endif
                                    del1 = dum +
     +                                     (del - dble(i))*(dum1 - dum)

                                    dum = one - dagTnd
                                    dum1 = dum/xresnd
                                    dum = (dum - dagxnd*Tnd)/cokF
                                    del = dagxnd/Tresnd
                                    c1 = znd(ix)*del + znd(ix-1)*dum1
                                    del = del + dum1
                                       if (del.lt.sig) then
                                        print *,
     +                                  'algorithm fails: prog ends'
                                        stop
                                       endif
                                    dum1 = (dum*del1 + c1)/del
                                       if ((invMO.le.zero .and.
     +                                      (dum1.lt.znd(ix-1)-sig1 .or.
     +                                       dum1.gt.znd(ix)+sig1)) .or.
     +                                     (invMO.gt.zero .and.
     +                                      dum1.lt.-sig1)) then
                                        print *,
     +                                  'algorithm fails: prog ends'
                                        print *,'try increasing resn'
                                        stop
                                       endif
                                    znd(ix) = dum1
                                    concf = dum

                                   endif

ccc                                when grid time furthest into ground
ccc                                run currently for processing is
ccc                                reached, calc dimensional exhaust
ccc                                characteristics as function exhaust
ccc                                age, with ages as separated by 1 s
ccc                                intervals: this is for purposes of
ccc                                output, or for initialisation of
ccc                                dispersion regime after upward
ccc                                rotation of aircraft is begun
                                   if (iTgs.eq.mTgs) then

ccc                                 initialise youngest appropriate age,
ccc                                 in integral s
                                    if (ix.eq.1) magnex = 1

ccc                                    when considering a grid location
                                       if (ix.le.mxnew) then

ccc                                     calc oldest appropriate age, in
ccc                                     integral s
                                        iage = idint(agend(ix)*ts)
                                           if (iage.lt.0 .or.
     +                                         iage.gt.iTgs/ngsps .or.
     +                                         iage.gt.iTrot) then
                                            print *,
     +                                      'alg fails: prog ends'
                                            stop
                                           endif

                                            do 350 i=magnex,iage
                                            j = i - (nage - iTfir)
                                            if (j.gt.0) xLol(j) = xL(j)
                                            dum = (xnd -
     +                                             (Tnd - agend(ix)/2)*
     +                                              agend(ix))*xs
                                            xL(i) = dum
                                            area(i) = (coke*dum)**2
                                            z(i) = znd(ix)*zs
                                            dum = concf*concs
                                            if (ix.le.mx)
     +                                      dum = dum/dagxnd
                                            conc(i) = dum

ccc                                            when upward rotation of
ccc                                            aircraft has been begun,
ccc                                            calc QB, dQBT and QBF so
ccc                                            dispersion may be
ccc                                            evaluated later
                                               if (iTlas.ge.iTrot) then
                                                j = mTgs - (i - 1)*ngsps
                                                   if (j.lt.1 .or.
     +                                                 j.gt.mTgs) then
                                                    print *,'alg fails:'
                                                    print *,'prog ends'
                                                    stop
                                                   endif
                                                QB(i)=
     +                                          MBnd(j)*dagxnd*QBs
                                                dQBT(i)=
     +                                          dMBTnd(j)*dagxnd*dQBTs
                                                QLF(i) =
     +                                          cokF*F0*area(i)*conc(i)
                                               endif

350                                        continue
                                        magnex = iage + 1
                                       endif

ccc                                    when there are no more grid
ccc                                    locations to be processed,
ccc                                    draw on characteristics of
ccc                                    exhaust-plume head
                                       if (ix.ge.mxnew) then

ccc                                     calc oldest appropriate age, in
ccc                                     integral s
                                        iage = mTgs/ngsps

                                            do 360 i=magnex,iage
                                            j = i - (nage - iTfir)
                                            if (j.gt.0) xLol(j) = xL(j)
                                            dum = (mx*xresnd + delXnd -
     +                                             (Tnd**2)/2)*xs
                                            xL(i) = dum
                                            area(i) = (coke*dum)**2
                                            z(i) = Zhnd*zs
                                            conc(i) = concf*concs

ccc                                            when upward rotation of
ccc                                            aircraft begins,
ccc                                            calc QB, dQBT and QBF so
ccc                                            dispersion may be
ccc                                            evaluated later
                                               if (iTlas.ge.iTrot) then
                                                QB(i) =
     +                                          MBnd(ngsps)*dagxnd*QBs
                                                dQBT(i) =
     +                                          dMBTnd(ngsps)*dagxnd*
     +                                           dQBTs
                                                QLF(i) =
     +                                          cokF*F0*area(i)*conc(i)
                                               endif
360                                         continue
                                       endif

                                   endif

ccc                             when grid locations remain to be
ccc                             processed, progress to next
                                if (ix.lt.mxnew) goto 200

500                         continue

ccc                        when upward rotation of aircraft begins,
ccc                        prepare to evaluate dispersion of exhausts
ccc                        released earlier
                           if (iTlas.gt.iTrot) then
                            iTold = iTrot
                            goto 120
                           endif

ccc                    when upward rotation of aircraft has been begun
                       else

ccc                     calc number of grid times to be processed
                        mTgs = (iTlas - iTold)*ng2sps

ccc                        when rotation had been begun by last grid
ccc                        time at which outputs were generated, store
ccc                        val of xL at grid time of last output
                           if (iTfir.ge.iTrot) then
                                do 600 iage=1,iTrot
                                xLol(iage) = xL(iage)
600                             continue
                           endif

ccc                         step through grid times
                            do 1000 iTgs=1,mTgs

ccc                             step through gridded ages, as separated
ccc                             by 1 s intervals
                                do 900 iage=1,iTrot

                                   if (area(iage).lt.sig .or.
     +                                 QLF(iage).lt.sig) then
                                    print *,'algorithm fails: prog ends'
                                    stop
                                   endif
                                dum = z(iage)*invMO

ccc                                when mean height exhausts have risen
ccc                                to buoyantly lies in lower ambient
ccc                                thermal layer
                                   if (dabs(dum).lt.one) then
                                    del = karho2

ccc                                when mean height exhausts have risen
ccc                                to buoyantly lies in upper ambient
ccc                                thermal layer (in stable
ccc                                stratification)
                                   else if (dum.ge.one) then
                                    del = karhu2

                                   else
                                    del = zero
                                   endif

ccc                                when stable ambient stratification
ccc                                is absent
                                   if (invMO.le.zero) then
                                       if (del.lt.sig) then
                                        dum = one
                                       else
                                        dum = dsinh(del)/del
                                       endif
                                    dum1 = dcosh(del)
                                    k = 1

                                   else
                                       if (del.lt.sig) then
                                        dum = one
                                       else
                                        dum = dsin(del)/del
                                       endif
                                    dum1 = dcos(del)
                                    k =-1
                                   endif

ccc                             calc dimensional exhaust characteristics
                                del1 = dum*dQBT(iage)/ng2sps +
     +                                 dum1*QB(iage)
                                dQBT(iage) =
     +                          dum1*dQBT(iage) +
     +                          k*ng2sps*(del**2)*dum*QB(iage)
                                c1 = QB(iage)
                                QB(iage) = del1
                                del = 2*ng2sps*rho
                                del1 = QLF(iage)/del
                                dum = 3*coke*del1/area(iage)**cthrhf
                                dum1 =
     +                          dsqrt(one + (c1/QLF(iage))**2) +
     +                          dsqrt(one + (QB(iage)/QLF(iage))**2)
                                dum = (one + dum*dum1)**twothd
                                c2 = area(iage)
                                area(iage) = dum*area(iage)
                                conc(iage) = conc(iage)/dum
                                del = 2*del
                                xL(iage) = xL(iage) +
     +                                     del1*(1/c2 + 1/area(iage))
                                z(iage) =
     +                          z(iage) + (c1/c2 +
     +                                     QB(iage)/area(iage))/del

900                             continue
1000                        continue

                       endif

                write (50,*) 'old time since start of ground run (s)'
                write (50,*) iTfir
                write (50,*) 'new time since start of ground run (s)'
                write (50,*) iTlas
                write (50,*) 'number ages tracked'

ccc             initialise number of gridded exhaust ages at which
ccc             dispersive characteristics are available
                ntr = 0
                ntrold = 0

ccc             set up output-loop invariants
                c1 = coke/rho

ccc                set nage to number of gridded ages
                   if (iTlas.gt.iTrot) then
                    nage = iTrot
                   else
                    nage = iTlas
                    c1 = c1*cokF*F0
                   endif

ccc                when old time preceded upward rotation of aircraft,
ccc                update itrac and include exhausts emitted after old
ccc                time
                   if (iTfir.lt.iTrot) then
                    j = nage - iTfir
                        do 1050 i=nage,1,-1
                           if (i.gt.j) then
                            itrac(i) = itrac(i-j)
                           else
                            itrac(i) = 1
                           endif
1050                    continue
                   endif

ccc             set c2 to velocity scale of free convection (in m s-1)
                dum = -iblht*invMO
                if (dum.gt.zero) c2 = therc8*(dum**third)*fricv

ccc                 step through set of exhaust ages three times, first
ccc                 to obtain nunber of gridded ages at which dispersive
ccc                 characteristics are available, second to output
ccc                 these characteristics, third to output reason where
ccc                 characteristics are not available
                    do 1500 igo=1,3

                       if (igo.eq.2) then
                        write (50,*) ntr
                        if (ntr.gt.0) write (50,'(a57)')
     +       '   t,  Xs,delta_xL,    xL,    zb,    zt,delta_y,     conc'
                       else if (igo.eq.3 .and. ntr.lt.ntrold) then
                        write (50,*)
     +                  'ages ceased tracking since old time'
                        write (50,'(a36)')
     +                  '   t,  Xs,                    reason'
                       endif

ccc                     step through gridded ages
                        do 1100 iage=1,nage

ccc                        when dispersion model was applicable at old
ccc                        time
                           if (itrac(iage).eq.1) then

                            if (igo.eq.1) ntrold = ntrold + 1

                               if (area(iage).lt.sig) then
                                print *,'algorithm fails: prog ends'
                                stop
                               endif

                            del = dsqrt(area(iage))/2

ccc                            when plume is in contact with ground
                               if (z(iage).lt.del) then
                                del1 = 2*z(iage)

                               else
                                del1 = z(iage) + del
                               endif

                               if (del1.lt.sig) then
                                print *,'algorithm fails: prog ends'
                                stop
                               endif

ccc                            when exhausts have risen buoyantly to a
ccc                            mean height greater than magnitude of
ccc                            Monin-Obukhov length in ambient
ccc                            convection, set velocity scale of ambient
ccc                            diffusion to that of free convection
                               if (invMO*z(iage).lt.-one) then
                                dum = c2

ccc                            set velocity scale of ambient diffusion
ccc                            to friction velocity
                               else
                                dum = fricv

                               endif

                               if (iTlas.le.iTrot) then
                                dum1 = c1*conc(iage)
                               else
                                dum1 = c1*dsqrt(QB(iage)**2 +
     +                                          QLF(iage)**2)/area(iage)
                               endif

ccc                            when exhaust dispersion reflects trapping
ccc                            at inversion capping boundary layer or
ccc                            action of ambient turbulence
                               if (del1.gt.dble(iblht) .or.
     +                             dum1.lt.dum) then
                                   if (igo.eq.3) then
                                    itrac(iage) = 0
                                    i = iage + iTlas - nage
                                    j = idnint(acc*((iTlas - i)**2)/2)
                                       if (del1.gt.dble(iblht) .and.
     +                                     dum1.ge.dum) then
                                        write (50,'(i4,i5,a27)')i,j,
     +                                  'base of inversion reached'
                                       else if (del1.le.
     +                                          dble(iblht) .and.
     +                                          dum1.lt.dum) then
                                        write (50,'(i4,i5,a27)') i,j,
     +                                  'ambient diffusion underway'
                                       else
                                        write (50,'(i4,i5,a55)') i,j,
     +          'ambient diffusion underway / base of inversion reached'
                                       endif
                                   endif

                               else

                                   if (igo.eq.1) then
                                    ntr = ntr + 1
                                   else if (igo.eq.2) then
                                    i = iage + iTlas - nage
                                    j = idnint(acc*((iTlas - i)**2)/2)
                                       if (iTfir.ge.iTrot) then
                                        k = 0
                                       else
                                        k = nage - iTfir
                                       endif
                                       if (iage.gt.k) then
                                        dum = xL(iage) - xLol(iage-k)
                                       else
                                        dum = xL(iage)
                                       endif

ccc                                    when plume is in contact with
ccc                                    ground
                                       if (z(iage).lt.del) then
                                        dum1 = zero
                                        del = area(iage)/del1

                                       else
                                        dum1 = z(iage) - del
                                        del = 2*del
                                       endif

                                    write
     +                              (50,'(i4,i5,f9.1,3f7.1,f8.1,e10.3)')
     +                              i,j,dum,xL(iage),dum1,del1,del,
     +                              conc(iage)
                                   endif
                               endif

                           endif

1100                    continue

1500                continue

                   if (ntrold.lt.1 .or. ntrold.lt.ntr) then
                    print *,'algorithm fails: prog ends'
                    stop
                   endif

                if (ntr.gt.0) goto 100

               endif

        close (50)

        end