      subroutine heights(dset,xlat,xlon,hireskey,convkey,
     &                   zradius,zareoid,zsurface,ier)

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Purpose:
c =======
c heights is a subroutine that can be used to convert vertical coordinates
c expressed as distance to the center of the planet, height above the
c areoid (zero datum) and height above the local surface. Given any of these,
c this routine computes the other two, according to the "convkey" flag.
c
c This routine works at two distinct resolution, according to flag "hireskey".
c At low resolution, distances are computed with respect to GCM (64x48) grid.
c At high resolution, MOLA 32 pix/deg topography (and high resolution areoid,
c computed from mgm1025 spherical harmonics expansion) are used.
c
c The main routine call_mcd does such conversions; these are provided here
c for users interested in a light and fast tool to do just that.
c
c "Low resolution" procedure based on MCD v4.1 "calc_height" routine
c
c Returned status/error code 'ier' are identical to call_mcd ones:
c       ier=0:  No problem
c       ier=4:  Wrong value fpr input flag 'hireskey'
c       ier=7:  Wrong value for input latitude 'xlat'
c       ier=15: Could not open a datafile
c       ier=16: Failed loading data from a database file
c
c       ier=24: Wrong value for input switch 'convkey'
c
c Additionnaly, if an error occurs, zradius,zareoid and zsurface are
c are set to -999.0
c
c Author:
c ======
c        E. Millour 12/2006
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

       implicit none
       include "netcdf.inc"
       include "constants_mcd.inc"

c Arguments:
c    inputs:
       character*(*) dset ! Path to MCD datafiles
                          ! if an empty string (ie: dset='')
                          ! then default dset='MCD_DATA/' path is assumed
       real xlon          ! East longitude (in degrees)
       real xlat          ! North latitude (in degrees)
       integer hireskey   ! resolution flag (=0 GCM resolution;=1 high res)
       integer convkey    ! switch to choose which z variable is known
                          ! and used to find the other two:
               ! convkey=1 zradius is known, compute zareoid and zsurface
               ! convkey=2 zareoid is known, compute zradius and zsurface
               ! convkey=3 zsurface is known, compute zradius and zareoid
c inputs/outputs:
       real zradius  ! distance to center of planet (m)
       real zareoid  ! altitude above areoid (m)
       real zsurface ! altitude above local surface (m)
c   outputs:
       integer ier ! status/error code (=0 if OK)

c local variables:
      integer lendataset ! length of 'dset' string
      character*140 path ! path given by 'dset'
      character*140 datafile ! full path to GCM 'mountain.nc' data file
      integer unet ! 'mounatin.nc' file ID
      character*16 name ! variable name
      integer varid     ! variable ID
      integer ierr1,ierr2 ! to store NetCDF routine error codes

      real lon  ! local value of longitude (degrees east)
      real lat  ! local value of latitude (degrees north)
      
      integer dlon(4)  !index, along longitudes, of database points
      ! dummy initialization to get rid of compiler warnings
      data dlon /1,1,1,1/
      integer dlat(4)  !index, along latitudes, of database points
      ! dummy initialization to get rid of compiler warnings
      data dlat /1,1,1,1/
      real    t        !weight (normalized longitudinal distance to point 1)
      real    u        !weight (normalized latitudinal distance to point 1)
      data u /0/   ! dummy initialization to get rid of compiler warnings
      data t /0/   ! dummy initialization to get rid of compiler warnings
      real    alon,alat
      integer i

      real orog(dimlon,dimlat)   ! to store GCM orography
      real areoid(dimlon,dimlat) ! to store GCM areoid
      save orog, areoid          ! saved for further use
      logical firstcall 
      data firstcall /.true./
      save firstcall

      real xorog ! orography at position lon,lat
      real xareoid ! distance to center of planet of areoid at position lon,lat
      
      
      ier=0

