anchor.diff:      enddo
            ! compute sigma coordinate quantities (bottom to top).
       do k=1,nlevs+1
  -        si(nlevs+2-k)= ak(k)/101300.0+bk(k) ! si are now sigmas
  +        si(nlevs+2-k)= (ak(k) - model_top_pressure)/(py_pdryini - model_top_pressure)+bk(k) ! si are now sigmas
       enddo
       do k=1,nlevs
           sl(k) = 0.5*(si(k)+si(k+1))
  
positive.code:     !call spectogrd(lnpsspec, psg)
      !psg = 1000.*exp(psg) ! convert to Pa
      !call grdtospec(log(psg), lnpsspec) ! back to spectral.
      !JOY adding read from namelist to eliminate using sigio to read ak,bk
      ak = (/ 200.00000000,566.898010,1290.53296,2210.97900&
      ,3376.51611,4844.03613,6678.60791,8913.76660,&
      11343.6543,13671.4268,15613.5645,16912.1309,&
      17364.6582,16959.9941,15812.9258,14089.5352,&
      11991.4277,9731.80664,7508.53223,5481.14404,&
      3757.14209,2389.20508,1381.52600,701.453003,&
      292.576996,86.4570007,11.6350002,8.99999961E-03,0.00000000 /)
  
      bk = (/ 0.00000000,0.00000000,0.00000000,0.00000000,&
              0.00000000,0.00000000,0.00000000,3.78680008E-04,&
              3.93276010E-03,1.43262697E-02,3.49500999E-02,&
              6.86749965E-02,0.117417611,0.180667043,0.256084293,&
              0.340293199,0.429195374,0.518456697,0.604054928,&
              0.682746828,0.752347112,0.811784863,0.860974848,&
              0.900580883,0.931749582,0.955872416,0.974401832,&
              0.988725841,1.00000000 /)
  
  
      do k=1,nlevs
  
negative.code:     !call init_griddata()
      integer, intent(in):: py_damping_level
      real(c_double), intent(in):: py_taustratdamp
      real(c_double), intent(in):: model_top_pressure
  
      ! initialize spherical harmonic lib
      !print *,'Initialising shtns'
      call shtns_init(nlons,nlats,ntrunc,nthreads=1,polar_opt=polar_opt)
  