tomk wrote:...and for divergence (wrf.div.f):
c--------------------------------------------------------
C NCLFORTSTART
SUBROUTINE DCOMPUTEDIV(DIV,U,V,MSFT,DX,DY,NX,NY,NZ)
IMPLICIT NONE
INTEGER NX,NY,NZ
DOUBLE PRECISION U(NX,NY,NZ),V(NX,NY,NZ)
DOUBLE PRECISION DIV(NX,NY,NZ),MSFT(NX,NY)
DOUBLE PRECISION DX,DY
C NCLEND
INTEGER JP1,JM1,IP1,IM1,I,J,K
DOUBLE PRECISION DSY,DSX,DUDX,DVDY
DOUBLE PRECISION MM
C Note all data must be on T-pts
DO K = 1,NZ
DO J = 1,NY
JP1 = MIN(J+1,NY)
JM1 = MAX(J-1,1)
DO I = 1,NX
IP1 = MIN(I+1,NX)
IM1 = MAX(I-1,1)
DSX = (IP1-IM1)*DX
DSY = (JP1-JM1)*DY
c Careful with map factors...
MM = MSFT(I,J)*MSFT(I,J)
c DVDX = (V(IP1,J,K)/MSFT(IP1,J) -
c + V(IM1,J,K)/MSFT(IM1,J))/DSX*MM
c DUDY = (U(I,JP1,K)/MSFT(I,JP1) -
c + U(I,JM1,K)/MSFT(I,JM1))/DSY*MM
DUDX = ( U(IP1,J,K)/MSFT(IP1,J) -
+ U(IM1,J,K)/MSFT(IM1,J) )/DSX*MM
DVDY = ( V(I,JP1,K)/MSFT(I,JP1) -
+ V(I,JM1,K)/MSFT(I,JM1) )/DSY*MM
DIV(I,J,K) = DUDX + DVDY
END DO
END DO
END DO
RETURN
END
Actually, in NCL source code (in $NCARG_ROOT/ni/src/lib/nfpfort/wrf_pvo.f), the calculation of absolute vorticity does not require on u,v unstaggering to t grid. But I have some problems understanding why map scale factor outside the parentheses is that of t-grid's.
Following is NCL source code segment:
DO k = 1,nz
DO j = 1,ny
jp1 = MIN(j+1, ny)
jm1 = MAX(j-1, 1)
DO i = 1,nx
ip1 = MIN(i+1, nx)
im1 = MAX(i-1, 1)
! PRINT *,jp1,jm1,ip1,im1
dsx = (ip1 - im1) * dx
dsy = (jp1 - jm1) * dy
mm = msft(i,j)*msft(i,j)
! PRINT *,j,i,u(i,jp1,k),msfu(i,jp1),u(i,jp1,k)/msfu(i,jp1)
dudy = 0.5D0*(u(i,jp1,k)/msfu(i,jp1) + u(i+1,jp1,k)/msfu(i+1,jp1) - &
u(i,jm1,k)/msfu(i,jm1) - u(i+1,jm1,k)/msfu(i+1,jm1))/dsy*mm
dvdx = 0.5D0*(v(ip1,j,k)/msfv(ip1,j) + v(ip1,j+1,k)/msfv(ip1,j+1) - &
v(im1,j,k)/msfv(im1,j) - v(im1,j+1,k)/msfv(im1,j+1))/dsx*mm
avort = dvdx - dudy + cor(i,j)
av(i,j,k) = avort*1.D5
END DO
END DO
END DO