C------------------------------------------------------------------------------C
C   R.J.Purser, National Centers for Environmental Prediction, Washington D.C. C
C   wd23jp@sun1.wwb.noaa.gov                                          1996     C
C                   SUBROUTINE XMTOXC                                          C
C    Transform from map coordinates to cartesian coordinates on unit-sphere    C
C                                                                              C
C --> XM:      coordinates (2 components) in [-1.,1.] in specified map panel   C
C <-- XC:      cartesian coordinates (3 components) of corresponding point     C
C <-- DXCDXM:  Jacobian of the transformation (a (3*2)-matrix)                 C
C --> IPAN:    index of the specified map panel (between 1 and 6)              C
C------------------------------------------------------------------------------C
      SUBROUTINE XMTOXC(XM,XC,DXCDXM,IPAN)
      PARAMETER(N=64,NP=N+1)
ccc      PARAMETER(N=128,NP=N+1)  !!!!!!DRAGAN
      COMMON/QPAN/X(3,-1:NP,-1:NP)
      DIMENSION XM(2),XC(3),DXCDXM(3,2)
     *,RPAN(3,3,6)
      DATA RPAN/0.,1.,0.,  0.,0.,1.,  1., 0.,0.
     *        ,-1.,0.,0.,  0.,0.,1.,  0., 1.,0.
     *        ,0.,-1.,0.,  0.,0.,1., -1., 0.,0.
     *        ,1., 0.,0.,  0.,0.,1.,  0.,-1.,0.
     *        ,0., 1.,0., -1.,0.,0.,  0., 0.,1.
     *        ,0., 1.,0.,  1.,0.,0.,  0.,0.,-1./

      NV=N+3
      NM=N-1
      AX=ABS(XM(1))*N
      AY=ABS(XM(2))*N
      IX2=AX
      IF(IX2.GE.N)IX2=NM
      IX1=IX2-1
      RX2=AX-IX2
      RX3=RX2-1.
      IY2=AY
      IF(IY2.GE.N)IY2=NM
      IY1=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
C  ACCUMULATE FOUR CONTRIBUTIONS FROM THE 3*3-TEMPLATE INTERPOLATIONS
C  CENTERED AT THIS ELEMENT'S CORNERS...
      CALL QIN3(X(1,IX1,IY1),RX2,RY2,W22,RY3,RX3
     *,XC,DXCDXM(1,1),DXCDXM(1,2),NV)
      CALL QIN3(X(1,IX2,IY1),RX3,RY2,W32,-RY3,-RX2
     *,XC,DXCDXM(1,1),DXCDXM(1,2),NV)
      CALL QIN3(X(1,IX1,IY2),RX2,RY3,W23,-RY2,-RX3
     *,XC,DXCDXM(1,1),DXCDXM(1,2),NV)
      IF(IX2.NE.NM.OR.IY2.NE.NM)THEN
        CALL QIN3(X(1,IX2,IY2),RX3,RY3,W33,RY2,RX2
     *  ,XC,DXCDXM(1,1),DXCDXM(1,2),NV)
      ELSE
C  ... BUT INVOKE SPECIAL "FRACTIONAL POWER" CODE WHEN ONE OF THESE ELEMENT-
C  CORNER POINTS IS THE CORNER OF THE MAP PANEL ITSELF:
        CALL FIN3(X(1,IX2,IY2),RX3,RY3,W33,RY2,RX2
     *  ,XC,DXCDXM(1,1),DXCDXM(1,2),NV)
      ENDIF

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

C  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,IPAN),XC)
      CALL ROTATE(RPAN(1,1,IPAN),DXCDXM(1,1))
      CALL ROTATE(RPAN(1,1,IPAN),DXCDXM(1,2))

      RETURN
      END

