    module gfs_init
    contains
!
!     ******************************************************************
!     ******************************************************************
!     ******                                                      ******
!     ******                SUBROUTINE INIT_VERT                  ******
!     ******                                                      ******
!     ******************************************************************
!     ******************************************************************
!
      SUBROUTINE init_vert(cl,sl,del,si,tov,rmu,rnu,rdel2, &
                           rk,slk,sl1,levshc)

      use resolution
      use constant
      use precision,only:D
      implicit none

      real(D) cl(levs),sl(levs),del(levs),si(levp1),tov(levs),rho(levs), &
           rmu(levs),rnu(levs),rdel2(levs),rk,slk(levs),sl1
      integer levshc

      integer k
      real(D) sigshc

      DO k=1,LEVS
        CL(k) = 1.E0 - SL(k)
      ENDDO

      DO k=1,LEVS
        DEL(k)= SI(k)-SI(k+1)
      ENDDO
      DO k=1,LEVS
        TOV(k)=300.E0
      ENDDO

      DO k=1,LEVM1
        RHO(K)=LOG(SI(K)/SI(K+1))
      ENDDO
      RHO(LEVS)=0.

      DO k=1,LEVS
        RMU(K)=1.-SI(K+1)*RHO(K)/DEL(K)
      ENDDO

      DO k=1,LEVM1
        RNU(K+1)=-1.+SI(K)*RHO(K)/DEL(K)
      ENDDO
      RNU(1)=0.

      DO k=1,LEVS
        RDEL2(k)=1./(2.*DEL(k))
      ENDDO

      RK = RD/CP
      sl1=SL(1)
      levshc=0
      sigshc=0.7
      DO k=1,LEVS
        SLK(K)=SL(K)**RK
        IF(SL(K).GT.SIGSHC) LEVSHC=K+1
      ENDDO

      RETURN
      END SUBROUTINE init_vert
!
!***********************************************************************
!
      SUBROUTINE GLATS(COLRAD,WGT,RCS2,SINLAT,COSLAT,colat1)
!
!***********************************************************************
!
      use resolution
      use precision,only:D
      implicit none
      real(D) COLRAG(LATG2),WGTG(LATG2),RCS2G(LATG2), &
           COLRAD(latg2),WGT (latg2),RCS2 (latg2), &
           SINLAT(latg2),COSLAT(latg2),colat1
      real(D) EPS,SI,RL2,SCALE,PI,DRADZ,RAD,DRAD,X,W,SN,RC,P2,P1
      integer L2,K1,K,ITER,J
      include 'machine.h'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! EPS IS MACHINE DEPENDENT (WORDSIZE IN BITS)
      EPS=PRECIS
!     PRINT 101
!101  FORMAT ('0 I   COLAT   COLRAD     WGT', 12X,'ITER  RES')
      SI = 1.0
      L2=2*LATG2
      RL2=L2
      SCALE = 2.0/(RL2*RL2)
      K1=L2-1
      PI = ATAN(SI)*4.E+00
      DRADZ = PI / 360.
      RAD = 0.0
      DO 1000 K=1,LATG2
!        print *, k 
        ITER=0
        DRAD=DRADZ
1       CALL POLY(L2,RAD,P2)
2       P1 =P2
        ITER=ITER+1
        RAD=RAD+DRAD
        CALL POLY(L2,RAD,P2)
!       print *, SIGN(SI,P1), SIGN(SI,P2) 
       IF(SIGN(SI,P1).EQ.SIGN(SI,P2)) GO TO 2
!       print *, DRAD, EPS
       IF(DRAD.LT.EPS)GO TO 3
        RAD=RAD-DRAD
        DRAD = DRAD * 0.25
       GO TO 1
3      CONTINUE
        COLRAG(K)=RAD
        CALL POLY(K1,RAD,P1)
        X = COS(RAD)
        W = SCALE * (1.0 - X*X)/ (P1*P1)
        WGTG(K) = W
        SN = SIN(RAD)
        W=W/(SN*SN)
        RC=1./(SN*SN)
        RCS2G(K) = RC
        CALL POLY(L2,RAD,P1)
!       PRINT 102,K,PHI,COLRAG(K),WGTG(K),ITER,P1
!102    FORMAT(1H ,I2,2X,F6.2,2X,F10.7,2X,E13.7,2X,I4,2X,D13.7)
1000  CONTINUE
!mg      PRINT 100,LATG2
!100   FORMAT(1H ,'SHALOM FROM 0.0E0 GLATS FOR ',I3)
!yt
      colat1=colrag(1)
      DO j=1,latg2
        COLRAD(j) = COLRAG(j)
        WGT(j)    = WGTG  (j)
        RCS2(j)   = RCS2G (j)
        SINLAT(j) = COS(COLRAD(j))
        COSLAT(j) = SQRT(1.0 - SINLAT(j)*SINLAT(j))
      ENDDO

      RETURN
      END SUBROUTINE GLATS
