module proc1
	contains
!------------------------------------------------------------------------------C
!   R.J.Purser, National Meteorological Center, Washington D.C.  1995          C
!                   SUBROUTINE INMAP                                         C
!  Initialize the rotation matrix ROT3 needed to transform standard            C
!  earth-centered cartesian components to the alternative cartesian frame      C
!  oriented so as to put geographical point (ALT0,ALN0) on the projection      C
!  axis.                                                                       C
!------------------------------------------------------------------------------C
  SUBROUTINE INMAP(ALN0,ALT0,ROT3)
      DIMENSION  ROT3(3,3)
      PI=4*ATAN(1.)
      DR=PI/180.
      BLON0=DR*ALN0
      BLAT0=DR*ALT0
      CLON0=COS(BLON0)
      SLON0=SIN(BLON0)
      CLAT0=COS(BLAT0)
      SLAT0=SIN(BLAT0)
      ROT3(1,1)=-SLON0
      ROT3(1,2)=CLON0
      ROT3(1,3)=0.
      ROT3(2,1)=-SLAT0*CLON0
      ROT3(2,2)=-SLAT0*SLON0
      ROT3(2,3)=CLAT0
      ROT3(3,1)=CLAT0*CLON0
      ROT3(3,2)=CLAT0*SLON0
      ROT3(3,3)=SLAT0
  END SUBROUTINE INMAP

  SUBROUTINE NMAP(ROT3,X,T)
	 use mulmm_m
    DIMENSION X(3),T(3),U(3), ROT3(3,3)
      CALL MULMM(ROT3,X,U,3,3,1,3,3,3)
      DO I=1,3
      T(I)=U(I)
      ENDDO
  END SUBROUTINE NMAP

  SUBROUTINE NMAPT(ROT3,X,T)
  use mulmm_m
    DIMENSION X(3),T(3),U(3), ROT3(3,3)
      CALL MULTM(ROT3,X,U,3,3,1,3,3,3)
      DO I=1,3
      T(I)=U(I)
      ENDDO
   END SUBROUTINE NMAPT
end module proc1

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!
!                         WD23JP
!                                                         *****************
!                                                         *   SOCT3.FOR   *
!                                                         *  PURSER 1995  *
!                                                         *****************

!------------------------------------------------------------------------------C
!   R.J.Purser, NCEP, Washington D.C.   1996                                   C
!                   SUBROUTINE INSOCT2                                         C
!   Initialize the constants of common/CSTOCT/ used to describe grid and to    C
!   perform transformations between map coordinates and the sphere or disk     C
!                                                                              C
! --> FLON0,FLAT0 longitude and latitude of center of octagon-1                C
! --> SCEN      scale-enhancement parameter = radius of image of octagon-1 in  C
!               a stereographic projection centered on (FLAT0,FLON0) that maps C
!               the concentric hemisphere to the unit disk                     C
! --> TFILE     character*39 name of file from which table of location are got C
!------------------------------------------------------------------------------C
      SUBROUTINE INSOCT2(flon0,flat0,SCEN,TFILE)
		use proc1
      CHARACTER*39 TFILE
      PARAMETER(NGEN=6,NRAT1=1,NRAT2=3,N=NRAT2*2**NGEN,NP=N+1,NAMAG=1)
		real, dimension (3,3) :: rotm, rotm0
		real :: aln0, alt0
      COMMON/QPAN/X(3,-1:NP,-1:NP)
      COMMON/CSTOCT/MAG,M1G,ROTM,SC,SCI,A
!test
       print *,'flon0=',flon0
       print *,'flat0=',flat0
       print *,'scen =',scen
       print *,' '
       print *,'tfile=',tfile
!test
		 aln0 = flon0
		 alt0 = flat0

      CALL INMAP(aln0,alt0,ROTM0)
		ROTM=ROTM0
      MAG=NRAT1*NAMAG
      M1G=NRAT2*NAMAG
      SC=SCEN
      SCI=1./SC
      A=FLOAT(NRAT1)/NRAT2
      OPEN(UNIT=9,FILE=TFILE,STATUS='UNKNOWN',FORM='unformatted')
      READ(9)X
      CLOSE(UNIT=9)
      CALL INFIN3
      RETURN
      END