c 1. Check input arguments
c 1.1 Check latitude
      if (abs(xlat).gt.90.0) then
        if (output_messages) then
         write(out,*)'HEIGHTS Error: wrong value for latitude=',xlat
        endif
        ier=7
        goto 999
      else
        lat=xlat
      endif

c 1.2 Set longitude
      ! At GCM resolution, longitude must be in [-180:180]
      lon=mod(xlon,360.0)
      if (lon.lt.-180.) then ! in case lon in [-360:-180]
        lon=lon+360.
      endif
      if (lon.gt.180.0) then ! in case lon in [180:360]
        lon=lon-360.
      endif
      ! Note: at high resolution, 'mola' subroutine
      ! needs longitude to be in [0:360] but handles the conversion

c 1.3 Check hireskey
      if ((hireskey.lt.0).or.(hireskey.gt.1)) then
        if (output_messages) then
         write(out,*)'HEIGHTS Error: wrong value for parameter hireskey'
         write(out,*)'               hireskey=',hireskey
        endif
        ier=4
        goto 999
      endif

c 1.4 Check convkey
      if ((convkey.lt.1).or.(convkey.gt.3)) then
        if (output_messages) then
          write(out,*)'HEIGHTS Error: wrong value for parameter convkey'
          write(out,*)'               convkey=',convkey
        endif
        ier=24
        goto 999
      endif

c 1.5 set 'dset' to default 'MCD_DATA' , if necessary
c     find last non-blank character and length of 'dset' srting
      do lendataset=len(dset),1,-1
         if (dset(lendataset:lendataset).ne.' ') goto 100
      enddo
  100 continue

      if (lendataset.eq.0) then ! 'dset' not specified; set default path
        path='MCD_DATA/'
        lendataset=9
      else
        path=dset
      endif

c 2. If first call and low resolution, load GCM orography and areoid
c    from 'mountain.nc' file
      if (firstcall .and.(hireskey.eq.0)) then
        datafile=path(1:lendataset)//'mountain.nc'
        ! open datafile
        ier=NF_OPEN(datafile,NF_NOWRITE,unet)
        if (ier.ne.0) then
          ier=15
          if (output_messages) then
            write(out,*)'HEIGHTS Error: cannot open file ',datafile
          endif
          goto 999
        endif
        
        ! load orography from file
        name='orography'
        ierr1=NF_INQ_VARID(unet,name,varid)
        ierr2=NF_GET_VAR_REAL(unet,varid,orog)
        if ((ierr1.ne.NF_NOERR).or.(ierr2.ne.NF_NOERR)) then
          ier=16
          if (output_messages) then
            write(out,*)'HEIGHTS Error: failed loading ',name
          endif
          goto 999
        endif
        
        ! load areoid from file
        name='areoid'
        ierr1=NF_INQ_VARID(unet,name,varid)
        ierr2=NF_GET_VAR_REAL(unet,varid,areoid)
        if ((ierr1.ne.NF_NOERR).or.(ierr2.ne.NF_NOERR)) then
          ier=16
          if (output_messages) then
            write(out,*)'HEIGHTS Error: failed loading ',name
          endif
          goto 999
        endif
        
        firstcall=.false.
      endif