!
!***********************************************************************
!
      SUBROUTINE POLY(N,RAD,P)
!
!***********************************************************************
!
      use precision,only:D
      implicit none
      real(D) X,Y1,Y2,G,Y3,P,RAD
      integer I,N
      X = COS(RAD)
      Y1 = 1.0
      Y2=X
      DO 1 I=2,N
      G=X*Y2
      Y3=G-Y1+G-(G-Y1)/FLOAT(I)
      Y1=Y2
      Y2=Y3
1     CONTINUE
      P=Y3
      RETURN
      END SUBROUTINE POLY
!
!***********************************************************************
!
      SUBROUTINE PLN2I(PLNT,PDERT,PLNTOP,RCS2,COLRAD)
!
!***********************************************************************
!
      use resolution
      use constant
      use precision,only:D
      implicit none
      real(D) PLNT(LNTD)
      real(D) PLNV(LNUT)
      real(D) COLRAD
      real(D) PLNTOP(0:JCAP)
!
!SELA ADDED FROM GOZRIN
      real(D) PDERT(LNTD)
      real(D) RCS2
!SELA ADDED FROM GOZRIN
!
      real(D) W(0:JCAP),EPSI,REPSI,X,Y
!
      COMMON /PLN2TI/ EPSI(LNUT),REPSI(LNUT), &
                      X(0:JCAP),Y(0:JCAP)
!
      real(D) XA,XB
      COMMON /GOZCOM/ XA(LNTD),XB(LNTD)
!
      INTEGER OFFSRI
      INTEGER OFFVRI
      INTEGER N,N0,IND,I,L
      REAL(D) SINLA,COS2,PROD,WCSA
!
!     OFFSRI IS THE SCALAR REAL OR IMAGINARY TRIANGLE INDEX IN WORDS.
      OFFSRI(N,L) = JCAP*L + N + 1 + L*(1-L)/2
!
!     OFFVRI IS THE VECTOR REAL OR IMAGINARY TRIANGLE INDEX IN WORDS.
      OFFVRI(N,L) = JCAP*L + N + 1 + L*(3-L)/2
!
!     DATA IFIR /0/
!           PART BETWEEN GUARDS MADE INTO SR GPLN2I.
!           7 DEC 1990      M. ROZWODOSKI
!
      SINLA  = COS(COLRAD)
      COS2   = 1.0 - SINLA * SINLA
      PROD   = 1.0
      DO 600 L=0,JCAP
             W(L) = 0.5*PROD
!            IF (PROD .LT. 1.0E-75)  PROD=0.0
             PROD = PROD*COS2*Y(L)
  600 CONTINUE
      DO 620 L=0,JCAP
             PLNV(OFFVRI(L,L)) = SQRT(W(L))
  620 CONTINUE
!
      DO 700 L=0,JCAP
             PLNV(OFFVRI(L+1,L)) = PLNV(OFFVRI(L,L)) * X(L) * SINLA
  700 CONTINUE
!
      DO 740 N0=2,JCAP+1
      DO 720 L=0,JCAP+1-N0
             PLNV(OFFVRI(N0+L  ,L)) =                            &
          (  PLNV(OFFVRI(N0+L-1,L)) * SINLA                      &
           - PLNV(OFFVRI(N0+L-2,L)) *  EPSI(OFFVRI(N0+L-1,L))  ) &
                                    * REPSI(OFFVRI(N0+L  ,L))
  720 CONTINUE
  740 CONTINUE
!
      DO 860 L=0,JCAP
      DO 840 N=L,JCAP
             PLNT(OFFSRI(N,L)) = PLNV(OFFVRI(N,L))
  840 CONTINUE
  860 CONTINUE
!yt   Keep values for vectors
      DO l=0,JCAP
        PLNTOP(l)=PLNV(OFFVRI(JCAP1,l))
      ENDDO
!
!SELA ADDED FROM GOZRIN
!............................................................
!     COMPUTE PLN DERIVATIVES.
      WCSA=RCS2/REARTH
!
      ind=0
      DO l=0,JCAP
        ind=ind+1
        PDERT(ind)=0.
        DO n=l+1,JCAP
          ind=ind+1
          PDERT(ind)=XA(ind)*PLNT(ind-1)
        ENDDO
      ENDDO
!
      ind=0
      DO l=0,JCAP
        DO n=l,JCAP-1
          ind=ind+1
          PDERT(ind)=PDERT(ind)+XB(ind)*PLNT(ind+1)
        ENDDO
        ind=ind+1
        PDERT(ind)=PDERT(ind)+XB(ind)*PLNTOP(l)
      ENDDO
