Bug in the computation of buoyancy

[arango]

Originally posted on: 09-10-02 17:25

Dear Fellows,

Scott Dursky brought to my attention today a somewhat serious bug in the computation of the Brunt-Vaisala frequency in "rho_eos.F". This will only affect applications using the nonlinear equation of state.

Code: Select all

bulk = bulk0 + bulk1*z + bulk2*z*z

den = (den1 * bulk) / (bulk + 0.1*z)


            |-----w-----| k+1
            |           |
            |           |      bulk0(k+1)
    den_up  |    t,s    | k+1  bulk0(k+1)*z_w(k),
            |           |      bulk2(k+1)*z_w(k)**2
            |           | 
            |-----w-----| k    z_w(k), bvf(k) 
            |           |
            |           |      bulk0(k)
    den_dn  |    t,s    | k    bulk1(k)*z_w(k), 
            |           |      bulk2(k)*z_w(k)**2
            |           |
            |-----w-----| k-1 

The upper potential density is computed by lowering the water parcel adiabatically to z_w(k):

Code: Select all

      bulk_up=bulk0(i,k+1)-
     &        z_w(i,j,k)*(bulk1(i,k+1)-
     &        bulk2(i,k+1)*z_w(i,j,k))
      den_up=(den1(i,k+1)*bulk_up)/(bulk_up+0.1_r8*z_w(i,j,k))

The lower potential density is computed by rising the water parcel adiabatically to z_w(k):

Code: Select all

      bulk_dn=bulk0(i,k )-
     &        z_w(i,j,k )*(bulk1(i,k )-
     &        bulk2(i,k )*z_w(i,j,k ))
      den_dn=(den1(i,k )*bulk_dn)/(bulk_dn+0.1_r8*z_w(i,j,k ))

Then, the Brunt-Vaisala frequency becomes:

Code: Select all

      bvf(i,j,k)=-g*(den_up-den_dn)/
     &           (0.5_r8*(den_up+den_dn)*
     &           (z_r(i,j,k+1)-z_r(i,j,k)

The computation of "bulk_up" and "den_up" is WRONG in the distributed versions of ROMS. It is very weird that we didn't catch this one before.

At this point I don't know the gravity of this bug. It will definitely alter vertical mixing in near neutral conditions.

Thank you very much to Scott Dursky for finding this crucial bug.

I am sending the new routine to everybody via e-mail.

Best, H