[Mitgcm-support] KPP (c30) bug

[mitgcm-support at dev.mitgcm.org]

There is a bug in c30 version of kpp_routines!

KRi_range(Nr+1) is not defined but referenced.
Have a look at:
      do i = 1, imt

         v(i,Nrp1) = v(i,Nr)

         do k = 1, Nr
            KRi_range(k) = p5 + SIGN(p5,v(i,k))
            KRi_range(k) = KRi_range(k) *
     &                     ( p5 + SIGN(p5,(Riinfty-v(i,k))) )
         end do

         zwork  = KRi_range(1) * v(i,1)
         v(i,1) = p2 * v(i,1) +
     &            KRi_range(1) * KRi_range(2) * v(i,2)
         zflag  = p2 + KRi_range(1) * KRi_range(2)
         v(i,1) = v(i,1) / zflag

         do k = 2, Nr
            km1 = k - 1
            kp1 = k + 1
            zflag = v(i,k)
            v(i,k) = p2 * v(i,k) +
     &           KRi_range(k) * KRi_range(kp1) * v(i,kp1) +
crg                                       ^^^
     &           KRi_range(k) * zwork
            zwork = KRi_range(k) * zflag
            zflag = p2 + KRi_range(k)*(KRi_range(kp1)+KRi_range(km1))
crg                                              ^^^
            v(i,k) = v(i,k) / zflag
         end do

      end do

It is fixed by using kp1 = min( k + 1, Nr ) in KRi_range(kp1).
I also modified the routine for TAMC.
It needs to be verified.

      _RS     zwork, zflag, ztmp1, ztmp
      integer KRi_range(0:Nr)
      integer i, k, km1, kp1

#ifdef ALLOW_AUTODIFF_TAMC
C--   dummy assignment to end declaration part for TAMC
      i = 0

C--   HPF directive to help TAMC
CHPF$ INDEPENDENT
#endif /* ALLOW_AUTODIFF_TAMC */
      do i = 1, imt

         KRi_range(0) = 0
         do k = 1, Nr
            if (v(i,k) .lt. 0.) then
               KRi_range(k) = 0
            else
               if (Riinfty-v(i,k) .lt. 0.) then
                  KRi_range(k) = 0
               else
                  KRi_range(k) = 1
               endif
            endif
         end do

         zwork  = 0.0

CADJ INIT z121tape = common, Nr
         do k = 1, Nr
CADJ STORE v(i,k), zwork = z121tape
            km1 = k - 1
            kp1 = min( k + 1, Nr )
            zflag = v(i,k)
            ztmp  = 2 * v(i,k)
     &            + KRi_range(k) * KRi_range(kp1) * v(i,kp1)
     &            + KRi_range(k) * zwork
            zwork = KRi_range(k) * zflag
            v(i,k) = ztmp
     &            / (2 + KRi_range(k)*(KRi_range(kp1)+KRi_range(km1)))
         end do

      end do


-- 
###################################################
 Ralf Giering
 FastOpt, Hamburg, Germany
 Tel. +49 4154 7938935
 Email: Ralf.Giering at FastOpt.de
 URL: http://www.FastOpt.de
###################################################