c 3. get values of orography and areoid at position lon,lat
      if (hireskey.eq.0) then ! low resolution
      ! find longitude of four nearest points (as in grid4 routine)
        if (lon.ge.lonmax) then
         dlon(1)=dimlon
         dlon(2)=1
         dlon(3)=1
         dlon(4)=dimlon
         t=(lon-lonmax)/deltalon
        else
         alon=lonmin
         do i=1,dimlon-1
            if ((lon.ge.alon).and.(lon.lt.alon+deltalon)) then
               dlon(4)=i
               dlon(3)=i+1
               dlon(1)=i
               dlon(2)=i+1
               t=(lon-alon)/deltalon
            endif
            alon=alon+deltalon
         enddo
        endif
      ! find latitude of four nearest points (as in grid4)
        if (lat.ge.latmax) then
         dlat(4)=1
         dlat(3)=1
         dlat(1)=1
         dlat(2)=1
         u=0.
        elseif (lat.le.latmin) then
         dlat(4)=dimlat
         dlat(3)=dimlat
         dlat(1)=dimlat
         dlat(2)=dimlat
         u=0.
        else
         alat=latmin
         do i=dimlat,2,-1
            if ((lat.ge.alat).and.(lat.lt.alat+deltalat)) then
               dlat(1)=i
               dlat(2)=i
               dlat(3)=i-1
               dlat(4)=i-1
               u=(lat-alat)/deltalat
            endif
            alat=alat+deltalat
         enddo       
        endif
      
      ! bilinear interpolation to find local orography and areoid
        xorog=(1.-t)*(1.-u)*orog(dlon(1),dlat(1))
     &       +t*(1.-u)*orog(dlon(2),dlat(2))
     &       +t*u*orog(dlon(3),dlat(3))
     &       +(1.-t)*u*orog(dlon(4),dlat(4))
        xareoid=(1.-t)*(1.-u)*areoid(dlon(1),dlat(1))
     &       +t*(1.-u)*areoid(dlon(2),dlat(2))
     &       +t*u*areoid(dlon(3),dlat(3))
     &       +(1.-t)*u*areoid(dlon(4),dlat(4))

      else ! high resolution mode
        ! orography interpolated from MOLA file
        call mola(path(1:lendataset),lat,lon,xorog,ier)
        if (ier.ne.0) then
          goto 999 ! ier is set in 'mola' and error message is issued there
        endif
        
        ! areoid computed from mgm1025 spherical harminics expansion
        call molareoid(path(1:lendataset),lon,lat,xareoid)
        if (ier.ne.0) then
          goto 999 ! ier is set in 'molareoid' and error message is issued there
        endif
        
      endif ! of if (hireskey.eq.0)
      
c 4. convert between heights
      if (convkey.eq.1) then
         zareoid=zradius-xareoid
         zsurface=zareoid-xorog
      else if (convkey.eq.2) then
         zradius = zareoid+xareoid
         zsurface = zareoid-xorog
      else ! convkey.eq.3
         zareoid=zsurface+xorog
         zradius=zareoid+xareoid
      endif

      return

c     error handling
  999 zradius=-999.0
      zareoid=-999.0
      zsurface=-999.0
      return
      
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine mola(dset,latitude,longitude,alt,ierr)
!Give the MOLA altitude (alt), given lat and lon coordinates
!Using bilinear interpolation from 32 pixels/degree MOLA file
!
!  12/2016 swiched some internal computations to double precision EM.
!
      implicit none

      include "netcdf.inc"
      include "constants_mcd.inc"

! Arguments
! inputs
      character*(*) dset ! Path to MCD datafiles 
      real latitude  ! north latitude (degrees)
      real longitude ! east longitude (degrees)
! outputs
      real alt     ! above areoid altitude of surface (meters)
      integer ierr ! returned status code (==0 if OK)

! Local variables

      logical firstcall
      data firstcall/.true./
      save firstcall
      character*140 molafile ! MOLA datafile 
c      data molafile/'mola_32.2.nc'/
c      data molafile/'mola16.nc'/
c      data molafile/'mola32.nc'/
      double precision resol
c      parameter(resol=16) ! MOLA pixel/degree resolution
      parameter(resol=32)
c      real invresol
c      parameter(invresol=1./resol)
      integer jjm, iim   ! # of longitude and latitude MOLA data values
      parameter(jjm=180*resol, iim=2*jjm)
      integer*2 altmola(iim,jjm) ! MOLA altitude dataset
      save altmola

      integer*2 mintopo_check,maxtopo_check ! known min and max of MOLA dataset