C------------------------------------------------------------------------------C
C   R.J.Purser, National Centers for Environmental Prediction, Washington D.C. C
C   wd23jp@sun1.wwb.noaa.gov                                          1996     C
C                   SUBROUTINE XCTOXM                                          C
C    Transform from map coordinates to cartesian coordinates on unit-sphere    C
C                                                                              C
C --> XC:      cartesian coordinates (3) of point on unit sphere               C
C <-- XM:      map coordinates (2) of this point                               C
C <-- DXMDXC:  generalized Jacobian of the transformation (a (2*3)-matrix)     C
C <-- IPAN:    index of the map panel (between 1 and 6) containing the point   C
C------------------------------------------------------------------------------C
      SUBROUTINE XCTOXM(XC,XM,DXMDXC,IPAN)
      DIMENSION XC(3),XM(2),DXMDXC(2,3),DXCDXM(3,2),G2(2,2),XCT(3)
     *,RPAN(3,3,6)
      DATA RPAN/0.,1.,0.,  0.,0.,1.,  1., 0.,0.
     *        ,-1.,0.,0.,  0.,0.,1.,  0., 1.,0.
     *        ,0.,-1.,0.,  0.,0.,1., -1., 0.,0.
     *        ,1., 0.,0.,  0.,0.,1.,  0.,-1.,0.
     *        ,0., 1.,0., -1.,0.,0.,  0., 0.,1.
     *        ,0., 1.,0.,  1.,0.,0.,  0.,0.,-1./
      DATA NIT/20/,SMAX/1.E-12/
      AXC=ABS(XC(1))
      AYC=ABS(XC(2))
      AZC=ABS(XC(3))
      IF(AXC.GT.AYC)THEN
       IF(AXC.GT.AZC)THEN
        IF(XC(1).GT.0.)GOTO 1    !        |X|
                       GOTO 3    !      biggest
       ENDIF
      ELSEIF(AYC.GT.AZC)THEN
        IF(XC(2).GT.0.)GOTO 2    !        |Y|
                       GOTO 4    !      biggest
      ENDIF
        IF(XC(3).GT.0.)GOTO 5    !        |Z|
                       GOTO 6    !      biggest

1     IPAN=1
      GOTO 7
2     IPAN=2
      GOTO 7
3     IPAN=3
      GOTO 7
4     IPAN=4
      GOTO 7
5     IPAN=5
      GOTO 7
6     IPAN=6
7     CALL MULTM(RPAN(1,1,IPAN),XC,XCT,3,3,1,3,3,3)
      D=1./XCT(3)
      XM(1)=XCT(1)*D
      XM(2)=XCT(2)*D
      CALL XMTOXC(XM,XCT,DXCDXM,IPAN)
      DO IT=1,NIT
       CALL MULTM(DXCDXM,DXCDXM,G2,2,3,2,3,3,2)
       DO J=1,3
       DO I=1,2
        DXMDXC(I,J)=DXCDXM(J,I)
       ENDDO
       ENDDO
       DET=G2(1,1)*G2(2,2)-G2(1,2)*G2(2,1)
       IF(ABS(DET).LE.1.E-14)THEN
