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