***********************************************************************
*     SPHERICAL HARMONICS TRANSFORMATION LIBRARY
***********************************************************************
*     Copyright (C) 2000-2016 GFD Dennou Club. All rights reserved.
*-----------------------------------------------------------------------
      SUBROUTINE SHTLIB

*      REAL Q((JM+1)*(4*JM+5*MM+14)+(MM+1)*(MM+1)+MM+2+6*IM+15)
*      REAL FUN(-JM:JM,M:MM)
*      REAL S((MM+1)*(MM+1)),W(-JM:JM,-MM:MM),G(-IM:IM,-JM:JM)
*      REAL SM(M:MM),WM(-JM:JM),WZ(-JM:JM),WR(-JM:JM),WI(-JM:JM)
*      REAL WJ(-MM:MM),GJ(-IM:IM)
*      REAL A((MM+1)*(MM+1)),B((MM+1)*(MM+1))
      REAL Q(*),FUN(*),S(*),W(*),G(*),SM(*),WM(*),WZ(*),WR(*),WI(*)
      REAL WJ(*),GJ(*),A(*),B(*)
      SAVE

      RETURN
*----------------------------------------------------------------------
*     INITIALIZATION
*----------------------------------------------------------------------
      ENTRY SHTINT(MM,JM,IM,Q)

      IF(IM .LT. MM+1) THEN
        CALL MSGDMP('E','SHTLIB','IM MUST BE IM >= MM+1')
      END IF
      IF(JM .LT. (MM+1)/2) THEN
        CALL MSGDMP('E','SHTLIB','JM MUST BE JM >= (MM+1)/2')
      END IF

      NY =JM+1
      NYW=JM+1
      NYS=JM+1
      NYC=JM+1
      NX =JM+1
      NXW=JM+1
      NXS=JM+1
      NXC=JM+1
      NZ =JM*(JM+1)*4
      NR =(MM+1)*(MM+1)
      NPY=2*(JM+1)*(MM+1)
      NPX=2*(JM+1)*(MM+1)
      NPM=(MM+2)*(JM+1)
      NSD=MM+2
      NH =2*IM
      NWFFT=4*IM+15

      NWY=2*(JM+1)
      NWX=2*(JM+1)
      NPMJ=MM+2

      IY=1
      IYW=IY +NY
      IYS=IYW+NYW
      IYC=IYS+NYS
      IX =IYC+NYC
      IXW=IX +NX
      IXS=IXW+NXW
      IXC=IXS+NXS
      IZ =IXC+NXC
      IR =IZ +NZ
      IPY=IR +NR
      IPX=IPY+NPY
      IPM=IPX+NPX
      ISD=IPM+NPM
      IH =ISD+NSD
      IWFFT=IH+NH

      IWY=IY
      IWX=IX
      IPMJ=IPM

      CALL SHINIT(MM,JM,IM,Q(IY),Q(IYS),Q(IYC),Q(IX),Q(IXS),Q(IXC),
     &            Q(IXW),Q(IZ),Q(IPY),Q(IPX),Q(IR),Q(IWFFT))

      RETURN
*----------------------------------------------------------------------
*     OPERATE LAPLACIAN
*----------------------------------------------------------------------
      ENTRY SHTLAP(MM,IND,A,B)

      CALL SHOLAP(MM,IND,A,B)

      RETURN
*----------------------------------------------------------------------
*     CALCULATE THE POSITION OF A SPECTRUM COEFFICIENT OF P_N^M
*----------------------------------------------------------------------
      ENTRY SHTNML(MM,N,M,LR,LI)

      CALL SHNM2L(MM,N,M,LR,LI)

      RETURN
*----------------------------------------------------------------------
*     CALCULATE LEGENDRE FUNCTIONS
*----------------------------------------------------------------------
      ENTRY SHTFUN(MM,JM,M,FUN,Q)

      CALL SHPFUN(MM,JM,M,FUN,Q(IPM),Q(IYS),Q(IYC),Q(IPY),Q(IR))

      RETURN
