[ncl-talk] calculating horizontal vorticity in x and y dim
BasitAli Khan
BasitAli.Khan at kaust.edu.sa
Sat Nov 29 15:07:22 MST 2014
Hi,
I want to calculate the first two terms of vorticity equation separately.
VOR = i(dw/dy – dv/dp) + j(du/dp – dw/dx) + k(dv/dx – du/dy)
I want to calculate the first and second term separately. I wrote the fortran code of x and y horizontal vorticity by utilizing the potential vorticity code of wrf_pvo.f.
Both codes runs fine but I am suspicious about the accuracy of my code. It appears that the sign of the vorticity is changed and the code calculates +ve vorticity as –ve and vice a versa.
I am copying the fortran code for both terms that I call in ncl using WRAPIT and will greatly appreciate if someone could have a look and advice if I have made any mistake. I am not cc'ing it to WRFHELP since talking to wrfhelp for such problems is generally waste of time.
For the first term (dw/dy - dv/dp) I wrote the following code:
C NCLFORTSTART
SUBROUTINE xvort(XV,U,V,W,PRS,MSFU,MSFV,MSFT,COR,DX,DY,
+ NX,NY,NZ,NXP1,NYP1)
IMPLICIT NONE
INTEGER NX,NY,NZ,NXP1,NYP1
DOUBLE PRECISION U(NXP1,NY,NZ),V(NX,NYP1,NZ), W(NX,NY,NZ)
DOUBLE PRECISION PRS(NX,NY,NZ), XV(NX,NY,NZ)
DOUBLE PRECISION MSFU(NXP1,NY),MSFV(NX,NYP1),MSFT(NX,NY)
DOUBLE PRECISION COR(NX,NY)
DOUBLE PRECISION DX,DY
C NCLEND
INTEGER KP1,KM1,JP1,JM1,IP1,IM1,I,J,K
DOUBLE PRECISION DSY,DSX,DP,DWDY,DVDP
DOUBLE PRECISION XVOR,DTHDX,DTHDY,MM
DO K = 1,NZ
KP1 = MIN(K+1,NZ)
KM1 = MAX(K-1,1)
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
MM = MSFT(I,J)*MSFT(I,J)
DWDY = 0.5D0* (W(I,JP1,K)/MSFU(I,JP1)+
+ W(I+1,JP1,K)/MSFU(I+1,JP1)-
+ W(I,JM1,K)/MSFU(I,JM1)-
+ W(I+1,JM1,K)/MSFU(I+1,JM1))/DSY*MM
DP = PRS(I,J,KP1) - PRS(I,J,KM1)
DVDP = 0.5D0* (V(I,J,KP1)+V(I,J+1,KP1)-V(I,J,KM1)-
+ V(I,J+1,KM1))/DP
XVOR = DWDY - DVDP
XV(I,J,K) = XVOR*1.D2
END DO
END DO
END DO
RETURN
END
________________________________
For the second term (du/dp – dw/dx), I wrote the following code:
C NCLFORTSTART
SUBROUTINE yvort(YV,U,V,W,PRS,MSFU,MSFV,MSFT,COR,DX,DY,
+ NX,NY,NZ,NXP1,NYP1)
IMPLICIT NONE
INTEGER NX,NY,NZ,NXP1,NYP1
DOUBLE PRECISION U(NXP1,NY,NZ),V(NX,NYP1,NZ), W(NX,NY,NZ)
DOUBLE PRECISION PRS(NX,NY,NZ), YV(NX,NY,NZ)
DOUBLE PRECISION MSFU(NXP1,NY),MSFV(NX,NYP1),MSFT(NX,NY)
DOUBLE PRECISION COR(NX,NY)
DOUBLE PRECISION DX,DY
C NCLEND
INTEGER KP1,KM1,JP1,JM1,IP1,IM1,I,J,K
DOUBLE PRECISION DSY,DSX,DP,DWDX,DUDP
DOUBLE PRECISION YVOR,DTHDX,DTHDY,MM
print*,'nx,ny,nz,nxp1,nyp1'
print*,DX,DY,NX,NY,NZ
DO K = 1,NZ
KP1 = MIN(K+1,NZ)
KM1 = MAX(K-1,1)
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
MM = MSFT(I,J)*MSFT(I,J)
DWDX = 0.5D0* (W(IP1,J,K)/MSFV(IP1,J)+
+ W(IP1,J+1,K)/MSFV(IP1,J+1)-
+ W(IM1,J,K)/MSFV(IM1,J)-
+ W(IM1,J+1,K)/MSFV(IM1,J+1))/DSX*MM
DP = PRS(I,J,KP1) - PRS(I,J,KM1)
DUDP = 0.5D0* (U(I,J,KP1)+U(I+1,J,KP1)-U(I,J,KM1)-
+ U(I+1,J,KM1))/DP
YVOR = DUDP - DWDX
YV(I,J,K) = YVOR*1.D2
END DO
END DO
END DO
RETURN
END
Best regards and thanks in advance for the help
b:)
----
Basit A. Khan, Ph.D.
Postdoctoral Research Fellow
Division of Physical Sciences & Engineering
Office# 3204, Level 3, Building 1,
King Abdullah University of Science & Technology
4700 King Abdullah Blvd, Box 2753, Thuwal 23955 –6900,
Kingdom of Saudi Arabia.
Office: +966(0)12 808 0276, Mobile: +966(0)5 0860 3617
E-mail: basitali.khan at kaust.edu.sa<mailto:basitali.khan at kaust.edu.sa>
Skype name: basit.a.khan
________________________________
This message and its contents including attachments are intended solely for the original recipient. If you are not the intended recipient or have received this message in error, please notify me immediately and delete this message from your computer system. Any unauthorized use or distribution is prohibited. Please consider the environment before printing this email.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20141129/9b3794f1/attachment.html
More information about the ncl-talk
mailing list