[Mitgcm-support] Re: KPP (c30) bug

[mitgcm-support at dev.mitgcm.org]

Ralf, thanks for picking up the bug.

I have verified that value of KRi_range(Nrp1) makes absolutely
no difference to forward results (see test code below), but I
presume that an uninitialized value will cause havock with TAMC.

I suggest fixing this problem by intitializing KRi_range to zero
once, right after local variable declaration:

      p0 = 0.0
      p25 = 0.25
      p5 = 0.5
      p2 = 2.0

      do k = 0, Nrp1
         KRi_range(k) = p0
      end do

Please send complete adjoint modified code when ready.  Detlef, Armin and I
will verify your changes and try to find a solution for shear instability
term.

Dimitris


P.S. Here is test code for the z121 filter.  Incidentally, we never actually
use this filter in any of the JPL integrations because it does not seem to
make any difference to shear instability-related problems.  Instead we apply
horizontal 121 filters.


      PROGRAM test121
      integer    Nr, Nrp1, mr
      parameter (Nr = 22, Nrp1 = 23)
      real diffus (0:Nrp1)
C     DATA diffus / 0,-0.70,-0.47,-0.44,-0.14,-0.02,0.06,0.01,0.12,
C    &     -0.21,0.15,0.42,0.41,0.42,0.86,0.58,0.58,0.60,0.66,0.82,
C    &     0.83,0.92,1.05,0.0 /
      DATA diffus / 0,-0.70,-0.47,-0.44,-0.14,-0.02,0.06,0.01,0.12,
     &     -0.21,0.15,0.42,0.41,0.42,0.86,0.58,0.58,0.60,0.66,0.55,
     &     0.33,0.44,.66,999999999 /

      print*,diffus

      do mr = 1, 1
         call z121 (
     U     diffus)
      end do

      print*,' '
      print*,diffus
      
      stop
      end

c*************************************************************************

      subroutine z121 (
     U     v )

c     Apply 121 smoothing in k to 2-d array V(i,k=1,Nr)
c     top (0) value is used as a dummy
c     bottom (Nrp1) value is set to input value from above.

c     Note that its important to exclude from the smoothing any points
c     that are outside the range of the K(Ri) scheme, ie.  > 0.8, or <0.0.
c     Otherwise, there is interfence with other physics, especially
c     penetrative convection.

      IMPLICIT NONE

      integer    Nr, Nrp1, mr
      parameter (Nr = 22, Nrp1 = 23)
      integer    imt
      parameter (imt=1)

c input/output 
c-------------
c v     : 2-D array to be smoothed in Nrp1 direction
      real v(imt,0:Nrp1)

c local
      real p25, p5, p2, zwork, zflag,Riinfty
      parameter (Riinfty=.7)
      real KRi_range(0:Nrp1)
      integer i, k, km1, kp1

      p25 = 0.25
      p5 = 0.5
      p2 = 2.0

      do i = 1, imt

         v(i,Nrp1) = v(i,Nr)
         KRi_range(Nrp1) = -1.e20

         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) +
     &           KRi_range(k) * zwork
            if (k.eq.Nr) then
c           print*,' '
c           print*,km1,k,kp1,zflag,v(i,k),v(i,kp1)
c           print*,KRi_range(k),KRi_range(kp1),zwork
c           endif
            zwork = KRi_range(k) * zflag
            zflag = p2 + KRi_range(k)*(KRi_range(kp1)+KRi_range(km1))
            v(i,k) = v(i,k) / zflag
            if (k.eq.Nr) then
c           print*,zwork,zflag,v(i,k)
c           print*,' '
            endif
         end do

      end do

      return
      end