bug in calculating average potential vorticity

[jpringle]

Dear all --

I think I have found a bug in how the potential vorticity calculation is computed. I was trying to make the output match my expectations for a motionless, stratified f-plane ocean, and found this (I think) error.

In set_avg.F, the function that calculates the potential vorticity is called as (make sure to scroll to the bottom of the code)

Code: Select all

        CALL vorticity_tile (ng, tile,                                  &
     &                       LBi, UBi, LBj, UBj,                        &
     &                       IminS, ImaxS, JminS, JmaxS,                &
# ifdef SOLVE3D
     &                       Kout, Nout,                                &
# else
     &                       Kout,                                      &
# endif
# ifdef MASKING
     &                       GRID(ng) % pmask,                          &
     &                       GRID(ng) % umask,                          &
     &                       GRID(ng) % vmask,                          &
# endif
     &                       GRID(ng) % fomn,                           &       !BUG?, note it is passing fomn
     ...

However, in the definition of the vorticity_tile routine in vorticity.F, the argument list is

Code: Select all

      SUBROUTINE vorticity_tile (ng, tile,                              &
     &                           LBi, UBi, LBj, UBj,                    &
     &                           IminS, ImaxS, JminS, JmaxS,            &
# ifdef SOLVE3D
     &                           kout, nout,                            &
# else
     &                           kout,                                  &
# endif
# ifdef MASKING
     &                           pmask, umask, vmask,                   &
# endif
     &                           f, h, om_u, on_v, pm, pn,              &    !note that it is called f here
...

I think that in set_avg.F, it "f" should be passed in instead of "fomn" because in the vorticity calculation "f" is used in calculations such as

Code: Select all

              cff=0.0625_r8*                                            &
     &            (pm(i-1,j-1)+pm(i-1,j)+pm(i,j-1)+pm(i,j))*            &
     &            (pn(i-1,j-1)+pn(i-1,j)+pn(i,j-1)+pn(i,j))
              fomn_p=0.25_r8*(f(i-1,j-1)+f(i-1,j)+f(i,j-1)+f(i,j))/cff

which suggest that vorticity does really want to receive f, for otherwise the calculation of fomn_p (fomn on psi) is wrong... If this change is made, the calculation matches pencil and paper estimates.

Cheers,
Jamie