anchor.diff:    ep2baro(:,:) = 0.
     ep3baro(:,:) = 0.
     ep4(:,:) = 0.
  -  dc = dz/prefac
  +  dc = dz/6745.348
   
  -  do k = 2,kmax-1
  +  do k = 2,96
       zk = dz*float(k-1)
       do i = 1,imax
         do j = 6,nd-1            ! 5N and higher latitude
  
positive.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, aa, omega, dz, hh, rr, cp, prefac, &
          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(in) :: aa, omega, dz, hh, rr, cp, prefac
    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)
  
  
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
  
  