*----------------------------------------------------------------------
*     FORWARD LEGENDRE TRANSFORMATION
*----------------------------------------------------------------------
      ENTRY SHTLFW(MM,JM,M,ISW,WM,SM,Q)

      CALL SHLFWU(MM,JM,M,ISW,WM,SM,
     &    Q(ISD),Q(IXS),Q(IXC),Q(IPM),Q(IWY),Q(IWX),Q(IPX),Q(IR),Q(IZ))

      RETURN
*----------------------------------------------------------------------
*     BACKWARD LEGENDRE TRANSFORMATION
*----------------------------------------------------------------------
      ENTRY SHTLBW(MM,JM,M,ISW,SM,WM,Q)

      CALL SHLBWU(MM,JM,M,ISW,SM,WM,
     &            Q(ISD),Q(IPM),Q(IYS),Q(IYC),Q(IPY),Q(IR))

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM SPECTRA INTO WAVE
*----------------------------------------------------------------------
      ENTRY SHTS2W(MM,JM,ISW,S,W,Q)

      CALL SHMS2W(MM,JM,ISW,S,W,
     &            Q(ISD),Q(IPM),Q(IYS),Q(IYC),Q(IPY),Q(IR))

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM SPECTRA INTO WAVE (FROM M=M1 TO M=M2)
*----------------------------------------------------------------------
      ENTRY SHTSWA(MM,JM,ISW,M1,M2,S,W,Q)

      CALL SHMSWA(MM,JM,ISW,M1,M2,S,W,
     &            Q(ISD),Q(IPM),Q(IYS),Q(IYC),Q(IPY),Q(IR))

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM SPECTRA INTO WAVE (FOR M=0)
*----------------------------------------------------------------------
      ENTRY SHTSWZ(MM,JM,ISW,S,WZ,Q)

      CALL SHMSWZ(MM,JM,ISW,S,WZ,
     &            Q(ISD),Q(IPM),Q(IYS),Q(IYC),Q(IPY),Q(IR))

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM SPECTRA INTO WAVE (FOR M>0)
*----------------------------------------------------------------------
      ENTRY SHTSWM(MM,JM,M,ISW,S,WR,WI,Q)

      CALL SHMSWM(MM,JM,M,ISW,S,WR,WI,
     &            Q(ISD),Q(IPM),Q(IYS),Q(IYC),Q(IPY),Q(IR))

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM SPECTRA INTO WAVE FROM M=M1 TO M=M2 AT A LATITUDE
*----------------------------------------------------------------------
      ENTRY SHTSWJ(MM,JM,ISW,J,M1,M2,S,WJ,Q)

      CALL SHMSWJ(MM,JM,ISW,J,M1,M2,S,WJ,
     &            Q(ISD),Q(IPMJ),Q(IYS),Q(IYC),Q(IPY),Q(IR))

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM WAVE INTO SPECTRA
*----------------------------------------------------------------------
      ENTRY SHTW2S(MM,JM,ISW,S,W,Q)

      CALL SHMW2S(MM,JM,ISW,S,W,
     &    Q(ISD),Q(IPM),Q(IXS),Q(IXC),Q(IPX),Q(IWY),Q(IWX),Q(IR),Q(IZ))

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM WAVE TO GRID
*----------------------------------------------------------------------
      ENTRY SHTW2G(MM,JM,IM,W,G,Q)

      CALL SHFW2G(MM,JM,IM,W,G,Q(IH),Q(IWFFT))

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM WAVE TO GRID FROM M=M1 TO M=M2
*----------------------------------------------------------------------
      ENTRY SHTWGA(MM,JM,IM,M1,M2,W,G,Q)

      CALL SHFWGA(MM,JM,IM,M1,M2,W,G,Q(IH),Q(IWFFT))

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM WAVE TO GRID FOR M>0
*----------------------------------------------------------------------
      ENTRY SHTWGM(MM,JM,IM,M,WR,WI,G,Q)

      CALL SHFWGM(JM,IM,M,WR,WI,G,Q(IH),Q(IWFFT))

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM WAVE TO GRID FOR M=0
*----------------------------------------------------------------------
      ENTRY SHTWGZ(JM,IM,WZ,G)

      CALL SHFWGZ(JM,IM,WZ,G)

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM WAVE TO GRID FROM M=M1 TO M=M2 AT A LATITUDE
*----------------------------------------------------------------------
      ENTRY SHTWGJ(MM,IM,M1,M2,WJ,GJ,Q)

      CALL SHFWGJ(MM,IM,M1,M2,WJ,GJ,Q(IH),Q(IWFFT))

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM GRID INTO WAVE
*----------------------------------------------------------------------
      ENTRY SHTG2W(MM,JM,IM,G,W,Q)

      CALL SHFG2W(MM,JM,IM,G,W,Q(IH),Q(IWFFT))

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM SPECTRA INTO GRID
*----------------------------------------------------------------------
      ENTRY SHTS2G(MM,JM,IM,ISW,S,W,G,Q)

      CALL SHMS2W(MM,JM,ISW,S,W,
     &            Q(ISD),Q(IPM),Q(IYS),Q(IYC),Q(IPY),Q(IR))
      CALL SHFW2G(MM,JM,IM,W,G,Q(IH),Q(IWFFT))

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM SPECTRA INTO GRID (FROM M=M1 TO M=M2)
*----------------------------------------------------------------------
      ENTRY SHTSGA(MM,JM,IM,ISW,M1,M2,S,W,G,Q)

      CALL SHMSWA(MM,JM,ISW,M1,M2,S,W,
     &            Q(ISD),Q(IPM),Q(IYS),Q(IYC),Q(IPY),Q(IR))
      CALL SHFWGA(MM,JM,IM,M1,M2,W,G,Q(IH),Q(IWFFT))

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM SPECTRA INTO GRID (FOR M=0)
*----------------------------------------------------------------------
      ENTRY SHTSGZ(MM,JM,IM,ISW,S,WZ,G,Q)

      CALL SHMSWZ(MM,JM,ISW,S,WZ,
     &            Q(ISD),Q(IPM),Q(IYS),Q(IYC),Q(IPY),Q(IR))
      CALL SHFWGZ(JM,IM,WZ,G)

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM SPECTRA INTO GRID (FOR M>0)
*----------------------------------------------------------------------
      ENTRY SHTSGM(MM,JM,IM,M,ISW,S,WR,WI,G,Q)

      CALL SHMSWM(MM,JM,M,ISW,S,WR,WI,
     &            Q(ISD),Q(IPM),Q(IYS),Q(IYC),Q(IPY),Q(IR))
      CALL SHFWGM(JM,IM,M,WR,WI,G,Q(IH),Q(IWFFT))

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM SPECTRA INTO GRID FROM M=M1 TO M=M2 AT A LATITUDE
*----------------------------------------------------------------------
      ENTRY SHTSGJ(MM,JM,IM,ISW,J,M1,M2,S,WJ,GJ,Q)

      CALL SHMSWJ(MM,JM,ISW,J,M1,M2,S,WJ,
     &            Q(ISD),Q(IPMJ),Q(IYS),Q(IYC),Q(IPY),Q(IR))
      CALL SHFWGJ(MM,IM,M1,M2,WJ,GJ,Q(IH),Q(IWFFT))

      RETURN
*----------------------------------------------------------------------
*     TRANSFORM GRID INTO SPECTRA
*----------------------------------------------------------------------
      ENTRY SHTG2S(MM,JM,IM,ISW,G,W,S,Q)

      CALL SHFG2W(MM,JM,IM,G,W,Q(IH),Q(IWFFT))
      CALL SHMW2S(MM,JM,ISW,W,S,
     &    Q(ISD),Q(IPM),Q(IXS),Q(IXC),Q(IPX),Q(IWY),Q(IWX),Q(IR),Q(IZ))

      END