c      parameter(mintopo_check=-8156,maxtopo_check=21191) ! mola_32.2.nc
c      parameter(mintopo_check=-8177,maxtopo_check=21191) ! mola16.nc
      parameter(mintopo_check=-8206,maxtopo_check=21191) ! mola32.nc
      double precision dlat, dlon ! , lontmp
      integer i,j ! ,count
      double precision topo(4) ! neighboring MOLA points (for bilinear interpolation)
      integer latsup,latinf,loninf,lonsup ! indexes of neighboring points
      integer*2 mintopo, maxtopo ! min and max of read dataset
      integer nid,nvarid ! NetCDF file and variable IDs
      double precision colat ! colatitude
      double precision lat,lon ! longitude and latitude, local values (in degrees)

C 1. Load MOLA dataset upon first call
C    (dataset is stored in a 'save' array)

      if (firstcall) then
         firstcall=.false.

C 1.1. Open MOLA file
         molafile=dset//'mola32.nc'
         if (output_messages) then
           write(out,*)"Loading MOLA topography from file ",
     &                  trim(molafile)
         endif
         ierr = NF_OPEN (molafile, NF_NOWRITE,nid)
         if (ierr.NE.NF_NOERR) then
            if (output_messages) then
              write(out,*)"Error in mola: Could not open file ",
     &                     trim(molafile)
            endif
            ierr=15 !set appropriate error code
            return
         endif

C 1.2. Load data
         ierr = NF_INQ_VARID (nid, "alt", nvarid)
         ! note that MOLA "alt" are given as "short" (16 bits integers)
         ierr = NF_GET_VAR_INT2(nid, nvarid, altmola)
         if (ierr.ne.NF_NOERR) then
           if (output_messages) then
            write(out,*)"Error in mola: <alt> not found"
           endif
           ierr=16 ! set appropriate error code
           return
         endif


c 1.3. Close MOLA file
         ierr=NF_CLOSE(nid)

c 1.4 Check that the MOLA dataset was correctly loaded

         mintopo=mintopo_check
         maxtopo=maxtopo_check
         do i=1,iim
            do j=1,jjm
               mintopo=min(mintopo,altmola(i,j))
               maxtopo=max(maxtopo,altmola(i,j))
            enddo
         enddo
         if ((mintopo.ne.mintopo_check).or.
     &       (maxtopo.ne.maxtopo_check)) then
           if (output_messages) then
            write(out,*)"***ERROR Mola file ",molafile,
     &                  " is not well read"
            write(out,*) "Minimum found: ", mintopo
            write(out,*) "Minimum should be:",mintopo_check
            write(out,*) "Maximum found: ", maxtopo
            write(out,*) "Maximum should be:",maxtopo_check
           endif
           ierr=16
           return
         endif
         if (output_messages) then
           write(out,*) "Done reading MOLA data"
         endif
      endif ! End of if(firstcall)

C 2. Check that input longitude and latitude make sense
      lat=latitude
      if((lat.gt.90).or.(lat.lt.-90)) then
        if (output_messages) then
          write(out,*)"Error in mola: Wrong value for latitude"
        endif
        stop
      endif

c longitude must range from 0 to 360
      lon=longitude
      do while(lon.gt.360.)
        lon=lon-360.
      enddo
      
      do while(lon.lt.0.)
        lon=lon+360.
      enddo
      
C 3. Identify the four neighboring points from MOLA dataset
!    These points are arranged as follows:  1 2
!                                           3 4

      colat=90-lat

      if (colat.lt.1./(2.*resol)) then
         latsup=1 
         latinf=1
         dlat=0
      else if (colat.gt.180.-1./(2.*resol)) then
         latsup=jjm
         latinf=jjm
         dlat=0
      else
         latsup=1+int((colat-1./(2.*resol))*resol)
         latinf=latsup+1
         dlat=1-(colat-(1./(2.*resol)+(latsup-1)/resol))*resol
      endif
      ! handle cases where dlat is slightly out of range due to roundoffs
      if ((dlat.lt.0).and.(dlat.ge.-1.e-5)) dlat=0
      if ((dlat.gt.1).and.(dlat.le.1.00001)) dlat=1.      
