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:     ! ******** Column average: (25) of SI-HN17 ********
  
      astarbaro(:,:) = astarbaro(:,:)+(astar1(:,:,k)    &
      + astar2(:,:,k))*exp(-zk/hh)*dc
      ua1baro(:,:) = ua1baro(:,:)+ua1(:,:)*exp(-zk/hh)*dc
      ua2baro(:,:) = ua2baro(:,:)+ua2(:,:)*exp(-zk/hh)*dc
      ep1baro(:,:) = ep1baro(:,:)+ep1(:,:)*exp(-zk/hh)*dc
      ep2baro(:,:) = ep2baro(:,:)+ep2(:,:)*exp(-zk/hh)*dc
      ep3baro(:,:) = ep3baro(:,:)+ep3(:,:)*exp(-zk/hh)*dc
      do j = 6,nd  ! ### yet to be multiplied by cosine
        ubaro(:,j) = ubaro(:,j)+uu(:,j+90,k)*exp(-zk/hh)*dc
        urefbaro(j) = urefbaro(j)+uref(j-5,k)*exp(-zk/hh)*dc
      enddo
    enddo
  
  
negative.code:   REAL :: z(kmax)
    REAL :: u(nd,kmax)
  
    a = 6378000.
    pi = acos(-1.)
    om = 7.29e-5
    dp = pi/180.
    dz = 500.
    h = 7000.
    r = 287.
    rkappa = r/1004.
  
    ! *** Default values for boundary ***
    !jb = 5
    !jd = 86 ! nd - lower bounding latitude
  
  