Fix out-of-bounds access on zk in cvmix_kpp_compute_StokesXi

[alperaltuntas]

• Fix a Fortran runtime error caused by accessing zk(kSL+1) when the KPP surface layer cell index kSL equals the number of vertical levels
• Guard the post-loop tauEbot normalization: use the cell-center spacing zk(ktmp) - zk(ktmp+1) when ktmp < size(zk), and fall back to the interface spacing zi(ktmp) - zi(ktmp+1) at the bottom cell.

Fixes: #112

[alperaltuntas]

I don't have the permissions to add reviewers, so I am mentioning @mnlevy1981 and @gustavo-marques

[mnlevy1981]

I'd like to talk to Bill L about the best way to fix this. I see that the block of code being modified

! Integrate from top of surface layer cell to Surface layer Depth
    delH     = SLDepth + zi(ktmp)
    sigbot   = CVmix_kpp_params_in%surf_layer_ext
    Gbot     = cvmix_kpp_composite_shape(sigbot)
    TauMAG   = ustar * ustar * Gbot / sigbot
    tauCG    = CempCGm * Gbot *  (taux0 * cosOmega - tauy0 * sinOmega)
    ! tauDG  = sqrt( TauMAG**2 - tauCG**2 ) ! G
    tauDG   = TauMAG                       ! E
    tauxbot  = tauDG * cosOmega  -  tauCG * sinOmega
    tauybot  = tauDG * sinOmega  +  tauCG * cosOmega
    tauEbot  = (tauxbot * delU + tauybot * delV) / (zk(ktmp) - zk(ktmp+1) )

Is similar to what happens in the do loop above it:

! Integrate Both Shear Production Terms from Surface to top of surface layer cell
    do ktmp = 1, kSL-1
      delU     = uE(ktmp) - uE(ktmp+1)
      delV     = vE(ktmp) - vE(ktmp+1)
      Omega_E2x= atan2( delV  , delU )
      cosOmega = cos(Omega_E2x)
      sinOmega = sin(Omega_E2x)

      delH = zi(ktmp) - zi(ktmp+1)
      dbot = dtop + delH
      sigbot = dbot / BLdepth
      Gbot     = cvmix_kpp_composite_shape(sigbot)
      TauMAG   = ustar * ustar * Gbot / sigbot
      tauCG    = CempCGm * Gbot *  (taux0 * cosOmega - tauy0 * sinOmega)
      ! tauDG    = sqrt( TauMAG**2 - tauCG**2 ) ! G
      tauDG    = TauMAG                       ! E
      tauxbot  = tauDG * cosOmega  -  tauCG * sinOmega
      tauybot  = tauDG * sinOmega  +  tauCG * cosOmega
      tauEbot  = (tauxbot * delU + tauybot * delV) / (zk(ktmp) - zk(ktmp+1) )

So I just want to make sure we are computing all of those terms correctly when kSl == size(zk)