! Note: dlat is the (normalized) "latitudinal distance" to point 3 
!       ie: dlat=0 if lat=latitude of point 3
!           dlat=1 if lat=latitude of point 1

      if ((lon.lt.1./(2.*resol)).or.(lon.ge.(360-1./(2.*resol)))) then
         loninf=iim
         lonsup=1
         if (lon.lt.1./(2.*resol)) then
            dlon=lon*resol+0.5
         else
            dlon=(lon-(1./(2.*resol)+(loninf-1)/resol))*resol
         endif
      else
         if (((lon-1./(2.*resol))*resol).ge.0) then
            loninf=1+int((lon-1./(2.*resol))*resol)
         else
            loninf=int((lon-1./(2.*resol))*resol)
         endif
         lonsup=loninf+1
         dlon=(lon-(1./(2.*resol)+(loninf-1)/resol))*resol
      endif
      ! handle cases where dlon is slightly out of range due to roundoffs
      if ((dlon.lt.0).and.(dlon.ge.-1.e-5)) dlon=0
      if ((dlon.gt.1).and.(dlon.le.1.00001)) dlon=1.
! Note: dlon is the (normalized) "longitudinal distance" to point 3
!       ie: dlon=0 if lon=longitude of point 3
!           dlon=1 if lon=longitude of point 4

      if((dlat.gt.1).or.(dlon.gt.1).or.(dlat.lt.0).or.(dlon.lt.0)) then
        if (output_messages) then
         write(out,*)"Error in mola: on dlat or dlon" 
         write(out,*) "dlat: ",dlat
         write(out,*) "lat: ",lat
         write(out,*) "dlon: ",dlon
         write(out,*) "lon: ",lon
        endif
        ierr=-5
        return
      endif

C 4. Interpolation

! Four nearest points are aranged as follows:  1 2
!                                              3 4
      topo(1)=dble(altmola(loninf,latsup))
      topo(2)=dble(altmola(lonsup,latsup))
      topo(3)=dble(altmola(loninf,latinf))
      topo(4)=dble(altmola(lonsup,latinf))
   
! Use bilinear interpolation to evaluate alt

      alt=(1.-dlon)*(1.-dlat)*topo(3)+(1.-dlat)*dlon*topo(4)
     &     +dlat*(1.-dlon)*topo(1)+dlat*dlon*topo(2)

      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

        subroutine molareoid(dset,lon,lat,rareoid)
C This subroutine returns the radial position (ie: distance to the center
C of Mars) of the reference areoid for a given position (given values of
C longitude and latitude)
C Based on "areoid.f" program by G. Neumann, which is available at
C ftp://ltpftp.gsfc.nasa.gov/projects/tharsis/MOLA/SOFTWARE/
C uses the gravity field coefficients file 'mgm1025' available in
C the same ftp directory
        implicit none
        
        include "constants_mcd.inc"
        
c inputs:
        character*(*) dset ! path to datafiles
        real lon ! East longitude (degrees)
        real lat ! North latitude (degrees)
c output:
        real rareoid ! distance (in m) of areoid to center of Mars
        
c COMMON: (shared with readcs.F and geoid.F)
        integer ndeg,ndeg2,nd2p3
        parameter (ndeg=90,ndeg2=2*ndeg,nd2p3=ndeg2+3)
c data structure of gravity field
        double precision v0,omega,ae,gm
        double precision clm(0:ndeg,0:ndeg),slm(0:ndeg,0:ndeg)
        common /gmm1/v0,omega,ae,gm,clm,slm
        integer lmin,lmax
        double precision root(nd2p3)
        double precision requator
        common /sqr/ lmin,lmax,root,requator

        double precision dlon    ! double precision version of lon
        double precision dlat    ! double precision version of lat
        double precision dareoid ! double precision version of rareoid
!        double precision pi,d2r
!        parameter (pi=3.141592653589792D0, d2r=pi/180.d0)
        character*140 mgm     ! gravity field coefficients file
!        data mgm /'mgm1025'/
        integer llmax
        data llmax /50/