!------------------------------------------------------------------------------C
!   R.J.Purser, National Centers for Environmental Prediction, Washington D.C. C
!   wd23jp@sun1.wwb.noaa.gov                                          1996     C
!                   SUBROUTINE XMTOXC                                          C
!    Transform from map coordinates to cartesian coordinates on unit-sphere    C
!                                                                              C
! --> XM:      coordinates (2 components) in [-1.,1.] in specified map panel   C
! <-- XC:      cartesian coordinates (3 components) of corresponding point     C
! <-- DXCDXM:  Jacobian of the transformation (a (3*2)-matrix)                 C
! --> IOCT:    index of the specified octagon (1 or 2)                         C
!------------------------------------------------------------------------------C
      SUBROUTINE XMTOXC(XM,XC,DXCDXM,IOCT)
      PARAMETER(NGEN=6,NRAT1=1,NRAT2=3,N=NRAT2*2**NGEN,NP=N+1,NA=NRAT1*2**NGEN)
      COMMON/QPAN/X(3,-1:NP,-1:NP)
      DIMENSION XM(2),XC(3),DXCDXM(3,2),RPAN(3,3,2)
      DATA RPAN/1.,0.,0.,  0.,1.,0.,   0.,0.,1. &
              ,-1.,0.,0.,  0.,1.,0.,   0.,0.,-1./
      NM=N-1
      AX=ABS(XM(1))*N
      AY=ABS(XM(2))*N
      IX2=AX
      IF(IX2.GE.N)IX2=NM
      IX1=IX2-1
      IX3=IX2+1
      RX2=AX-IX2
      RX3=RX2-1.
      IY2=AY
      IF(IY2.GE.N)IY2=NM
      IY1=IY2-1
      IY3=IY2+1
      RY2=AY-IY2
      RY3=RY2-1.
      W22=RX3*RY3
      W32=-RX2*RY3
      W23=-RX3*RY2
      W33=RX2*RY2
      DO I=1,3
       XC(I)=0.
       DXCDXM(I,1)=0.
       DXCDXM(I,2)=0.
      ENDDO
!  ACCUMULATE FOUR CONTRIBUTIONS FROM THE 3*3-TEMPLATE INTERPOLATIONS
!  CENTERED AT THIS ELEMENT'S CORNERS...
      CALL GIN3(X,RX2,RY2,W22,RY3,RX3,XC,DXCDXM(1,1),DXCDXM(1,2),IX2,IY2,NA,N)
      CALL GIN3(X,RX3,RY2,W32,-RY3,-RX2,XC,DXCDXM(1,1),DXCDXM(1,2),IX3,IY2,NA,N)
      CALL GIN3(X,RX2,RY3,W23,-RY2,-RX3,XC,DXCDXM(1,1),DXCDXM(1,2),IX2,IY3,NA,N)
      CALL GIN3(X,RX3,RY3,W33,RY2,RX2,XC,DXCDXM(1,1),DXCDXM(1,2),IX3,IY3,NA,N)

      DO I=1,3
       DXCDXM(I,1)=DXCDXM(I,1)*N
       DXCDXM(I,2)=DXCDXM(I,2)*N
      ENDDO

!  NORMALIZATION AND ORTHOGONALIZATION
      T1=0.
      T2=0.
      S=0.
      DO I=1,3
       T1=T1+XC(I)*DXCDXM(I,1)
       T2=T2+XC(I)*DXCDXM(I,2)
       S=S+XC(I)**2
      ENDDO
      S=1./SQRT(S)
      T1=T1*S
      T2=T2*S
      DO I=1,3
       XC(I)=XC(I)*S
       DXCDXM(I,1)=DXCDXM(I,1)-T1*XC(I)
       DXCDXM(I,2)=DXCDXM(I,2)-T2*XC(I)
      ENDDO
      IF(XM(1).LT.0.)THEN
        XC(1)=-XC(1)
        DXCDXM(2,1)=-DXCDXM(2,1)
        DXCDXM(3,1)=-DXCDXM(3,1)
        DXCDXM(1,2)=-DXCDXM(1,2)
      ENDIF
      IF(XM(2).LT.0.)THEN
        XC(2)=-XC(2)
        DXCDXM(2,1)=-DXCDXM(2,1)
        DXCDXM(1,2)=-DXCDXM(1,2)
        DXCDXM(3,2)=-DXCDXM(3,2)
      ENDIF
      CALL ROTATE(RPAN(1,1,IOCT),XC)
      CALL ROTATE(RPAN(1,1,IOCT),DXCDXM(1,1))
      CALL ROTATE(RPAN(1,1,IOCT),DXCDXM(1,2))
      RETURN
      END

      SUBROUTINE GIN3(V,XZ,YZ,W,DWX,DWY,VC,DVCX,DVCY,IX,IY,NA,N)
      DIMENSION V(3,-1:N+1,-1:*),VC(3),DVCX(3),DVCY(3),VT(3,-1:1,-1:1)
      NAP=NA+1
      IF(IX.EQ.N)THEN
       IF(IY.EQ.NA)THEN
        CALL FIN3(V(1,IX-1,IY-1),XZ,YZ,W,DWX,DWY,VC,DVCX,DVCY,N)
        RETURN
       ELSEIF(IY.EQ.NAP)THEN
        DO I=1,3
         DO JY=-1,1
         DO JX=-1,1
          VT(I,JX,JY)=V(I,IX+JX,IY+JY)
         ENDDO
         ENDDO
         VT(I,1,-1)=V(I,N,NA-1)
        ENDDO
        CALL QIN3(VT,XZ,YZ,W,DWX,DWY,VC,DVCX,DVCY,3)
        RETURN
       ENDIF
      ENDIF
      IF(IY.EQ.N)THEN
       IF(IX.EQ.NA)THEN
        CALL FIN3(V(1,IX-1,IY-1),XZ,YZ,W,DWX,DWY,VC,DVCX,DVCY,N)
        RETURN
       ELSEIF(IX.EQ.NAP)THEN
        DO I=1,3
         DO JY=-1,1
         DO JX=-1,1
          VT(I,JX,JY)=V(I,IX+JX,IY+JY)
         ENDDO
         ENDDO
         VT(I,-1,1)=V(I,NA-1,N)
        ENDDO
        CALL QIN3(VT,XZ,YZ,W,DWX,DWY,VC,DVCX,DVCY,3)
        RETURN
       ENDIF
      ENDIF
      CALL QIN3(V(1,IX-1,IY-1),XZ,YZ,W,DWX,DWY,VC,DVCX,DVCY,N+3)
      RETURN
      END

      SUBROUTINE QIN3(V,RX2,RY2,W,DWX,DWY,VC,DVCX,DVCY,NV)
