anchor.diff:          ep1(i,j) = ep1(i,j)+0.5*vv(i,j+90,k)**2    !F3a+b
           ep11 = 0.5*(pt(i,j+90,k)-tref(j-5,k))**2         !F3c
           zz = dz*float(k-1)
  -        ep11 = ep11*(r/hh)*exp(-rkappa*zz/hh)
  +        ep11 = ep11*(r/h)*exp(-rkappa*zz/h)
           ep11 = ep11*2.*dz/(tg(k+1)-tg(k-1))
           ep1(i,j) = ep1(i,j)-ep11                   !F3
           phip = dp*float(j)
  
positive.code:   REAL :: z(kmax)
    REAL :: u(nd,kmax)
  
    pi = acos(-1.)
    dp = pi/float(jmax-1)
    rkappa = rr/cp
  
  
    ! *** Default values for boundary ***
    !prefac = 6745.348
    !jb = 5
    !jd = 86 ! nd - lower bounding latitude
  
  
negative.code: SUBROUTINE compute_flux_dirinv(pv,uu,vv,pt,tn0,ts0,statn,stats,qref,uref,tref,fawa,ubar,tbar,&
          imax, JMAX, kmax, nd, nnd, jb, jd,&
          astarbaro,ubaro,urefbaro,ua1baro,ua2baro,ep1baro,ep2baro,ep3baro,ep4,astar1,astar2)
  
    INTEGER, INTENT(IN) :: imax, JMAX, kmax, nd, nnd, jb, jd
    REAL, INTENT(IN) :: pv(imax,jmax,kmax),uu(imax,jmax,kmax),vv(imax,jmax,kmax),pt(imax,jmax,kmax),&
            tn0(kmax),ts0(kmax),statn(kmax),stats(kmax),qref(nd,kmax),uref(jd,kmax),tref(jd,kmax),&
            fawa(nd,kmax),ubar(nd,kmax),tbar(nd,kmax)
    REAL, INTENT(OUT) :: astarbaro(imax,nd),ubaro(imax,nd),urefbaro(nd),ua1baro(imax,nd),ua2baro(imax,nd),&
            ep1baro(imax,nd),ep2baro(imax,nd),ep3baro(imax,nd),ep4(imax,nd),astar1(imax,nd,kmax),astar2(imax,nd,kmax)
  
  