Data in rotated grid coordinates
Ren Diandong
dd_ren at ROSSBY.METR.OU.EDU
Thu Nov 9 15:43:36 EST 2006
Hi,
PLEASE DO NOT HESITATE TO FURTHER CONSULT.
First, I have some comments about your forward code:
1)wrot(i,j)=-999.0 and did not see where else it is used (referenced).
see my inline comments for you to further put this part in...
2) the many sign is unnecessary, it is resulted fron solve atan(x)=y and
if you use the following conventioon, it is unnecessary: lat -pi/2~pi/2
lon (0-2pi)
I use tyhe pgf90 compiler and took liberty changed your forward code a
little bit.
Second, I am not sure of your compiler. So, I also attached my control
script I used here.
Finally, to illustrate the usage, I use a PSE projection (as mentioned
in the first email exchange) as example. But the design philosophy is
identical as the code here for you. I use the PSE projection merely
because I am familiar with it and my general system include this
projection as a member.
Attached figure Panel a is produced use GrADS background map and grid
line setting, using pdef in the *.ctl file. Panel b is produced solely
by calling multiple time of the backward-forward pare and only use GrADS
provided lat and lon (two implicit variables). The grid lines are draw
use draw line command.
Again, if you provide a complete (e.g., wrot() problem) forward code, I
am willing to improve the backward code. Or, you can follow this example
to insert the missing part. If wrot() is ornate, then the attached
routines suffice for general plotting purpose.
Best wishes,
Diandong Ren
OU/SoM
PS. Reference to the following for pgf90
FTN = pgf90
LOADER = pgf90
CC = cc
ZIP = gzip
TAR = tar
RM = rm
TOPDIR = .$(bash pwd)
INCL_DIR = include
DATA_DIR = data
On Thu, 2006-11-09 at 14:06 +0000, Gudrun Nina Petersen wrote:
> wrot(i,j)=-999.0
-------------- next part --------------
program main ()
write(*,*)"finished"
end program main
! *************************************************************************
! *************************************************************************
! *** A program for writing a pdef bilin file for GrADS data with rotated pole.
! *** Based on Jeff Chagnon's um_xy_nc.m that calculates the real latitudes
! *** and longitude variables for the UM in matlab.
! ***
! *** At the moment every thing is hard coded.
! *** ifort -o um_xy_nc um_xy_nc.f
! *************************************************************************
! *************************************************************************
subroutine RotFindLatlonForIj (Nx,Ny,sgn,polelat,polelon,x1,y1,dx,dy, &
LOCATION,LOCATIONij)
IMPLICIT NONE
! INTEGER,parameter :: Nx=326,Ny=420
INTEGER :: Nx,Ny,i,j,sgn
real, dimension (1:nx,1:ny,1:2) :: LOCATION
integer, dimension (1:nx,1:ny,1:2) :: LOCATIONij
! REAL,parameter :: polelat=36.00, polelon=140.00, x1=342.00, y1=-14.00, &
! dx=0.11,dy=0.11, pi=3.14159265,DR=pi/180.0,badvalue=-999.9
REAL :: polelat, polelon, x1, y1,dx,dy
real, parameter :: pi=3.14159265,DR=pi/180.0,badvalue=-999.9
REAL :: sin_phi_pole, cos_phi_pole
REAL :: E_x,E_y
REAL :: arg,arg2,a_phi,term1,term2,a_lambda
REAL, dimension (1:Nx,1:Ny) :: wrot
!-----------------------------------------------------------------------------
! This is to get the right angles if southern hemisphere
!-----------------------------------------------------------------------------
IF (polelat.ge.0) THEN
sin_phi_pole = sin(DR*polelat)
cos_phi_pole = cos(DR*polelat)
ELSE
sin_phi_pole = -sin(DR*polelat)
cos_phi_pole = -cos(DR*polelat)
ENDIF
DO i=1,Nx
DO j=1,Ny
LOCATIONij (i,j,1)=i
LOCATIONij (i,j,2)=j
ENDDO
ENDDO
DO i=1,Nx
DO j=1,Ny
wrot(i,j)=badvalue
E_x=DR*(x1+(LOCATIONij (i,j,1)-1)*dx)
E_y=DR*(y1+(LOCATIONij (i,j,2)-1)*dy)
! Scale longitude to range from -180 to 180
IF (E_x.gt.pi) THEN
E_x = E_x - 2*pi
ELSEIF (E_x.lt.-1*pi) THEN
E_x = E_x + 2*pi
ENDIF
! Compute latitude using equation (4.7)
arg=cos_phi_pole*COS(E_x)*COS(E_y) + SIN(E_y)*sin_phi_pole
arg=MIN(arg, 1.0)
arg=MAX(arg,-1.0)
a_phi=ASIN(arg)
LOCATION(i,j,2)=a_phi/DR
! Compute longitude using equation (4.8)
term1 =(COS(E_x)*COS(E_y)*sin_phi_pole -SIN(E_y)*cos_phi_pole)
term2=COS(a_phi)
IF (ABS(term2).lt.1e-5) THEN
a_lambda=0.0
ELSE
arg2=E_x-2*pi
IF (arg2.gt.0) THEN
sgn=1
ELSEIF (arg2.eq.0) THEN
sgn=0
ELSEIF (arg2.lt.0) THEN
sgn=-1
ELSE
PRINT*, 'What!!!!'
ENDIF
arg=term1/term2;
arg=MIN(arg, 1.0)
arg=MAX(arg,-1.0)
a_lambda=ACOS(arg)/DR
a_lambda=a_lambda*sgn
ENDIF
LOCATION(i,j,1)=a_lambda+polelon-180.0
ENDDO
ENDDO
end subroutine RotFindLatlonForIj
-------------- next part --------------
subroutine RotFindIjForLatlon (Nx,Ny,sgn,polelat,polelon,x1,y1,dx,dy, &
LOCATION,LOCATIONij)
IMPLICIT NONE
INTEGER :: Nx,Ny,i,j,sgn
real, dimension (1:nx,1:ny,1:2) :: LOCATION
integer, dimension (1:nx,1:ny,1:2) :: LOCATIONij
REAL :: polelat, polelon, x1, y1,dx,dy
real, parameter :: pi=3.14159265,DR=pi/180.0,badvalue=-999.9
REAL :: sin_phi_pole, cos_phi_pole
REAL :: E_x,E_y
REAL :: arg,arg2,a_phi,term1,term2,a_lambda
REAL, dimension (1:Nx,1:Ny) :: wrot
!-----------------------------------------------------------------------------
! This is to get the right angles if southern hemisphere
!-----------------------------------------------------------------------------
IF (polelat.ge.0) THEN
sin_phi_pole = sin(DR*polelat)
cos_phi_pole = cos(DR*polelat)
ELSE
sin_phi_pole = -sin(DR*polelat)
cos_phi_pole = -cos(DR*polelat)
ENDIF
DO i=1,Nx
DO j=1,Ny
!edit here for your rotat() arrayu
a_lambda=180.0+LOCATION(i,j,1)-polelon
a_phi=LOCATION(i,j,2)*DR
arg=sin(a_phi)
term2=COS(a_phi)
IF (ABS(term2).lt.1e-5) THEN
a_lambda=0.0
ELSE
a_lambda=a_lambda*sgn
arg=cos(DR*a_lambda)
term1 =term2*arg
ENDIF
E_y=asin( (arg-cos_phi_pole*(arg*cos_phi_pole+term1*sin_phi_pole) )/sin_phi_pole)
E_x=acos((arg*cos_phi_pole+term1*sin_phi_pole)/cos(E_y))
IF (E_x.gt.pi) THEN
E_x = E_x - 2*pi
ELSEIF (E_x.lt.-1*pi) THEN
E_x = E_x + 2*pi
ENDIF
LOCATIONij (i,j,1)=1+int((E_x/DR-x1)/dx)
LOCATIONij (i,j,2)=1+int((E_y/DR-y1)/dy)
ENDDO
ENDDO
end subroutine RotFindIjForLatlon
-------------- next part --------------
A non-text attachment was scrubbed...
Name: forNina.doc
Type: application/vnd.ms-word
Size: 655872 bytes
Desc: not available
Url : http://gradsusr.org/pipermail/gradsusr/attachments/20061109/98d30ff6/attachment.bin
More information about the gradsusr
mailing list