[ncl-talk] calculating horizontal vorticity in x and y dim

[BasitAli Khan]

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