geopotential_module_test Subroutine

    subroutine geopotential_module_test()

    use conversion_module
    use vector_module,  only: spherical_to_cartesian
    use random_module,  only: get_random_number

    implicit none

    character(len=*),parameter :: gravfile = './grav/GGM03C.GEO'  !! the coefficient file

    class(geopotential_model),pointer :: g
    type(geopotential_model_mueller)          ,target :: g_mueller
    type(geopotential_model_lear)             ,target :: g_lear
    type(geopotential_model_pines)            ,target :: g_pines
    type(geopotential_model_normalized_pines) ,target :: g_normalized_pines
    type(geopotential_model_kuga_carrara),target :: g_kuga_carrara

    real(wp),dimension(3) :: a1,a2,a3,a4,a5,rvec
    logical :: status_ok
    integer :: lat,lon,i,j,nmax,mmax
    real(wp) :: h,err1,err2,err3,err4,rlon,rlat,tmp
    character(len=20) :: name
    real :: tstart, tstop

    real(wp),dimension(3),parameter :: r = [0.1275627320e+05_wp, &
                                           0.1275627320e+05_wp, &
                                           0.1275627320e+05_wp ]   !! test case [km]

    integer,parameter :: n=9    !! degree
    integer,parameter :: m=9    !! order

    integer,parameter :: n_repeat = 1000000  !! number of time to repeat speed test

    write(*,*) ''
    write(*,*) '---------------'
    write(*,*) ' geopotential_module_test'
    write(*,*) '---------------'
    write(*,*) ''

    write(*,*) ''
    write(*,*) '*** number_of_coefficients routine ***'
    write(*,*) ''

    !tests:
    write(*,*) ''
    write(*,*) 6,6,number_of_coefficients(6,6)    ! 25
    write(*,*) 3,2,number_of_coefficients(3,2)    ! 6
    write(*,*) 4,0,number_of_coefficients(4,0)    ! 8
    write(*,*) 4,1,number_of_coefficients(4,1)    ! 9
    write(*,*) 0,0,number_of_coefficients(0,0)    ! -99 (error)
    write(*,*) ''

    write(*,*) ''
    write(*,*) '*** test case for the four methods ***'
    write(*,*) ''

    write(*,*) ''
    write(*,*) 'reading file: '//gravfile
    write(*,'(A,*(E30.16,1X))') 'r =',r

    call g_mueller%initialize(gravfile,n,m,status_ok)
    if (.not. status_ok) stop 'Error'
    call g_mueller%get_acc(r,n,m,a1)

    call g_lear%initialize(gravfile,n,m,status_ok)
    if (.not. status_ok) stop 'Error'
    call g_lear%get_acc(r,n,m,a2)

    call g_pines%initialize(gravfile,n,m,status_ok)
    if (.not. status_ok) stop 'Error'
    call g_pines%get_acc(r,n,m,a3)

    call g_kuga_carrara%initialize(gravfile,n,m,status_ok)
    if (.not. status_ok) stop 'Error'
    call g_kuga_carrara%get_acc(r,n,m,a4)

    call g_normalized_pines%initialize(gravfile,n,m,status_ok)
    if (.not. status_ok) stop 'Error'
    call g_normalized_pines%get_acc(r,n,m,a5)

    write(*,*) ''
    write(*,'(A,*(E30.16,1X))') 'mueller      =',a1
    write(*,'(A,*(E30.16,1X))') 'lear         =',a2
    write(*,'(A,*(E30.16,1X))') 'pines        =',a3
    write(*,'(A,*(E30.16,1X))') 'kuga_carrara =',a4
    write(*,'(A,*(E30.16,1X))') 'norm_pines   =',a3
    write(*,*) ''
    write(*,'(A,*(E30.16,1X))') 'mueller-lear difference      =',norm2(a1-a2)
    write(*,'(A,*(E30.16,1X))') 'mueller-pines difference     =',norm2(a1-a3)
    write(*,'(A,*(E30.16,1X))') 'mueller-normpines difference =',norm2(a1-a5)
    write(*,'(A,*(E30.16,1X))') 'lear-pines difference        =',norm2(a2-a3)
    write(*,'(A,*(E30.16,1X))') 'lear-kuga_carrara difference =',norm2(a2-a4)
    write(*,*) ''

    write(*,*) ''
    write(*,*) '*** accuracy test ***'
    write(*,*) ''

    h = 6778.0_wp    !radius magnitude
    status_ok = .true.
    do lat = -90, 90
        do lon = 0, 360

            rvec = spherical_to_cartesian(h,lon*deg2rad,lat*deg2rad)

            call g_mueller%get_acc(rvec,n,m,a1)          ! mueller
            call g_lear%get_acc(rvec,n,m,a2)             ! lear
            call g_pines%get_acc(rvec,n,m,a3)            ! pines
            call g_normalized_pines%get_acc(rvec,n,m,a5) ! normalized pines
            call g_kuga_carrara%get_acc(rvec,n,m,a4)     ! kuga/carrara

            err1 = norm2( a2 - a1 )
            err2 = norm2( a3 - a2 )
            err3 = norm2( a4 - a1 )
            err4 = norm2( a5 - a1 )
            if (err2>1.0e-15_wp .or. err1>1.0e-15_wp .or. err3>1.0e-15_wp .or. err4>1.0e-15_wp) then
                write(*,*) lat,lon,norm2(a1),norm2(a2),norm2(a2),norm2(a3),err1,err2,err3,err4
                status_ok = .false.
            end if
            if (abs(lat)==90) exit    !only do poles once

        end do
    end do
    call g_mueller%destroy()
    call g_lear%destroy()
    call g_pines%destroy()
    call g_normalized_pines%destroy()
    call g_kuga_carrara%destroy()
    if (status_ok) write(*,*) 'All tests passed.'

    write(*,*) ''
    write(*,*) '*** speed test ***'
    write(*,*) ''

    nmax = 10
    mmax = 10

    do i=1,5

        select case(i)
        case(1)
            g => g_mueller
            name = 'Mueller'
        case(2)
            g => g_lear
            name = 'Lear'
        case(3)
            g => g_pines
            name = 'Pines'
        case(4)
            g => g_kuga_carrara
            name = 'Kuga/Carrara'
        case(5)
            g => g_normalized_pines
            name = 'Normalized Pines'
        end select
        call g%initialize(gravfile,nmax,mmax,status_ok)
        if (.not. status_ok) stop 'Error'

        call random_seed()
        tmp = zero
        call cpu_time(tstart)
        do j=1,n_repeat

             h    = get_random_number(6778.0_wp, 10000.0_wp)
             rlon = get_random_number(0.0_wp, 360.0_wp)
             rlat = get_random_number(-90.0_wp, 90.0_wp)

             rvec = spherical_to_cartesian(h,rlon*deg2rad,rlat*deg2rad)
             call g%get_acc(rvec,nmax,mmax,a1)
             tmp = tmp + norm2(a1)

        end do
        call cpu_time(tstop)

        call g%destroy()

        write(*,'(A10,1X,E30.16,1X,F13.6,1X,A)') trim(name), tmp, tstop-tstart, 'sec'

    end do

    end subroutine geopotential_module_test