c flag for coefficients' initialization
        logical firstcall
        data firstcall /.true./
        
        if (firstcall) then
c initialize gravity coefficients
          lmax = llmax        
          ! build datafile path+name
          mgm=dset//'mgm1025'
          call readcs(mgm) ! read coefficients from file
          
          if((lmax.le.0).or.(lmax.gt.ndeg)) then
           if (output_messages) then
            write(out,*)'MOLAREOID Error: Gravity file in wrong format!'
            write(out,*)' lmin,max=',lmin,lmax
           endif
           stop
          endif
        firstcall=.false.
        endif ! of if (firstcall)

        dlon=dble(lon)
        dlat=dble(lat)
        call geoid(dlon,dlat,dareoid)
        rareoid=real(dareoid)
        
        end


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

        subroutine readcs(mgm)

        implicit none

        include "constants_mcd.inc"

c input
        character*(*) mgm ! file to read coefficients from

         character*64 title ! first line of file
         character*6 gcoef
c data structure of gravity field
         integer ndeg,ndeg2,nd2p3
          parameter (ndeg=90,ndeg2=2*ndeg,nd2p3=ndeg2+3)
          double precision c(0:ndeg,0:ndeg),s(0:ndeg,0:ndeg)
          double precision  v0,omega,ae,gm !,clm,slm
          double precision  plm(ndeg)
          common /gmm1/ v0,omega,ae,gm,c,s
          integer lmin,lmax
          double precision root(nd2p3)
          double precision r
          common /sqr/ lmin,lmax, root,r
c
        integer lcmax,mmax
        integer k,l,m
        double precision coef
        double precision x,xi,sum
        
          ae= 3396000.d0
          gm =42828.37d9
          omega=0.70882181d-4
c this value makes the equatorial mean radius equal to 3396 km.
        do k=1,nd2p3
         root(k)=sqrt(dble(k))
        enddo
c  initialize
        do l=0, ndeg
         do m=0,l
          c(l,m)=0.
          s(l,m)=0.
         enddo
        enddo
        open(unit=11,file=mgm,status='old',err=999)
        if (output_messages) then
          write(out,*)'Loading gravity field coefficients from file ',
     &                  trim(mgm)
        endif 
!        write(*,'(a,a)')' GCOEF potential model: ',mgm
        read (11,'(a)') title
!        write(out,*) title
        read (11,1000) gcoef,lcmax,mmax,gm,ae
!        write(out,*) gcoef,lcmax,mmax,gm,ae
!        write(out,*)
1        continue
        read(11,1000,end=99) gcoef,l,m,coef
1000   format(a6,8x,i3,i3,d24.14,d15.9)

c        read(1,1000,end=99) l,m,clm,slm
c1000   format(6x,i3,i3,1p,2d18.9)
        if(l.lt.0 .or. m.lt.0 .or. l.gt.ndeg .or. m.gt.l) goto 999
        if(gcoef.eq.'GCOEFC') c(l,m)=coef
        if(gcoef.eq.'GCOEFS') s(l,m)=coef
c        c(l,m)=clm
c        s(l,m)=slm
        goto 1
99        continue
        close(unit=11)
        r = 3396000.d0 !  MOLA potential surface
        xi = ae/r
        m=0
        x=0.
        sum=0.
        call lgndr(lmax,m,x,plm,root)
        do l=2, lmax
          sum = sum + c(l,m) *xi**l *plm(l-m+1)
        enddo
        sum = sum+1.
c centrifugal potential        
        v0=(gm*sum/r + 0.5*omega**2 * r**2)
!        write(out,*)'v0=', v0
        return

999        continue
        if (output_messages) then
         write(out,*)'READCS Error: file not found or incorrect data',
     &                lmax
         write(out,*)'title: ',title
         write(out,*)'gcoef,lcmax,mmax,gm,ae:',gcoef,lcmax,mmax,gm,ae
        endif
        end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

        subroutine geoid(dlon,dlat,rg)