!  APPLY QUADRATIC INTERPOLATION AT NON-CORNER ELEMENT OF TABLE
      DIMENSION V(3,NV,*),VC(3),DVCX(3),DVCY(3),PX(3),PY(3), &
       DPX(3),DPY(3),YV(3),YDVX(3),YDVY(3),XV(3,3),XDVX(3,3)
      RX1=RX2+1.
      RX3=RX2-1.
      PX(1)= (RX2*RX3)*.5
      PX(2)=-(RX1*RX3)
      PX(3)= (RX1*RX2)*.5
      DPX(1)= (RX2+RX3)*.5
      DPX(2)=-(RX1+RX3)
      DPX(3)= (RX1+RX2)*.5
      RY1=RY2+1.
      RY3=RY2-1.
      PY(1)= (RY2*RY3)*.5
      PY(2)=-(RY1*RY3)
      PY(3)= (RY1*RY2)*.5
      DPY(1)= (RY2+RY3)*.5
      DPY(2)=-(RY1+RY3)
      DPY(3)= (RY1+RY2)*.5
      DO I=1,3
       YV(I)=0.
       YDVX(I)=0.
       YDVY(I)=0.
      ENDDO
      DO IY=1,3
       DO I=1,3
        XV(I,IY)=0.
        XDVX(I,IY)=0.
       ENDDO
       DO IX=1,3
        DO I=1,3
         XV(I,IY)=XV(I,IY)+PX(IX)*V(I,IX,IY)
         XDVX(I,IY)=XDVX(I,IY)+DPX(IX)*V(I,IX,IY)
        ENDDO
       ENDDO
       DO I=1,3
        YV(I)=YV(I)+PY(IY)*XV(I,IY)
        YDVY(I)=YDVY(I)+DPY(IY)*XV(I,IY)
        YDVX(I)=YDVX(I)+PY(IY)*XDVX(I,IY)
       ENDDO
      ENDDO
      DO I=1,3
       VC(I)=VC(I)+W*YV(I)
       DVCX(I)=DVCX(I)+W*YDVX(I)+DWX*YV(I)
       DVCY(I)=DVCY(I)+W*YDVY(I)+DWY*YV(I)
      ENDDO
      RETURN
      END

      SUBROUTINE INFIN3