!
      DO i=1,LNT
        PDERT(i)=PDERT(i)*wcsa
      ENDDO
!
!............................................................
!SELA ADDED FROM GOZRIN
      RETURN
      END SUBROUTINE PLN2I
!
!***********************************************************************
!
      SUBROUTINE GPLN2I
!
!***********************************************************************
!
      use resolution
      use precision,only:D
      implicit none
!C
      real(D) EPSI,REPSI,X,Y,TEMP
      integer N,L
      COMMON /PLN2TI/ EPSI(LNUT),REPSI(LNUT), &
                      X(0:JCAP),Y(0:JCAP)
!CC
!C    OFFVRI IS THE VECTOR REAL OR IMAGINARY TRIANGLE INDEX IN WORDS.
      INTEGER OFFVRI
      OFFVRI(N,L) = JCAP*L + N + 1 + L*(3-L)/2
!C
      DO 240 L=0,JCAP
              EPSI(OFFVRI(L,L))=0.0
             REPSI(OFFVRI(L,L))=0.0
             DO 220 N=L+1,JCAP+1
                    TEMP=((N*N-L*L)/(4.*N*N-1.))
                     EPSI(OFFVRI(N,L))=SQRT(    TEMP)
                    REPSI(OFFVRI(N,L))=SQRT(1.0/TEMP)
  220        CONTINUE
  240 CONTINUE
!C
      DO 400 L=0,JCAP
             X(L) = (L+1)*2+1
  400 CONTINUE
      DO 420 L=0,JCAP
             Y(L) = X(L)/(X(L)-1.)
  420 CONTINUE
      DO 440 L=0,JCAP
             X(L) = SQRT(X(L))
  440 CONTINUE
!C
      RETURN
      END SUBROUTINE GPLN2I
!
!***********************************************************************
!
      SUBROUTINE GGOZRM
!
!***********************************************************************
!
      use resolution
      use precision,only:D
      implicit none
      real(D) epsi(LNT),epstop(JCAP1),XA(LNTD),XB(LNTD)
      integer ind,n,l
!
      COMMON /GOZCOM/XA,XB
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@c
!
      ind=0
      DO l=0,JCAP
        ind=ind+1
        epsi(ind)=0.
        DO n=l+1,JCAP
          ind=ind+1
          epsi(ind)=sqrt( (n*n-l*l)/(4.*n*n-1.) )
        ENDDO
        n=JCAP1
        epstop(l+1)=sqrt( (n*n-l*l)/(4.*n*n-1.) )
      ENDDO
!
      ind=0
      DO l=0,JCAP
        DO n=l,JCAP
          ind=ind+1
          xa(ind)=(n+1.)*epsi(ind)
        ENDDO
      ENDDO
!
      ind=0
      DO l=0,JCAP
        DO n=l,JCAP-1
          ind=ind+1
          xb(ind)=-n*epsi(ind+1)
        ENDDO
        ind=ind+1
        xb(ind)=-JCAP*epstop(l+1)
      ENDDO
!
      RETURN
      END SUBROUTINE GGOZRM
!
!***********************************************************************
!
      SUBROUTINE ToEvenOdd(num,fln)
!
!***********************************************************************
!
      use resolution
      use precision,only:D
      implicit none
      integer num,i,l,k,np,ni
      real(D) fln(lntd,num),buff(JCAP1)

      DO k=1,num
        i=1
        DO l=0,JCAP-2
          np=(JCAP1-l+1)/2
          ni=(JCAP1-l)/2
          CALL dcopy(JCAP1-l,fln(i,k),1,buff,1)
          CALL dcopy(np,buff(1),2,fln(i,k),1)
          CALL dcopy(ni,buff(2),2,fln(i+np,k),1)
          i=i+JCAP1-l
        ENDDO
      ENDDO

      RETURN
      END SUBROUTINE ToEvenOdd
!
!***********************************************************************
!
      SUBROUTINE FromEvenOdd(num,fln)
!
!***********************************************************************
!
      use resolution
      use precision,only:D
      implicit none
      integer num,i,l,k,np,ni
      real(D) fln(lntd,num),buff(JCAP1)

      DO k=1,num
        i=1
        DO l=0,JCAP-2
          np=(JCAP1-l+1)/2
          ni=(JCAP1-l)/2
          CALL dcopy(JCAP1-l,fln(i,k),1,buff,1)
          CALL dcopy(np,buff      ,1,fln(i,k),2)
          CALL dcopy(ni,buff(np+1),1,fln(i+1,k),2)
          i=i+JCAP1-l
        ENDDO
      ENDDO

      RETURN
      END SUBROUTINE FromEvenOdd
!
      end module gfs_init