C  DON'T TRY TO GET JACOBIAN AT THE CORNER SINGULARITIES WHERE DET VANISHES!!!
        CALL ZERM(DXMDXC,2,3,2)
        RETURN
       ENDIF
       CALL LINVMM(G2,DXMDXC,2,3,2,2)
       S=0.
       DO I=1,3
        XCT(I)=XCT(I)-XC(I)
        S=S+XCT(I)**2
       ENDDO
       CALL MSBMM(DXMDXC,XCT,XM,2,3,1,2,3,2)
       IF(S.LE.SMAX)RETURN
       CALL XMTOXC(XM,XCT,DXCDXM,IPAN)
      ENDDO
      STOP' ITERATIONS IN XCTOXM INSUFFICIENT FOR CONVERGENCE'
      END

      SUBROUTINE RELAXA(XO,XB,XC,XD,XF,XG,ORC,AODD)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XO(3),XA(3),XB(3),XC(3),XD(3),XE(3),XF(3),XG(3),XH(3)
     *,JOFI(3)
      DATA JOFI/3,2,1/
      DO I=1,3
       J=JOFI(I)
       XA(I)=XC(J)
       XE(I)=XF(J)
       XH(I)=XG(J)
      ENDDO
      CALL RELAX(XO,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD)
      RETURN
      END

      SUBROUTINE RELAXB(XO,XA,XC,XD,XG,XH,ORC,AODD)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XO(3),XA(3),XB(3),XC(3),XD(3),XE(3),XF(3),XG(3),XH(3)
     *,JOFI(3)
      DATA JOFI/1,3,2/
      DO I=1,3
       J=JOFI(I)
       XB(I)=XD(J)
       XE(I)=XH(J)
       XF(I)=XG(J)
      ENDDO
      CALL RELAX(XO,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD)
      RETURN
      END

      SUBROUTINE RELAXC(XO,XA,XB,XD,XE,XH,ORC,AODD)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XO(3),XA(3),XB(3),XC(3),XD(3),XE(3),XF(3),XG(3),XH(3)
     *,LOFI(3)
      DATA LOFI/-1,1,1/
      DO I=1,3
       L=LOFI(I)
       XC(I)=XA(I)*L
       XF(I)=XE(I)*L
       XG(I)=XH(I)*L
      ENDDO
      CALL RELAX(XO,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD)
      RETURN
      END

      SUBROUTINE RELAXD(XO,XA,XB,XC,XE,XF,ORC,AODD)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XO(3),XA(3),XB(3),XC(3),XD(3),XE(3),XF(3),XG(3),XH(3)
     *,LOFI(3)
      DATA LOFI/1,-1,1/
      DO I=1,3
       L=LOFI(I)
       XD(I)=XB(I)*L
       XG(I)=XF(I)*L
       XH(I)=XE(I)*L
      ENDDO
      CALL RELAX(XO,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD)
      RETURN
      END

      SUBROUTINE RELAX(X,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (DELTA=.01D0,DELTAI=1.D0/DELTA)
      DIMENSION X(3),XA(3),XB(3),XC(3),XD(3),XE(3),XF(3),XG(3),XH(3)
     *,XU(3),XV(3),U(3),V(3)
      SX=0.
      SY=0.
      DO I=1,3
       U(I)=(XA(I)-XC(I))*.5
       V(I)=(XB(I)-XD(I))*.5
       SX=SX+X(I)*U(I)
       SY=SY+X(I)*V(I)
      ENDDO
      DO I=1,3
       U(I)=U(I)-X(I)*SX
       V(I)=V(I)-X(I)*SY
      ENDDO
      SX=0.
      SY=0.
      DO I=1,3
       XU(I)=X(I)+U(I)*DELTA
       XV(I)=X(I)+V(I)*DELTA
       SX=SX+XU(I)**2
       SY=SY+XV(I)**2
      ENDDO
      SX=1./SQRT(SX)
      SY=1./SQRT(SY)
      DO I=1,3
       XU(I)=XU(I)*SX
       XV(I)=XV(I)*SY
      ENDDO
      CALL REL(X ,XA,XB,XC,XD,XE,XF,XG,XH,U,V,AODD,DX,DY)
      CALL REL(XU,XA,XB,XC,XD,XE,XF,XG,XH,U,V,AODD,DXX,DYX)
      CALL REL(XV,XA,XB,XC,XD,XE,XF,XG,XH,U,V,AODD,DXY,DYY)
      DXX=(DXX-DX)*DELTAI
      DXY=(DXY-DY+DYX-DX)*DELTAI*.5
      DYY=(DYY-DY)*DELTAI
      DETI=ORC/(DXX*DYY-DXY*DXY)
      SX=(DYY*DX-DXY*DY)*DETI
      SY=(DXX*DY-DXY*DX)*DETI
      S=0.
      DO I=1,3
       X(I)=X(I)-U(I)*SX-V(I)*SY
       S=S+X(I)**2
      ENDDO
      S=1./SQRT(S)
      DO I=1,3
       X(I)=X(I)*S
      ENDDO
      RETURN
      END

      SUBROUTINE REL(X,XA,XB,XC,XD,XE,XF,XG,XH,U,V,AODD,DX,DY)
C
C           XF  XB  XE
C             QF  QE
C           XC  X   XA
C             QG  QH
C           XG  XD  XH
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION X(3),XA(3),XB(3),XC(3),XD(3),XE(3),XF(3),XG(3),XH(3)
     *,U(3),V(3)
      DX=0.
      DY=0.
      UAB=0.
      UBC=0.
      UCD=0.
      UDA=0.
      UE=0.
      UF=0.
      UG=0.
      UH=0.
      VAB=0.
      VBC=0.
      VCD=0.
      VDA=0.
      VE=0.
      VF=0.
      VG=0.
      VH=0.
      SAB=0.
      SBC=0.
      SCD=0.
      SDA=0.
      SE=0.
      SF=0.
      SG=0.
      SH=0.
      DO I=1,3
       UI=U(I)
       VI=V(I)
       XI=X(I)
       D1=XE(I)-XI
       D2=XB(I)-XA(I)
       D3=XA(I)-XI
       UE=UE+UI*D1
       UAB=UAB+UI*D2
       DX=DX+UI*D3
       VE=VE+VI*D1
       VAB=VAB+VI*D2
       DY=DY+VI*D3
       SAB=SAB+D2*D2
       SE=SE+D1*D2

       D1=XF(I)-XI
       D2=XC(I)-XB(I)
       D3=XB(I)-XI
       UF=UF+UI*D1
       UBC=UBC+UI*D2
       DX=DX+UI*D3
       VF=VF+VI*D1
       VBC=VBC+VI*D2
       DY=DY+VI*D3
       SBC=SBC+D2*D2
       SF=SF+D1*D2

       D1=XG(I)-XI
       D2=XD(I)-XC(I)
       D3=XC(I)-XI
       UG=UG+UI*D1
       UCD=UCD+UI*D2
       DX=DX+UI*D3
       VG=VG+VI*D1
       VCD=VCD+VI*D2
       DY=DY+VI*D3
       SCD=SCD+D2*D2
       SG=SG+D1*D2

       D1=XH(I)-XI
       D2=XA(I)-XD(I)
       D3=XD(I)-XI
       UH=UH+UI*D1
       UDA=UDA+UI*D2
       DX=DX+UI*D3
       VH=VH+VI*D1
       VDA=VDA+VI*D2
       DY=DY+VI*D3
       SDA=SDA+D2*D2
       SH=SH+D1*D2
      ENDDO
      DX=DX+AODD*(UE*SAB+UF*SBC+UG*SCD+UH*SDA
     *           -UAB*SE-UBC*SF-UCD*SG-UDA*SH)
      DY=DY+AODD*(VE*SAB+VF*SBC+VG*SCD+VH*SDA
     *           -VAB*SE-VBC*SF-VCD*SG-VDA*SH)
      RETURN
      END

      SUBROUTINE QIN3(V,RX2,RY2,W,DWX,DWY,VC,DVCX,DVCY,NV)
C  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
C  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,NV)
C  APPLY FRACTIONAL-POWER TECHNIQUE IN CORNER ELEMENT INTERPOLATION
      IMPLICIT COMPLEX(C)
      COMMON/CORNEL/C3O8,R4O3,DCO(7,7)
      DIMENSION V(3,NV,*),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,2,2)
        DVCX(I)=DVCX(I)+DWX*V(I,2,2)
        DVCY(I)=DVCY(I)+DWY*V(I,2,2)
       ENDDO
       RETURN
      ENDIF
!      PRINT'('' ENTERING NONTRIVIAL FIN3'')'

      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,2,2)
       VN(2,I)=V(I,3,2)
       VN(3,I)=V(I,1,3)
       VN(4,I)=V(I,1,2)
       VN(5,I)=V(I,1,1)
       VN(6,I)=V(I,2,1)
       VN(7,I)=V(I,3,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)
      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

