subroutine pqk(almd,aphd,inside,p,q,k) c c *** Routine to transform alm, aph (longitude, latitude) into c transformed (rotated) longitude, latitude, check whether the c values obtained are within the model region, and, if they c are, calculate values of p,q (coordinates within a square c formed by connecting four neighboring height points) and k. c c *** Grid constants: c c wbd,sbd - Western and southern boundaries on ll or tll grid c in degrees c tlm0d,tph0d - Angles of rotation of the ll coordinate system c in the direction of lambda and phi respectively c in order to obtain the tll coordinate system c dlmd,dphd - Mesh sides in degrees c implicit none c include 'ecommons.h' c real*4 almd,aphd,p,q . ,rim1,rdlmd,rdphd . ,ebd,nbd . ,rlmd,srlm,crlm,saph,caph,cc . ,anum,denom . ,x1,x2,y1,y2 c integer*4 k,i2,j2,i3,j3,jr CTO REAL*4,PARAMETER :: d2r=.017453293 CTO c logical inside c_______________________________________________________________________________ c c *** Define constants. c rim1=im1 rdlmd=1./dlmd rdphd=1./dphd ebd=wbd+2*(im1)*dlmd nbd=sbd+ (jm1)*dphd c c *** almd, aphd is lon, lat in degrees. c *** Transform into rotated longitude, latitude. c rlmd=almd-tlm0d srlm=sind(rlmd) crlm=cosd(rlmd) saph=sind(aphd) caph=cosd(aphd) cc=crlm*caph c anum=srlm*caph denom=ctph0*cc+stph0*saph c CTO almd=atan2d(anum,denom) CTO aphd=asind(ctph0*saph-stph0*cc) almd=atan2(anum,denom) almd=almd*(1./d2r) aphd=asin(ctph0*saph-stph0*cc) aphd=aphd*(1./d2r) c *** Is almd, aphd inside the (one line reduced) model domain? c if (almd .le. wbd .or. almd .ge. ebd .or. . aphd .le. sbd .or. aphd .ge. nbd) then inside=.false. return endif inside=.true. c c *** x1, y1 is a coordinate system with dlmd, dphd as length units. c x1=(almd-wbd)*rdlmd y1=(aphd-sbd)*rdphd c c *** x2, y2 rotated for +45 deg. & translated for im1. c x2=0.5*( x1+y1) y2=0.5*(-x1+y1)+rim1 c c *** i2, j2 are coordinates of southern points of grid boxes. c i2=int(x2) j2=int(y2) c c *** Remaining parameters need to know the position c of the transformed point within the model region. c p=x2-i2 q=y2-j2 c c *** Index k corresponds to i2, j2. c jr=j2-im1 i3=i2-jr j3=i2+jr k=min(imjm,max(1,j3*im-j3/2+(i3+2)/2)) c return end