read_geopotential_file Subroutine

    subroutine read_geopotential_file(me,filename,nmax,mmax,status_ok)

    implicit none

    class(geopotential_model),intent(inout) :: me
    character(len=*),intent(in)        :: filename
    integer,intent(in)                 :: nmax
    integer,intent(in)                 :: mmax
    logical,intent(out)                :: status_ok

    character(len=100) :: c1,c2,fmt    ! for reading the format statements
    integer  :: iunit,istat,i1,i2,nc,i,ii,n,m
    real(wp) :: d1,d2,d3,d4,f1,t1,t2

    call me%destroy()

    if (nmax>0 .and. mmax>0 .and. mmax<=nmax) then

        status_ok = .true.
        nc = number_of_coefficients(nmax,mmax)
        me%nmax = nmax
        me%mmax = mmax
        me%filename = trim(filename)

        open(newunit=iunit,file=me%filename,status='OLD',iostat=istat)

        if (istat==0) then

            !size the coefficient arrays:
            select type (me)
            class is (geopotential_model_vector_coeff)

                allocate(me%c(nc))    !they are stored compressed into arrays
                allocate(me%s(nc))    !
                me%c = zero
                me%s = zero

            class is (geopotential_model_matrix_coeff)

                !Note: will get all the nmax x nmax coefficients, even if mmax < nmax

                allocate(me%cnm(nmax,0:nmax))    !probably could replace with 2:nmax
                allocate(me%snm(nmax,0:nmax))
                me%cnm = zero
                me%snm = zero

            class default

                write(*,*) 'ERROR: INVALID geopotential_model CLASS!'
                status_ok = .false.
                return

            end select

            !Read the file:

            ! (2A10,2E20.10)                                               J2-DOT = -26x10-12
            !EGM 96                  398600.44150E+09          6378136.30
            ! (A6,2I3,2D19.12,2D13.6,F4.0)
            !RECOEF  2  0-0.484165371736E-03 0.000000000000E+00 0.356106E-10 0.000000E+00 -1.
            !HONKCR  2  0-0.484169544736E-03 0.000000000000E+00 0.356106E-10 0.000000E+00 -1.
            !IERS    2  1-0.186987640000E-09 0.119528010000E-08 0.100000E-29 0.100000E-29 -1.
            !RECOEF  2  2 0.243914352398E-05-0.140016683654E-05 0.537392E-10 0.543533E-10 -1.
            !RECOEF  3  0 0.957254173792E-06 0.000000000000E+00 0.180942E-10 0.000000E+00 -1.
            !RECOEF  3  1 0.202998882184E-05 0.248513158716E-06 0.139652E-09 0.136459E-09 -1.

            read(iunit,'(A)',iostat=istat) c1    !this is the FMT statement for the next line
            call get_format_statement(c1,fmt)    !

            read(iunit,trim(fmt),iostat=istat) c1,c2,d1,d2

            me%name = trim(c1)//trim(c2)
            me%mu = d1/1000.0_wp**3    !km3/s2
            me%re = d2/1000.0_wp       !km

            read(iunit,'(A)',iostat=istat) c1    !this is the FMT statement for the next lines
            call get_format_statement(c1,fmt)    !

            do i=1,nc    !...until end of file or all the coefficients have been read...

                read(iunit,trim(fmt),iostat=istat) c1,i1,i2,d1,d2,d3,d4,f1

                if (istat > 0)  then

                    write(*,*) 'Error reading file:' //trim(filename)
                    call me%destroy()
                    status_ok = .false.
                    exit

                else if (istat < 0) then

                    ! end of file:
                    if (i>nc) then
                        write(*,*) 'Error: not enough coefficients in file.'
                        call me%destroy()
                        status_ok = .false.
                        exit
                    end if

                end if

                select type (me)
                class is (geopotential_model_vector_coeff)
                    me%c(i) = d1
                    me%s(i) = d2
                class is (geopotential_model_matrix_coeff)
                    me%cnm(i1,i2) = d1
                    me%snm(i1,i2) = d2
                end select

               end do

            close(iunit,iostat=istat)

        else
            write(*,*) 'Error reading file: '//trim(filename)
            call me%destroy()
            status_ok = .false.
        end if

    else
        write(*,*) 'Error: invalid n,m values: ',nmax,mmax
        call me%destroy()
        status_ok = .false.
    end if

    !unnormalize the coefficients if necessary:
    if (status_ok) then

        select type (me)
        class is (geopotential_model_vector_coeff)

            !for this one, the coefficients are stored in arrays

            ii = 0    !counter
            do n = 2, me%nmax    ! based on Lear's CONVERT routine
                t1 = 2*n+1
                ii=ii+1
                me%c(ii) = sqrt(t1)*me%c(ii)
                me%s(ii) = zero
                do m = 1, n
                    ii=ii+1
                    t2 = sqrt(FL(n-m)*t1*two / FL(n+m))
                    me%c(ii) = t2*me%c(ii)
                    me%s(ii) = t2*me%s(ii)
                end do
            end do

        class is (geopotential_model_matrix_coeff)

            !for this one, the coefficients are stored in matrices

            select type (me)
            class is (geopotential_model_kuga_carrara)
                !this model uses the normalized coefficients
                return
            class is (geopotential_model_normalized_pines)
                !this model uses the normalized coefficients
                return
            end select

            call convert(me%nmax,me%cnm,me%snm)

        end select

    end if

    end subroutine read_geopotential_file