      SUBROUTINE GEOTHICK(dz,kdim,temp,maxlevel,maxcalc,ik,rlat,klev)

C     GEOTHICK CALCULATES THE GEOPOTENTIAL THICKNESS
C     OF A GIVEN PROFILE BASED ON ITS DENSITY,
C     DEPTH, AND POSITION

C     GEOPOTENTIAL THICKNESS WILL BE IN UNITS OF DYNAMIC CM

c     dz - depths
c     temp - array holding input density information/
c            output geopotential thickness information
c     kdim - maximum number of depths
c     maxcalc - maximum number of parameters in temp
c     ik - parameter number
c     rlat - latitude of profile
c     klev - number of depth levels

      parameter (kdimmax=6000,bmiss=-1.E10)

      dimension dz(kdim),temp(maxlevel,maxcalc)
      dimension press(kdimmax)
      if ( dz(1) .eq. 0. ) then
       press(1) =0.
      else
       press(1)=fpress(dz(1),rlat,0.,0.,3)
      endif

c     calculate pressure based on depth and position

      kl=klev
      if ( klev .ge. kdim ) kl=kdim-1

      do 200 k=2,kl+1
200     press(k)=fpress(dz(k),rlat,dz(k-1),press(k-1),3)
      write(6,*) 'press(k)',press(k)
 
      do 20 k=1,kl

       if ( temp(k,ik) .gt. bmiss ) then

c     calculate layer thickness

        if ( dz(k)  .gt. .0001 ) then
         x0=press(k-1)+((press(k)-press(k-1))/2.)
        else
         x0=0.
        endif
        xz=press(k)+((press(k+1)-press(k))/2.)

        xthick=xz-x0
         
c     calculate standard geopotential

        astd= 1./(fdense(0.,35.,press(k))+1000.)

c      calculate geopotential

        a= 1./(temp(k,ik)+1000.)

c      calculate geopotential thickness

c       write(6,*) 'Level',k,temp(k,ik),fdense(0.,35.,press(k)),
c    * a, astd
        temp(k,ik)=(a-astd)*xthick*10.E4
c       write(6,*) press(k-1),press(k),press(k+1),x0,xz,xthick
c       write(6,*) temp(k,ik)

       endif

20    continue
   
      return
      end