c marsgeoid.f
c calculate radii of surface of constant gravity potential on rotating body
c given spherical harmonic potential of body at lat,lon,elevation
c
c W = V +Phi = gm/r [1+ (R/r)**n [harmonics]] -1/2 omega**2 r**2

        implicit none

        double precision dlon,dlat
        double precision rg
        
        double precision pi,d2r
        parameter (pi=3.141592653589792D0, d2r=pi/180.d0)
c data structure of gravity field
        integer ndeg,ndeg1,ndeg2,nd2p3
        parameter (ndeg=90,ndeg1=ndeg+1,ndeg2=2*ndeg,nd2p3=ndeg2+3)
        double precision v0,omega,ae,gm
        double precision c(0:ndeg,0:ndeg),s(0:ndeg,0:ndeg)
        common /gmm1/v0,omega,ae,gm,c,s
        integer lmin,lmax
        double precision root(nd2p3)
        double precision r
        common /sqr/ lmin,lmax,root,r
        double precision plm(ndeg1)
        double precision tol
        data tol /0.125d0/  ! tolerence on computed value of rg

        integer i,m,l
        double precision rlon,rlat
        double precision x,cslt,xi,sum,cslm
        double precision diff
c save r
        rg = r

        rlon=dlon*d2r
        rlat=dlat*d2r
        x = sin(rlat)
        cslt= cos(rlat)

c        csln= cos(rlon)
c        snln= sin(rlon)

         do i=1,8 ! usually 3 iterations suffice

          xi = ae/r
          sum = 0.
          do m=0, lmax
           call lgndr(lmax,m,x,plm,root)
           do l=m, lmax
             cslm =  c(l,m)*cos(m*rlon) + s(l,m)*sin(m*rlon)
             sum = sum + cslm*xi**l *plm(l-m+1)
           enddo
          enddo
          sum=sum+1.
c centrifugal potential        
          rg=(gm*sum+0.5*(omega**2)*(r**3)*(cslt**2))/v0
!          write(out,*) i,r,rg
          diff = r-rg
          r = rg
          if( abs(diff).lt. tol) goto 400
        enddo ! do i=1,8

 400        continue
c        write(*,500) dlon,dlat,r
c500        format(2f9.2, f12.3)
        end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      SUBROUTINE LGNDR(LMAX,M,X,PLM,SQR)
C Return vector of plm's of degree from m to lmax, for order m
      implicit none
      integer LMAX
      integer M
      double precision X
      double precision PLM(*),SQR(*)
      integer IFACT,I,LL
      double precision OM, PMM, PMMP1, PLL, SOMX2
C  needs sqrt of integers from 1 to 2*L+3
C  SIGN MODIF. OF NUM.REC. ALGORITHM, TO MATCH GEOPHYS. CONVENTION.
C  AND NORMALIZATION TO MEAN SQUARE UNITY.

      PMMP1=0 !dummy initialization to get rid of compiler warnings

      OM = 1.
      PMM=SQR(2*M+1)
      IF(M.gt.0) then
        PMM=SQR(2)*PMM
        SOMX2=SQRT((1.-X)*(1.+X))
        IFACT=1
        do I=1,M
          OM=-OM
          PMM=-PMM*(SQR(IFACT)/SQR(IFACT+1))*SOMX2
          IFACT=IFACT+2
        enddo
      ENDIF
      PLM(1)=OM*PMM
      if(LMAX.gt.M) then
        PMMP1=X*SQR(2*M+3)*PMM
        PLM(2)=OM*PMMP1
      endif
      if (LMAX.gt.M+1) then
         DO LL=M+2,LMAX
           PLL=SQR(2*LL+1)*SQR(2*LL-1)/(SQR(LL+M)*SQR(LL-M))
     & * (X*PMMP1-PMM*SQR(LL+M-1)*SQR(LL-M-1)/(SQR(2*LL-1)*SQR(2*LL-3)))
           PLM(LL-M+1)=OM*PLL
           PMM=PMMP1
           PMMP1=PLL
         ENDDO
      endif
      END

