NCL: need help to alculate Relative Vorticity & Divergence

[sahidul]

NCL Users,

I am having WRF forecast files(wrfout_d0* files), saved in 1 hour interval having total 88
files. We want contour plot for "Relative Vorticity" and "Divergence" every hourly.

Is there any function or Script for calculating "Relative Vorticity" and "Divergence" ?

Kindly help me to plot these variables.

Thanking you,

[dbh409]

There are several built-in scripts for NCL that calculate RV and Div, depending on what method you wish to use to calculate it. Here are a couple using finite differencing:
http://www.ncl.ucar.edu/Document/Functions/Built-in/uv2vr_cfd.shtml
http://www.ncl.ucar.edu/Document/Functions/Built-in/uv2vrdvF.shtml

If you want one WRF specific, that's a different story, and I'm not completely sure where to send you, though I know I've seen something like that before. One method for relative vorticity is to use the function wrf_avo and subtract coriolis at every grid point.

[tomk]

I was also wondering about this and emailed wrfhelp about it. They confirmed that for a horizontal vector field with components (A,B) in the x- and y-directions on the lambert grid, then the divergence is given by
m^2 * [ d/dx(A/m) + d/dy(B/m) ]
and the curl is
m^2 * [ d/dx(B/m) - d/dy(A/m) ],
where m is the map factor.

You can use my fortran routines (with wrapit) to compute vorticity and divergence, if you like.

[tomk]

Here's the fortran routine for curl (wrf_curl.f):


c--------------------------------------------------------
C NCLFORTSTART
SUBROUTINE DCOMPUTECURL(CRL,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 CRL(NX,NY,NZ),MSFT(NX,NY)
DOUBLE PRECISION DX,DY
C NCLEND
INTEGER JP1,JM1,IP1,IM1,I,J,K
DOUBLE PRECISION DSY,DSX,DVDX,DUDY
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)
DVDX = (V(IP1,J,K)/MSFT(IP1,J) -
+ V(IM1,J,K)/MSFT(IM1,J))/DSX*MM
DUDY = (U(I,JP1,K)/MSFT(I,JP1) -
+ U(I,JM1,K)/MSFT(I,JM1))/DSY*MM

CRL(I,J,K) = DVDX - DUDY
END DO
END DO
END DO
RETURN
END

[tomk]

...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

[kindahero]

@tomk

thanks for sharing the fortran code.. really helpful

[kindahero]

@tomk
it would be even better., if you can post the associated ncl script.

[ronbeag]

Unless I'm mistaken, you can also get WRF to output divergence as a standard output variable...

In Registry.EM_COMMON around line 1327 (ARW version 3.4.1) there's the variable "div", described as "DIVERGENCE". So you can either edit the Registry file so that this variable is written to history, or use 'Run-Time IO' (see ARW manual page 5-23) while running WRF.

[Basit]

Hi Tom,
I also want to calculate relative vorticity form my WRF model output. I have few questions about the wrf_cur.f

1. Your have written a subroutine, why it is not a 'function' ? because it is calculating vorticity but won't return anything.
2. The first argument in this subroutine is CRL that i believe is the curl(rel. vorticity). But this is what we are trying to calculate. So, why and what should we pass for the first argument when I would call this subroutine from my ncl script?
3. I believe you have already tested this code and it works correctly?

Regards,

BAS.

[MeteoBoy4]

...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