      SUBROUTINE CALC_GPH(GP,THROW,SHTJ,PRESSG,PHIS,RGAS,ILG,ILEV,
     1                    LEV,IL1,IL2,ILEVP1)
C
C     * JASON COLE - APRIL 10, 2019  UPDATED TO WORK ON XC40, REMOVE IBM
C     *                              INTRINSICS AND USE GENERIC ROUTINES.
C     * JASON COLE - MARCH 30, 2010. NEW ROUTINE FOR GCM15H FOR
C     *                              CALCULATION OF GEOPOTENTIAL HEIGHT
C     *                              AT MODEL LEVELS.
C     *                              BASED ON PAHGT
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
C
      REAL, INTENT(OUT), DIMENSION(ILG,ILEV)   :: GP
      REAL, INTENT(IN),  DIMENSION(ILG,ILEVP1) :: THROW
      REAL, INTENT(IN),  DIMENSION(ILG,LEV)    :: SHTJ
      REAL, INTENT(IN),  DIMENSION(ILG)        :: PRESSG,PHIS
      REAL, ALLOCATABLE, DIMENSION(:,:)        :: TERM1,TERM2,PF
C
C---------------------------------------------------------------------
C     * CALCULATE GEOPOTENTIAL HEIGHT
C
      ALLOCATE(PF(ILG,ILEV))
      DO L=1,ILEV
        PF(IL1:IL2,L)=SHTJ(IL1:IL2,L+1)*PRESSG(IL1:IL2)
      ENDDO
      ILEVM=ILEV-1
      GP(IL1:IL2,ILEV) = 0.
      ILT=IL2-IL1+1
      ALLOCATE(TERM1(ILT,ILEVM))
      ALLOCATE(TERM2(ILT,ILEVM))
      DO L=1,ILEVM
        TERM1(1:ILT,L)=PF(IL1:IL2,L+1)/PF(IL1:IL2,L)
      ENDDO
      DO L=1,ILEVM
        DO I=1,ILT
          TERM2(I,L) = LOG(TERM1(I,L))
        ENDDO
      ENDDO

      DO L=ILEVM,1,-1
        GP(IL1:IL2,L)= GP(IL1:IL2,L+1)
     1               + RGAS*THROW(IL1:IL2,L+2)*TERM2(1:ILT,L)
      ENDDO

! Add surface geopotential
      DO L=1,ILEV
        DO IL = IL1, IL2
         GP(IL,L)=GP(IL,L)+PHIS(IL)
       END DO ! IL
      END DO ! L

      DEALLOCATE (TERM1)
      DEALLOCATE (TERM2)
      DEALLOCATE (PF)
C
      END SUBROUTINE CALC_GPH
