jm-c
commented
2 years ago I don't think that this is a bug. To evaluate the stratification (N^2) one needs to consider the vertical derivative of density at constant pressure; this is done by using the same Kref in the 2 calculations of density (the 2 calls to FIND_RHO_2D). Note that this is similar to what is done in do_oceanic_phys.F to compute "sigmaR" (and in fact we could and probably should use "sigmaR" in pp81_ri_number.F to improve the efficiency of pkg/pp81).
Hey, I've been trying to use the pp81 to parameterize internal waves, but I've been getting the maximum viscosity in a lot of places where I shouldn't have, mostly the second layer in my 3d model. After debugging and seeing a lot of places where the density differences between layers shouldn't be inverted. I double checked by calculating the density using temperature and salt from the state variables just to make sure this is a density calculation issue.
I noticed the following lines in the file pp81_ri_number.F:
CALL FIND_RHO_2D( I iMin, iMax, jMin, jMax, K, I theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj), O rhoKm1, I Km1, bi, bj, myThid )CALL FIND_RHO_2D( I iMin, iMax, jMin, jMax, K, I theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj), O rhoK, I K, bi, bj, myThid )Looking at the function FIND_RHO_2D, it says that the last K value (4th variable from the end) is just there for debugging reasons, and the first K (5th variable from the start) is what is used as Kref. Maybe the bug here is to switch the first K to become Km1? I tried this and it has lowered all the viscosity to background levels. I'm not sure yet if good or bad but definitely different. Since pp81 uses FIND_RHO_2D and not FIND_RHO_SCALAR I thought that this problem may be different than the previous opened issue.
Thanks for your hard work and any insights!