!  INITIALIZE PARAMETERS NEEDED BY FIN3
      IMPLICIT COMPLEX(C)
      COMMON/CORNEL/C3O8,R4O3,DCO(7,7)
      DIMENSION IPAR(7)
      DATA IPAR/1,1,-1,1,-1,1,1/
      R1O3=1./3
      C3O8=CMPLX(0.,1.)
      C3O8=-1./CSQRT(C3O8)
      R4O3=4./3
      CALL ZERM(DCO,7,7,7)
      RHO=2.**R1O3
      RHOS=RHO*RHO
      S=SQRT(.75)
      SIG=R1O3/(RHO+2.)
      RHOP=RHO+1.
      RHOM=RHO-1.
      DCO(1,1)=1.
      DCO(2,2)=4.*SIG
      DCO(2,3)=SIG*RHOS/2
      DCO(2,4)=-2.*SIG
      DCO(2,5)=-SIG*RHOS
      DCO(3,3)=SIG*S*RHOS
      DCO(3,4)=4.*SIG*S
      DCO(4,1)=-15.*SIG
      DCO(4,2)=4.*SIG*RHOP
      DCO(4,3)=-SIG*(RHOS-1.)
      DCO(4,4)=2.*SIG*(2.-RHO)
      DCO(4,5)=SIG*(2.*RHOS+1.)
      DCO(5,3)=2.*SIG*RHOS*S
      DCO(5,4)=-4.*SIG*RHO*S
      DCO(6,1)=DCO(4,1)
      DCO(6,2)=-4.*SIG*RHOM
      DCO(6,3)=RHOS/6
      DCO(6,4)=2./3
      DCO(6,5)=-SIG*(2.*RHOS-1.)
      DCO(7,1)=-RHO*(2.-RHO)/4.
      DCO(7,2)=SIG*RHO
      DCO(7,3)=-SIG/2
      DCO(7,4)=DCO(7,2)
      DCO(7,5)=DCO(7,3)
      DO I=2,7
       DCO(I,6)=DCO(I,4)*IPAR(I)
       DCO(I,7)=DCO(I,3)*IPAR(I)
      ENDDO
      RETURN
      END

      SUBROUTINE FIN3(V,XZ,YZ,W,DWX,DWY,VC,DVCX,DVCY,N)
		use mulmm_m
!  APPLY FRACTIONAL-POWER TECHNIQUE IN CORNER ELEMENT INTERPOLATION
      IMPLICIT COMPLEX(C)
      COMMON/CORNEL/C3O8,R4O3,DCO(7,7)
      DIMENSION V(3,-1:N+1,-1:*),VC(3),DVCX(3),DVCY(3)
      DIMENSION P(7),PX(7),PY(7),VCT(3),VCTX(3),VCTY(3),VN(7,3),DV(7,3)
      DATA P/7*0./,PX/7*0./PY/7*0./
      IF(XZ.EQ.0..AND.YZ.EQ.0.)THEN
       DO I=1,3
        VC(I)=VC(I)+W*V(I,0,0)
        DVCX(I)=DVCX(I)+DWX*V(I,0,0)
        DVCY(I)=DVCY(I)+DWY*V(I,0,0)
       ENDDO
       RETURN
      ENDIF
      WR4O3=W*R4O3
      CZ=CMPLX(XZ,YZ)
      C=-(CZ*C3O8)**R4O3
      CR=C/CZ
      WR=REAL(CR)*WR4O3
      WQ=AIMAG(CR)*WR4O3
      X=REAL(C)
      Y=AIMAG(C)
      DO I=1,3
       VN(1,I)=V(I,0,0)
       VN(2,I)=V(I,1,0)
       VN(3,I)=V(I,-1,1)
       VN(4,I)=V(I,-1,0)
       VN(5,I)=V(I,-1,-1)
       VN(6,I)=V(I,0,-1)
       VN(7,I)=V(I,1,-1)
       CALL MULMM(DCO,VN(1,I),DV(1,I),7,7,1,7,7,7)
      ENDDO
      P(1)=1.
      P(2)=X
      P(3)=Y
      P(4)=.5*X**2
      P(5)=X*Y
      P(6)=.5*Y**2
      P(7)=X**3-3.*X*Y**2
      PX(2)=1.
      PX(4)=X
      PX(5)=Y
      PX(7)=3.*(X**2-Y**2)
      PY(3)=1.
      PY(5)=X
      PY(6)=Y
      PY(7)=-6*X*Y
      DO I=1,3
       VCT(I)=0.
       VCTX(I)=0.
       VCTY(I)=0.
       DO J=1,7
        VCT(I)=VCT(I)+P(J)*DV(J,I)
        VCTX(I)=VCTX(I)+PX(J)*DV(J,I)
        VCTY(I)=VCTY(I)+PY(J)*DV(J,I)
       ENDDO
       VC(I)=VC(I)+W*VCT(I)
       DVCX(I)=DVCX(I)+WR*VCTX(I)+WQ*VCTY(I)+DWX*VCT(I)
       DVCY(I)=DVCY(I)-WQ*VCTX(I)+WR*VCTY(I)+DWY*VCT(I)
      ENDDO
      RETURN
      END

      SUBROUTINE ROTATE(R3,XX)
		use mulmm_m
      DIMENSION R3(3,3),XX(3),VV(3)
      CALL MULMM(R3,XX,VV,3,3,1,3,3,3)
      DO I=1,3
      XX(I)=VV(I)
      ENDDO
      RETURN
      END