text:
      MODULE dijkl1_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.4G  10:47:09  03/09/06
      SUBROUTINE dijkl1 (C, N, NATI, W, CIJ, WCIJ, CKL, XY)
      USE vast_kind_param,ONLY: DOUBLE
      use meci_C, only : nmos
      INTEGER, INTENT(IN) :: N
      INTEGER, INTENT(IN) :: NATI
      REAL(DOUBLE), DIMENSION(N,*), INTENT(IN) :: C
      REAL(DOUBLE), DIMENSION(*) :: W
      REAL(DOUBLE), DIMENSION(*), INTENT(OUT) :: CIJ
      REAL(DOUBLE), DIMENSION(*), INTENT(INOUT) :: WCIJ
      REAL(DOUBLE), DIMENSION(*), INTENT(INOUT) :: CKL
      REAL(DOUBLE), DIMENSION(NMOS,NMOS,NMOS,NMOS), INTENT(OUT) :: XY
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
      SUBROUTINE fftfax_g(n, ifax, trigs)
      USE stel_kinds
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER n
      INTEGER, DIMENSION(13) :: ifax
      REAL(rprec), DIMENSION(*) :: trigs
#if !defined(CRAY) || defined(LONESTAR) || defined(MCURIE)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: mode=3, i
!-----------------------------------------------
!
! mode 3 is used for REAL/half-complex transforms.  it is possible
! to do complex/complex transforms with other values of mode, but
! documentation of the details were not available when this routine
! was written.
!
      CALL fax (ifax, n, mode)
      i = ifax(1)
      IF (ifax(i+1)>5 .or. n<=4) ifax(1) = -99
      IF (ifax(1) <= 0) STOP 'IFAX(1) <= 0 in fftfax'
      CALL fftrig_g (trigs, n, mode)
#else
      CALL fftfax (n, ifax, trigs)
#endif
      END SUBROUTINE fftfax_g

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: ORNL-Fusion/LIBSTELL


text:
  program test_kerr

      use geodesics
      use class_four_Vector
      use kerr, only: calc_polar_psi, comoving_ortho_debug, comoving_ortho, calc_kappapw
      use grtrans_inputs
      use phys_constants, GC=>G
      use fluid_model

      implicit none

      real(kind=8), dimension(:), allocatable :: r, th
      real(kind=8) :: q2,l,alpha,beta,su,sm
      real(kind=8), dimension(:), allocatable :: s2xi,c2xi,c2psi,s2psi,ang,rshift,cosne, &
           aat,aar,aath,aaph,kht,khr,khth,khph,bhatt,bhatr,bhatth,bhatph, cosne2
      real(kind=8), dimension(:,:), allocatable :: aahat,bbhat
      type (four_Vector), dimension(:), allocatable :: uhat,bhat,khat,aa,kdifft
      real(kind=8) :: maxang, maxangdiff, maxgeodiff,maxkbdiff
      integer :: gunit, status
      type (geo) :: g
      type (geokerr_args) :: gargs
      type (fluid_args) :: fargs
      type (fluid) :: f
      character(len=500) :: file='inputs.in'

      call read_inputs(file)
      call assign_fluid_args(fargs,fdfile,fhfile,fgfile,fsim,fnt,findf,fnfiles,fjonfix, &
            fnw,fnfreq_tab,fnr,foffset,fdindf,fmagcrit,frspot,fr0spot,fn0spot,ftscl,frscl, &
            fwmin,fwmax,ffmin,ffmax,frmax,fsigt,ffcol,fmdot,mbh,fnscl,fnnthscl,fnnthp,fbeta,&
            fbl06,fnp,ftp,frin,frout,fthin,fthout,fphiin,fphiout,fscalefac)

      write(6,*) 'inputs read'
      call load_fluid_model(fname,spin,fargs)
!      call initialize_pixels(gargs,.true.,standard,mu0,spin,uout,rcut,
!     & nrotype,a1,a2,b1,b2,nro,nphi,nup)
      write(6,*) 'fluid loaded', standard, mu0(1), spin, uout, uin, rcut, nrotype
      write(6,*) 'fluid: ',a1,a2,b1,b2,nro,nphi,nup
      call initialize_geokerr_args(gargs,nro*nphi)
!       write(6,*) 'init pixels: ',allocated(gargs%muf)
      call initialize_pixels(gargs,.true.,standard,mu0(1),phi0,spin,uout,uin, &
        rcut,nrotype,a1,a2,b1,b2,nro,nphi,nup)
      write(6,*) 'pixels'
      call initialize_geodesic(g,gargs,6258,status)
      write(6,*) 'geo'
      call initialize_fluid_model(f,fname,spin,g%npts)
      write(6,*) 'init fluid', size(g%x), g%gk%a
      call get_fluid_vars(g%x,g%k,g%gk%a,f)
      write(6,*) 'kdotk: ',maxval(g%k*g%k)
      write(6,*) g%gk%alpha(1), g%gk%beta(1)
!      r=g%x%data(2); th=g%x%data(3)
      allocate(s2xi(g%npts)); allocate(c2xi(g%npts)); allocate(ang(g%npts))
      allocate(rshift(g%npts)); allocate(aahat(g%npts,3))
      allocate(kdifft(g%npts-1)); allocate(aa(g%npts))
      allocate(s2psi(g%npts)); allocate(c2psi(g%npts)); allocate(cosne(g%npts))
      allocate(aat(g%npts)); allocate(aar(g%npts)); allocate(cosne2(g%npts))
      allocate(aath(g%npts)); allocate(aaph(g%npts)); allocate(bhatt(g%npts))
      allocate(bhatr(g%npts)); allocate(bhatth(g%npts)); allocate(bhatph(g%npts))
      allocate(kht(g%npts)); allocate(khr(g%npts)); allocate(khth(g%npts))
      allocate(khph(g%npts)); allocate(bhat(g%npts)); allocate(khat(g%npts))
      call comoving_ortho_debug(g%x%data(2),g%x%data(3),g%gk%a,g%gk%alpha(1), &
       g%gk%beta(1),g%gk%mu0,f%u,f%b,g%k,s2xi,c2xi,ang,rshift,cosne2,aahat,aat,aar, &
       aath,aaph,kht,khr,khth,khph,bhatt,bhatr,bhatth,bhatph)
      khat%data(1)=kht; khat%data(2)=khr; khat%data(3)=khth; khat%data(4)=khph
      bhat%data(1)=bhatt; bhat%data(2)=bhatr; bhat%data(3)=bhatth; bhat%data(4)=bhatph
      aa%data(1)=aat; aa%data(2)=aar; aa%data(3)=aath; aa%data(4)=aaph
      call calc_polar_psi(g%x%data(2),cos(g%x%data(3)),g%gk%q2,g%gk%a,g%gk%alpha, &
           g%gk%beta,rshift,g%gk%mu0,g%k,c2psi,s2psi,cosne)
! Test difference between psi, chi:
      maxang=maxval(abs(c2psi-c2xi))
      write(6,*) 'test kerr xi: ',c2xi
      write(6,*) 'test kerr psi: ',c2psi
      maxkbdiff=maxval(abs(cosne2-cosne))
! Test chi:
      maxangdiff=maxval(abs(c2xi*c2xi+s2xi*s2xi-1d0))
      write(6,*) 'diffs: ',maxang,maxangdiff,maxval(abs(g%k*g%k))
!      write(6,*) 'chi: ', c2xi
!      write(6,*) 'psi: ',c2psi
!      write(6,*) 'ang: ',cos(ang)
!      write(6,*) 'cosne: ',cosne
!      write(6,*) 'cosne2: ',g%k*f%b / (g%k*f%u)
!      write(6,*) 'cosne2: ',g%k(1)%data
!      write(6,*) 'cosne2: ',f%b(1)%data
!      write(6,*) 'cosne2: ',f%u(1)%data
!      write(6,*) 'cosne3: ',sqrt(g%gk%q2)*rshift/g%x%data(2), g%x%data(2), rshift, g%gk%q2
 ! write to file
! the definition in calc_kappapw is K1 - i K2 while in my transport_perpk it is K1 - i K2, so compare -K2 to imaginary part
      write(6,*) 'kappa: ',calc_kappapw(g%gk%a,g%x%data(2),cos(g%x%data(3)),g%k,aa), &
           g%gk%alpha+g%gk%a*sqrt(1.-g%gk%mu0**2.), -(-g%gk%beta)
      open(unit=12,file='unit_test_kerr.out')
      write(12,*) maxang, maxangdiff, maxval(abs(g%k*g%k)), &
           maxval(abs(khat%data(2)*aahat(:,1)+khat%data(3)*aahat(:,2) &
           +khat%data(4)*aahat(:,3))), &
           maxval(abs(g%gk%beta-aimag(calc_kappapw(g%gk%a, &
           g%x%data(2),cos(g%x%data(3)),g%k,aa))))
      close(unit=12)
      open(unit=12,file='test_kerr.out')
      write(12,*) g%gk%alpha(1), g%gk%beta(1), g%gk%a, g%gk%mu0, g%gk%sm, g%gk%su, g%npts
      write(12,*) g%x%data(1)
      write(12,*) g%x%data(2)
      write(12,*) g%x%data(3)
      write(12,*) g%x%data(4)
      write(12,*) g%tprarr
      write(12,*) g%tpmarr
      write(12,*) g%k%data(1)
      write(12,*) g%k%data(2)
      write(12,*) g%k%data(3)
      write(12,*) g%k%data(4)
      write(12,*) f%u%data(1)
      write(12,*) f%u%data(2)
      write(12,*) f%u%data(3)
      write(12,*) f%u%data(4)
      write(12,*) f%b%data(1)
      write(12,*) f%b%data(2)
      write(12,*) f%b%data(3)
      write(12,*) f%b%data(4)
      write(12,*) ang
      write(12,*) rshift
      write(12,*) s2xi
      write(12,*) c2xi
      write(12,*) aa%data(2)
      write(12,*) aa%data(3)
      write(12,*) aa%data(4)
      write(12,*) aahat(:,1)
      write(12,*) aahat(:,2)
      write(12,*) aahat(:,3)
      write(12,*) bhat%data(1)
      write(12,*) bhat%data(2)
      write(12,*) bhat%data(3)
      write(12,*) bhat%data(4)
      write(12,*) khat%data(1)
      write(12,*) khat%data(2)
      write(12,*) khat%data(3)
      write(12,*) khat%data(4)
      close(unit=12)
      call del_geodesic(g); call del_fluid_model(f)
      deallocate(s2xi); deallocate(c2xi); deallocate(ang); deallocate(rshift)
      deallocate(aahat); deallocate(aa); deallocate(bhat); deallocate(khat)
      deallocate(c2psi); deallocate(s2psi); deallocate(kdifft); deallocate(cosne)
      deallocate(aat); deallocate(aar); deallocate(aath); deallocate(aaph)
      deallocate(bhatt); deallocate(bhatr); deallocate(bhatth); deallocate(bhatph)
      deallocate(kht); deallocate(khr); deallocate(khth); deallocate(khph)
      deallocate(cosne2)

      end program test_kerr

max_stars_repo_name: achael/grtrans

max_issues_repo_name: achael/grtrans

max_forks_repo_name: achael/grtrans


text:
 General Info

Im an undergraduate student majoring in linguistics whos nearly done at UCD (graduates in 2006).

I have studied Japanese for more than four years now and am minoring in it.  I studied in Japan in the summer of 2005.

I was born in Fort Worth, Texas, but lived in Redding, California for most of my life before coming to Davis.  I transfered to Davis from Shasta College (in Redding).

 Pictures

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!##############################################################################
!# ****************************************************************************
!# <name> LinearSolverAutoInitialise </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module contains routines to initialise a linear solver for a given
!# situation. On one hand, it provides initialisation routines for standard
!# solver settings. On the other hand, there is a parser included which can
!# read user-defined text files from hard disc that configure the setting of
!# a solver.
!#
!# The following routines can be found here:
!#
!# 1.) linsolinit_initFromFile
!#     -> Initialise a linear solver node by reading parameters from a
!#        INI/DAT file
!#
!# 2.) linsolinit_initParams
!#     -> Read parameters of a linear solver from a section of a DAT file
!#        and write them into a solver structure.
!#
!# </purpose>
!##############################################################################

module linearsolverautoinitialise

!$ use omp_lib
  use fsystem
  use genoutput
  use spatialdiscretisation
  use linearsolver
  use multilevelprojection
  use filtersupport
  use paramlist

  implicit none

  private

  public :: linsolinit_CC2DMultigrid
  public :: linsolinit_initFromFile
  public :: linsolinit_initParams

contains

  ! ***************************************************************************

!<subroutine>

  subroutine linsolinit_CC2DMultigrid(p_rsolverNode, isolverType, nlevels, &
                        niterationsMin, niterationsMax,daccuracyAbs,daccuracyRel,  &
                        nprecSteps, nsmoothingSteps, &
                        nmaxCoarseGridSteps, dcoarseGridAccuracyAbs, &
                        dcoarseGridAccuracyRel)

  !<description>

  ! This routine builds a solver node for the linear solver for CC2D-like
  ! applications. The following combinations are supported:
  ! 1.) Multigrid solver, VANKA smoother, VANKA coarse grid solver
  ! 2.) BiCGStab-solver, Multigrid preconditioner,
  !     VANKA smoother, VANKA coarse grid solver
  ! The caller must attach the matrices of all levels to the solver manually
  ! by calling linsol_setMatrices. Afterwards the linear solver can be started.

  !</description>

  !<input>

  ! Type of solver structure which should be build up.
  ! 1 = Multigrid solver, VANKA smoother, VANKA coarse grid solver
  ! 2 = BiCGStab-solver, Multigrid preconditioner, VANKA smoother,
  !     VANKA coarse grid solver
  integer, intent(in)               :: isolverType

  ! Number of levels
  integer, intent(in)               :: nlevels

  ! Minimum number of solver iterations
  integer, intent(in)               :: niterationsMin

  ! Maximum number of solver iterations
  integer, intent(in)               :: niterationsMax

  ! If solver conbination is BiCGStab with MG preconditioning, number of
  ! steps multigrid should perform. Otherwise ignored
  integer, intent(in)               :: nprecSteps

  ! Absolute accuracy of the solver
  real(DP), intent(in)                   :: daccuracyAbs

  ! Relativew accuracy of the solver
  real(DP), intent(in)                   :: daccuracyRel

  ! Number of pre- and postsmoothing steps
  integer, intent(in)               :: nsmoothingSteps

  ! Absolute accuracy on the coarse grid
  real(DP), intent(in)                   :: dcoarseGridAccuracyAbs

  ! Relative accuracy on the coarse grid
  real(DP), intent(in)                   :: dcoarseGridAccuracyRel

  ! Maximum iterations on the coarse grid
  integer, intent(in)               :: nmaxCoarseGridSteps

  ! A list of spatial discretisation structures for the three equations
  ! that are supported by CC2D: x-velocity, y-velocity, pressure
  type(t_spatialDiscretisation), dimension(3) :: RspatialDiscr

!</input>

!<output>

  ! The solver node that identifies the solver. The caller must attach
  ! level-dependent data (matrix information) to it with the standard
  ! routines in the module LinearSolver. Afterwards the structure
  ! can be used to solve the problem.

  type(t_linsolNode),pointer :: p_rsolverNode

!</output>

!</subroutine>

  ! local variables
  type(t_linsolNode),pointer :: p_rmgSolver
  integer               :: ilevel
  type(t_linsolNode),pointer :: p_rpreSmoother
  type(t_linsolNode),pointer :: p_rpostSmoother
  type(t_linsolNode),pointer :: p_rcoarseGridSolver
  type(t_linsolMGLevelInfo), pointer :: p_rlevelInfo
  type(t_interlevelProjectionBlock) :: rprojection

  ! Create the solver node - either BiCGStab or Multigrid.
  ! If we create BiCGStab, attach multigrid as preconditioner.

  call linsol_initMultigrid (p_rmgSolver)

  if (isolverType .eq. 1) then

    p_rsolverNode => p_rmgSolver

  else

    call linsol_initBiCGStab (p_rsolverNode,p_rmgSolver)

    ! Configure MG preconditioning as a fixed number of MG steps
    ! without checking the residuum.

    p_rmgSolver%nminIterations = nprecSteps
    p_rmgSolver%nmaxIterations = nprecSteps
    p_rmgSolver%iresCheck      = NO

  end if

  ! Configure the solver
  p_rsolverNode%nminIterations = niterationsMin
  p_rsolverNode%nmaxIterations = niterationsMax
  p_rsolverNode%depsRel = daccuracyRel
  p_rsolverNode%depsAbs = daccuracyAbs

  ! Initialise a standard interlevel projection structure for all levels.
  call mlprj_initProjectionDirect (rprojection,RspatialDiscr)

  ! Continue to configure MG by accessing p_rmgSolver.
  ! Loop through the levels.

  do ilevel = 1,nlevels

    ! On the lowest level create a coarse grid solver structure
    if (ilevel .eq. 1) then
      call linsol_initVANKA (p_rcoarseGridSolver)
      p_rcoarseGridSolver%depsRel = dcoarseGridAccuracyRel
      p_rcoarseGridSolver%depsAbs = dcoarseGridAccuracyAbs
      p_rcoarseGridSolver%nmaxIterations = nmaxCoarseGridSteps
    else
      nullify(p_rcoarseGridSolver)
    end if

    ! Create pre- and postsmoother structure on the current level
    call linsol_initVANKA (p_rpreSmoother)
    call linsol_initVANKA (p_rpostSmoother)

    ! Configure the structures to form a smoother. A smoother is a solver
    ! that iterates a finite number of times without respecting the
    ! residuum.

    call linsol_convertToSmoother (p_rpreSmoother,nsmoothingSteps)
    call linsol_convertToSmoother (p_rpostSmoother,nsmoothingSteps)

    ! Create the level, attrach it to the solver and proceed to the
    ! next level
    call linsol_addMultigridLevel (p_rlevelInfo,p_rmgSolver, rprojection,&
                    p_rpresmoother,p_rpostsmoother,p_rcoarseGridSolver)

  end do

  end subroutine

  ! ***************************************************************************

!<subroutine>

  recursive subroutine linsolinit_initFromFile (p_rsolverNode, rparamList, &
                                                ssolverName, nlevels, RfilterChain, &
                                                rinterlevelProjection)

!<description>
  ! This routine creates a new linear solver node p_rsolverNode of the
  ! heap. The configuration on this solver is set up according to the
  ! parameters in the parameter list rparamList (which can be read e.g.
  ! from an INI/DAT file). The string ssolverName specifies the name of a
  ! section in the parameter list rparamList that serves as 'head' of the
  ! solver, i.e. that defines the 'main' linear solver.
  ! The routine automatically creates all sub-solvers (preconditioners,
  ! smoothers,...) by evaluating the parameters in rparamList,
!</description>

!<input>
  ! The parameter list that contains the whole solver configuration.
  type(t_parlist), intent(in) :: rparamList

  ! The name of a section in rparamList that contains the configuration of
  ! the main linear solver.
  character(LEN=*), intent(in) :: ssolverName

  ! Number of levels in the discretisation.
  integer, intent(in) :: nlevels

  ! OPTIONAL: A filter chain (i.e. an array of t_filterChain
  ! structures) if filtering should be applied to the vector during the
  ! iteration. If not, no filtering will be used.
  ! The filter chain (i.e. the array) must exist until the system is solved!
  ! The filter chain must be configured for being applied to defect vectors.
  type(t_filterChain), dimension(:), intent(in), target, optional :: RfilterChain

  ! OPTIONAL: An interlevel projection structure that configures the projection
  ! between the solutions on a finer and a coarser grid. The structure
  ! must have been initialised with mlprj_initProjection.
  !
  ! Note that this structure is level-independent (as long as the
  ! spatial discretisation structures on different levels are 'compatible'
  ! what they have to be anyway), so the same structure can be used
  ! to initialise all levels!
  !
  ! The structure must be present if any kind of Multigrid solver is used.
  ! If the application knowns that there is no MG solver, the parameter
  ! can be ommitted.
  type(t_interlevelProjectionBlock), optional :: rinterlevelProjection
!</input>

!<output>
  ! A pointer to a new linear solver node on the heap, initialised
  ! by the configuration in the parameter list.
  type(t_linsolNode), pointer :: p_rsolverNode
!</output>

!</subroutine>

    ! local variables
    type(t_parlstSection), pointer :: p_rsection
    character(LEN=SYS_STRLEN) :: spreconditioner,spresmoother, spostsmoother
    character(LEN=SYS_STRLEN) :: scoarsegridsolver,sString
    integer :: isolverType,isolverSubtype
    type(t_linsolNode), pointer :: p_rpreconditioner,p_rpresmoother,p_rpostsmoother
    type(t_linsolNode), pointer :: p_rcoarsegridsolver
    integer :: i1,ilev,ikrylowDim
    real(DP) :: d1
    type(t_filterChain), dimension(:), pointer :: p_Rfilter
    type(t_linsolMGLevelInfo), pointer     :: p_rlevelInfo
    type(t_linsolMG2LevelInfo), pointer     :: p_rlevelInfo2

    nullify(p_Rfilter)
    if (present(RfilterChain)) then
      p_Rfilter => RfilterChain
    end if

    ! Check that there is a section called ssolverName - otherwise we
    ! cannot create anything!

    call parlst_querysection(rparamList, ssolverName, p_rsection)

    if (.not. associated(p_rsection)) then
      call output_line ('Cannot create linear solver; no section ''' &
          // trim(ssolverName) // '''!', OU_CLASS_ERROR, OU_MODE_STD, &
          'linsolinit_initFromFile')
      call sys_halt()

    end if

    ! Ok, we have the section where we can get parameters from.
    ! Get the solver type we should set up.
    ! Let us hope the parameter 'isolverType' exists! That one is
    ! mandatory; if not, the get-routine will stop.
    call parlst_getvalue_int (p_rsection, 'isolverType', isolverType)

    ! Try to get the solver subtype from the parameter list.
    ! This allows switching between different variants of the same
    ! basic algorithm (e.g. VANKA)
    call parlst_getvalue_int (p_rsection, 'isolverSubtype', isolverSubtype,0)

    ! Many solvers support preconditioners - we try to fetch the
    ! name of the preconditioner in-advance to prevent code
    ! duplication
    call parlst_getvalue_string (p_rsection, 'spreconditioner', sString,'')
    spreconditioner = ''
    if (sString .ne. '') read (sString,*) spreconditioner

    ! Now comes a biiig select for all the different types of solvers
    ! that are supported by this routine.
    select case (isolverType)

    case (LINSOL_ALG_DEFCORR)
      ! Defect correction
      !
      ! Initialise a solver node for the preconditioner - if there is one.
      nullify(p_rpreconditioner)
      if (spreconditioner .ne. '') then
        call linsolinit_initFromFile (p_rpreconditioner,rparamList,&
                                      spreconditioner,nlevels,RfilterChain)
      end if
      ! Init the solver node
      call linsol_initDefCorr (p_rsolverNode,p_rpreconditioner,p_Rfilter)

    case (LINSOL_ALG_JACOBI)
      ! Jacobi solver
      !
      ! Init the solver node
      call linsol_initJacobi (p_rsolverNode)

    case (LINSOL_ALG_SOR)
      ! SOR/GS solver
      !
      ! Init the solver node. domega is set to 1.0 as standard
      ! to activate GS. The actual domega is initialised later if specified
      ! in the DAT file.
      call linsol_initSOR (p_rsolverNode,1.0_DP)

    case (LINSOL_ALG_SSOR)
      ! SSOR solver
      !
      ! Init the solver node.
      ! Parameters:
      !  iscale = 1   -> Scale preconditioned vector according to actual formula
      !                  in the literature
      !         = 0   -> no scaling, original FEAT implementation
      call parlst_getvalue_int (p_rsection, 'bscale', i1, 0)
      call linsol_initSSOR (p_rsolverNode,bscale = i1 .eq. 1)

    case (LINSOL_ALG_BICGSTAB)
      ! BiCGStab solver
      !
      ! Initialise a solver node for the preconditioner - if there is one.
      nullify(p_rpreconditioner)
      if (spreconditioner .ne. '') then
        call linsolinit_initFromFile (p_rpreconditioner,rparamList,&
                                      spreconditioner,nlevels,RfilterChain)
      end if

      ! Init the solver node
      call linsol_initBiCGStab (p_rsolverNode,p_rpreconditioner,p_Rfilter)

    case (LINSOL_ALG_GMRES)
      ! GMRES solver
      !
      ! Initialise a solver node for the preconditioner - if there is one.
      nullify(p_rpreconditioner)
      if (spreconditioner .ne. '') then
        call linsolinit_initFromFile (p_rpreconditioner,rparamList,&
                                      spreconditioner,nlevels,RfilterChain)
      end if

      ! Try to get the solver subtype from the parameter list.
      ! This allows switching between right- and left preconditioned BiCGStab.
      call parlst_getvalue_int (p_rsection, 'isolverSubtype', isolverSubtype,0)

      ! Krylow space dimension
      call parlst_getvalue_int (p_rsection, 'ikrylovDim', ikrylowDim,40)

      ! Apply Gram Schmidt twice
      call parlst_getvalue_int (p_rsection, 'itwiceGS', i1,0)

      ! Init the solver node
      call linsol_initGMRES (p_rsolverNode,ikrylowDim,p_rpreconditioner,p_Rfilter,&
        i1 .eq. 1)

    case (LINSOL_ALG_UMFPACK4)
      ! UMFPACK4 solver
      !
      ! Init the solver node
      call linsol_initUMFPACK4 (p_rsolverNode)

    case (LINSOL_ALG_ILU0)
      ! ILU0 solver
      !
      ! Init the solver node
      call linsol_initILU0 (p_rsolverNode)

    case (LINSOL_ALG_MILUS1x1)
      ! (M)ILU solver
      !
      ! Parameters:
      !  ifillinLevel = 0 / 1 / 2 / ...   -> Fill-in level for factorisation
      !  drelax       = 0.0               -> Build ILU(s)
      !               > 0.0               -> Build MILU(s)
      !
      ! Get fill-in level
      call parlst_getvalue_int (p_rsection, 'ifill', i1, 0)

      ! Get MILU relaxsation parameter
      call parlst_getvalue_double (p_rsection, 'drelax', d1, 0.0_DP)

      ! Init the solver node
      call linsol_initMILUs1x1 (p_rsolverNode,i1,d1)

    case (LINSOL_ALG_VANKA)
      ! VANKA solver
      !
      ! Init the solver node
      call linsol_initVANKA (p_rsolverNode,1.0_DP,isolverSubtype)

    case (LINSOL_ALG_MULTIGRID)
      ! Multigrid solver
      !
      ! Parameters:
      !  icycle         = 0               -> F-cycle
      !                 = 1               -> V-cycle
      !                 = 2               -> W-cycle
      !  dalphaMin      >= 0.0            -> minimum damping parameter; standard = 1.0
      !  dalphaMin      >= 0.0            -> maximum damping parameter; standard = 1.0
      !  spreSmootherName                 -> Name of the presmoother section
      !  spostSmootherName                -> Name of the postsmoother section
      !  scoarseGridSolver                -> Name of the coarse grid solver section
      !
      ! Ok, that is the most complicated one :-)
      ! At first, initialise the solver:
      call linsol_initMultigrid (p_rsolverNode,p_Rfilter)

      ! Then, get solver specific data.
      call parlst_getvalue_int (p_rsection, 'icycle', &
                                p_rsolverNode%p_rsubnodeMultigrid%icycle,&
                                p_rsolverNode%p_rsubnodeMultigrid%icycle)

      ! Coarse grid correction parameters
      call parlst_getvalue_int (p_rsection, 'ccorrectionTypeAlpha', &
           p_rsolverNode%p_rsubnodeMultigrid%rcoarseGridCorrection%ccorrectionType,&
           p_rsolverNode%p_rsubnodeMultigrid%rcoarseGridCorrection%ccorrectionType)

      call parlst_getvalue_double (p_rsection, 'dalphaMin', &
           p_rsolverNode%p_rsubnodeMultigrid%rcoarseGridCorrection%dalphaMin,&
           p_rsolverNode%p_rsubnodeMultigrid%rcoarseGridCorrection%dalphaMin)

      call parlst_getvalue_double (p_rsection, 'dalphaMax', &
           p_rsolverNode%p_rsubnodeMultigrid%rcoarseGridCorrection%dalphaMax,&
           p_rsolverNode%p_rsubnodeMultigrid%rcoarseGridCorrection%dalphaMax)

      ! Do we have a presmoother?
      call parlst_getvalue_string (p_rsection, 'spreSmootherName', sString,'')
      spresmoother = ''
      if (sString .ne. '') read (sString,*) spresmoother

      ! A postsmoother?
      call parlst_getvalue_string (p_rsection, 'spostSmootherName', sString,'')
      spostsmoother = ''
      if (sString .ne. '') read (sString,*) spostsmoother

      ! A coarse grid solver?
      call parlst_getvalue_string (p_rsection, 'scoarseGridSolver', sString,'')
      scoarsegridsolver = ''
      if (sString .ne. '') read (sString,*) scoarsegridsolver

      ! We need an interlevel projection structure!
      if (.not. present(rinterlevelProjection)) then
        call output_line ('Cannot create linear solver; ' // &
            'no interlevel projection structure!', OU_CLASS_ERROR, OU_MODE_STD, &
            'linsolinit_initFromFile')
        call sys_halt()
      end if

      ! Initialise the coarse grid solver - if there is one. There must be one!
      if (scoarsegridsolver .eq. '') then
        call output_line ('Cannot create linear solver; ' // &
            'no coarse grid solver for MG!', OU_CLASS_ERROR, OU_MODE_STD, &
            'linsolinit_initFromFile')
        call sys_halt()
      end if
      call linsolinit_initFromFile (p_rcoarsegridsolver,rparamList,&
                                    scoarsegridsolver,nlevels,RfilterChain)

      ! Build all levels. Level 1 receives the coarse grid solver.
      do ilev = 1,nlevels

        ! Where we have a coarse grid solver (level 1), we do not need smoothers.
        if (associated(p_rcoarsegridsolver)) then
          nullify(p_rpresmoother)
          nullify(p_rpostsmoother)
        else
          ! Is there a presmoother?
          if (spresmoother .ne. '') then
            call linsolinit_initFromFile (p_rpresmoother,rparamList,&
                                          spresmoother,nlevels,RfilterChain)
          else
            nullify(p_rpresmoother)
          end if

          ! Is there a postsmoother?
          if (spostsmoother .ne. '') then
            ! Check if pre- and postsmoother are identical.
            call sys_toupper (spresmoother)
            call sys_toupper (spostsmoother)
            if (spresmoother .ne. spostsmoother) then
              call linsolinit_initFromFile (p_rpostsmoother,rparamList,&
                                            spostsmoother,nlevels,RfilterChain)
            else
              ! Let the pointer point to the presmoother
              p_rpostsmoother => p_rpresmoother
            end if

          else
            nullify(p_rpostsmoother)
          end if
        end if

        ! Add the level to Multigrid
        call linsol_addMultigridLevel (p_rlevelInfo,p_rsolverNode, &
                    rinterlevelProjection, &
                    p_rpresmoother,p_rpostsmoother,p_rcoarseGridSolver)

        ! Reset the coarse grid solver pointer to NULL.
        ! That way, the coarse grid solver is only attached to level 1.
        nullify(p_rcoarsegridsolver)

      end do

    case (LINSOL_ALG_MULTIGRID2)
      ! Multigrid solver
      !
      ! Parameters:
      !  icycle         = 0               -> F-cycle
      !                 = 1               -> V-cycle
      !                 = 2               -> W-cycle
      !  dalphaMin      >= 0.0            -> minimum damping parameter; standard = 1.0
      !  dalphaMin      >= 0.0            -> maximum damping parameter; standard = 1.0
      !  spreSmootherName                 -> Name of the presmoother section
      !  spostSmootherName                -> Name of the postsmoother section
      !  scoarseGridSolver                -> Name of the coarse grid solver section

      ! At first, initialise the solver:
      call linsol_initMultigrid2(p_rsolverNode, nlevels, p_Rfilter)

      ! Then, get solver specific data.
      call parlst_getvalue_int (p_rsection, 'icycle', &
                                p_rsolverNode%p_rsubnodeMultigrid2%icycle,&
                                p_rsolverNode%p_rsubnodeMultigrid2%icycle)

      ! Coarse grid correction parameters
      call parlst_getvalue_int (p_rsection, 'ccorrectionTypeAlpha', &
           p_rsolverNode%p_rsubnodeMultigrid2%rcoarseGridCorrection%ccorrectionType,&
           p_rsolverNode%p_rsubnodeMultigrid2%rcoarseGridCorrection%ccorrectionType)

      call parlst_getvalue_double (p_rsection, 'dalphaMin', &
           p_rsolverNode%p_rsubnodeMultigrid2%rcoarseGridCorrection%dalphaMin,&
           p_rsolverNode%p_rsubnodeMultigrid2%rcoarseGridCorrection%dalphaMin)

      call parlst_getvalue_double (p_rsection, 'dalphaMax', &
           p_rsolverNode%p_rsubnodeMultigrid2%rcoarseGridCorrection%dalphaMax,&
           p_rsolverNode%p_rsubnodeMultigrid2%rcoarseGridCorrection%dalphaMax)

      ! Do we have a presmoother?
      call parlst_getvalue_string (p_rsection, 'spreSmootherName', sString,'')
      spresmoother = ''
      if (sString .ne. '') read (sString,*) spresmoother

      ! A postsmoother?
      call parlst_getvalue_string (p_rsection, 'spostSmootherName', sString,'')
      spostsmoother = ''
      if (sString .ne. '') read (sString,*) spostsmoother

      ! A coarse grid solver?
      call parlst_getvalue_string (p_rsection, 'scoarseGridSolver', sString,'')
      scoarsegridsolver = ''
      if (sString .ne. '') read (sString,*) scoarsegridsolver

      ! Initialise the coarse grid solver - if there is one. There must be one!
      if (scoarsegridsolver .eq. '') then
        call output_line ('Cannot create linear solver; ' // &
            'no coarse grid solver for MG!', OU_CLASS_ERROR, OU_MODE_STD, &
            'linsolinit_initFromFile')
        call sys_halt()
      end if

      ! set up a coarse grid solver (always gridlevel 1)
      call linsol_getMultigrid2Level(p_rsolverNode, 1, p_rlevelInfo2)

      call linsolinit_initFromFile(p_rlevelInfo2%p_rcoarseGridSolver, rparamList,&
                                   scoarsegridsolver, nlevels, RfilterChain)

      ! Build all the  other levels
      do ilev = 2,nlevels

        ! add this multigrid level
        call linsol_getMultigrid2Level(p_rsolverNode, ilev, p_rlevelInfo2)

        ! Is there a presmoother?
        if (spresmoother .ne. '') then
          call linsolinit_initFromFile(p_rlevelInfo2%p_rpresmoother, rparamList, &
                                       spresmoother, nlevels, RfilterChain)
        else
          nullify(p_rlevelInfo2%p_rpresmoother)
        end if

        ! Is there a postsmoother?
        if (spostsmoother .ne. '') then
          ! Check if pre- and postsmoother are identical.
          call sys_toupper (spresmoother)
          call sys_toupper (spostsmoother)
          if (spresmoother .ne. spostsmoother) then
            ! they are not identical, so read in the post smoother
            call linsolinit_initFromFile(p_rlevelInfo2%p_rpostsmoother, rparamList, &
                                         spostsmoother, nlevels, RfilterChain)
          else
            ! otherwise, let the pointer point to the presmoother
            p_rlevelInfo2%p_rpostsmoother => p_rlevelInfo2%p_rpresmoother
          end if
        else
          nullify(p_rlevelInfo2%p_rpostsmoother)
        end if

      end do

    case default
      call output_line ('Cannot create linear solver; ' // &
            'unsupported solver type isolverType = ' // &
            trim(sys_si(isolverType,8)), &
             OU_CLASS_ERROR, OU_MODE_STD, 'linsolinit_initFromFile')
      call sys_halt()

    end select ! isolverType

    ! Ok, we should have a solver node p_rsolverNode now, initialised with
    ! standard parameters. The next task is to change all parameters
    ! in the solver node to those which appear in the given configuration.
    ! If a parameter does not exist in the given configuration, the default
    ! value (that which is already stored in the solver configuration)
    ! must be used.
    !
    ! So parse the given parameters now to initialise the solver node:

    call linsolinit_initParams (p_rsolverNode,rparamList,ssolverName)

    ! Up to now, we initialised for a linear solver. In case this solver is
    ! used as smoother in a Multigrid algorithm, this initialisation
    ! is not comppletely correct - we have to transform the solver into
    ! one that performs a fixed number of iterations.
    !
    ! We identify a smoother by the existence of the parameter nsmoothingSteps:

    if (parlst_queryvalue (p_rsection, 'nsmoothingSteps') .ne. 0) then

      call parlst_getvalue_int (p_rsection, 'nsmoothingSteps', &
                                i1, p_rsolverNode%nminIterations)

      ! Convert the solver node into a smoother:
      call linsol_convertToSmoother (p_rsolverNode,i1)

    end if

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine linsolinit_initParams (rsolverNode,rparamList,ssection,csolverType)

!<description>
  ! This section reads the standard parameters for a linear solver from the
  ! section ssection in the parameter list rparamList and changes
  ! specified entries in the solver node rsolverNode. Variables whose
  ! values are not specified in the parameter list are left unchanged
  ! in rsolverNode.
  !
  ! The parameters that can be specified in the parameter list have exactly
  ! the same names as in the linear solver structure rsolverNode.
  ! Note that not all parameters in rsolverNode are changed by this routine.
  ! 'Critical' parameters like solver type etc. are left unchanged.
  !
  ! csolverType allows to specify a special solver type whose parameters should
  ! be initialised. If this is desired, the routine should be called twice --
  ! once with csolverType=LINSOL_ALG_UNDEFINED to initialise the main parameters
  ! and once with csolverType=LINSOL_ALG_xxxx to initialise the special,
  ! solver dependent parameters.
!</description>

!<input>
  ! The parameter list that contains the whole solver configuration.
  type(t_parlist), intent(in) :: rparamList

  ! The name of a section in rparamList that contains the configuration of
  ! the linear solver.
  character(LEN=*), intent(in) :: ssection

  ! OPTIONAL: Type of solver structure that should be initialised.
  ! If unspecified or if set to LINSOL_ALG_UNDEFINED, the parameters in the
  ! main solver structure are initialised.
  ! Otherwise, csolverType must be a LINSOL_ALG_xxx constant that specifies
  ! a special solver type whose parameters should be initialised (e.g.
  ! the 'ikrylovDim' parameter of the GMRES method).
  integer, intent(in), optional :: csolverType
!</input>

!<output>
  ! A linear solver node whose parameters should be changed according to the
  ! parameters in rparamList.
  type(t_linsolNode), intent(inout) :: rsolverNode
!</output>

!</subroutine>

    ! local variables
    integer :: csolver
    type(t_parlstSection), pointer :: p_rsection
    integer :: i1,ikrylowDim,nlevels
    character(LEN=SYS_STRLEN) :: sstring
    integer, dimension(3) :: Ikdim
    type(t_linsolDeflGMRESLevelInfo), pointer :: p_rlevelInfoMK

    csolver = LINSOL_ALG_UNDEFINED
    if (present(csolverType)) csolver = csolverType

    ! Check that there is a section called ssolverName - otherwise we
    ! cannot create anything!

    call parlst_querysection(rparamList, ssection, p_rsection)

    if (.not. associated(p_rsection)) then
      call output_line ('Cannot create linear solver; no section ''' &
          // trim(ssection) // '''!', OU_CLASS_ERROR, OU_MODE_STD, &
          'linsolinit_initParams')
      call sys_halt()
    end if

    ! Now comes a biiig select for all the different types of solvers
    ! that are supported by this routine.
    select case (csolver)

    case (LINSOL_ALG_UMFPACK4)
      ! UMFPACK4 solver
      !
      ! Init the debug matrix output.
      ! Parameters:
      !  imatrixDebugOutput = 2   -> Write the matrix in Maple output
      !                     = 1   -> Write the matrix in text output
      !                     = 0   -> Do not write out the matrix
      call parlst_getvalue_int (p_rsection, 'imatrixDebugOutput', i1, -1)
      if (i1 .ne. -1) then
        rsolverNode%p_rsubnodeUMFPACK4%imatrixDebugOutput = i1
      end if

      call parlst_getvalue_string (p_rsection, "smatrixName", sstring, "", bdequote=.true.)
      if (sstring .ne. "") then
        rsolverNode%p_rsubnodeUMFPACK4%smatrixName = sstring
      end if

    case (LINSOL_ALG_SSOR)
      ! SSOR solver
      !
      ! Init the solver node.
      ! Parameters:
      !  iscale = 1   -> Scale preconditioned vector according to actual formula
      !                  in the literature
      !         = 0   -> no scaling, original FEAT implementation
      call parlst_getvalue_int (p_rsection, 'iscale', i1, -1)
      if (i1 .ne. -1) then
        rsolverNode%p_rsubnodeSSOR%bscale = i1 .eq. 1
      end if

    case (LINSOL_ALG_GMRES)
      ! GMRES solver

      ! Krylow space dimension
      call parlst_getvalue_int (p_rsection, 'ikrylovDim', ikrylowDim,40)

      ! Apply Gram Schmidt twice
      call parlst_getvalue_int (p_rsection, 'btwiceGS', i1,-1)

      rsolverNode%p_rsubnodeGMRES%ikrylovDim = ikrylowDim
      if (i1 .ne. -1) then
        rsolverNode%p_rsubnodeGMRES%btwiceGS = i1 .eq. 1
      end if

    case (LINSOL_ALG_MILUS1x1)
      ! (M)ILU solver
      !
      ! Parameters:
      !  ifillinLevel = 0 / 1 / 2 / ...   -> Fill-in level for factorisation
      !  drelax       = 0.0               -> Build ILU(s)
      !               > 0.0               -> Build MILU(s)
      !
      ! Get fill-in level
      call parlst_getvalue_int (p_rsection, 'ifill', &
        rsolverNode%p_rsubnodeMILUs1x1%ifill, &
        rsolverNode%p_rsubnodeMILUs1x1%ifill)

      ! Get MILU relaxsation parameter
      call parlst_getvalue_double (p_rsection, 'drelax', &
        rsolverNode%p_rsubnodeMILUs1x1%drelax, &
        rsolverNode%p_rsubnodeMILUs1x1%drelax)

    case (LINSOL_ALG_MULTIGRID)
      ! Multigrid solver
      !
      ! Parameters:
      !  icycle         = 0               -> F-cycle
      !                 = 1               -> V-cycle
      !                 = 2               -> W-cycle
      !  dalphaMin      >= 0.0            -> minimum damping parameter; standard = 1.0
      !  dalphaMin      >= 0.0            -> maximum damping parameter; standard = 1.0

      ! Then, get solver specific data.
      call parlst_getvalue_int (p_rsection, 'icycle', &
                                rsolverNode%p_rsubnodeMultigrid%icycle,&
                                rsolverNode%p_rsubnodeMultigrid%icycle)

      ! Coarse grid correction parameters
      call parlst_getvalue_int (p_rsection, 'ccorrectionTypeAlpha', &
           rsolverNode%p_rsubnodeMultigrid%rcoarseGridCorrection%ccorrectionType,&
           rsolverNode%p_rsubnodeMultigrid%rcoarseGridCorrection%ccorrectionType)

      call parlst_getvalue_double (p_rsection, 'dalphaMin', &
           rsolverNode%p_rsubnodeMultigrid%rcoarseGridCorrection%dalphaMin,&
           rsolverNode%p_rsubnodeMultigrid%rcoarseGridCorrection%dalphaMin)

      call parlst_getvalue_double (p_rsection, 'dalphaMax', &
           rsolverNode%p_rsubnodeMultigrid%rcoarseGridCorrection%dalphaMax,&
           rsolverNode%p_rsubnodeMultigrid%rcoarseGridCorrection%dalphaMax)

    case (LINSOL_ALG_MULTIGRID2)
      ! Multigrid solver
      !
      ! Parameters:
      !  icycle         = 0               -> F-cycle
      !                 = 1               -> V-cycle
      !                 = 2               -> W-cycle
      !  dalphaMin      >= 0.0            -> minimum damping parameter; standard = 1.0
      !  dalphaMin      >= 0.0            -> maximum damping parameter; standard = 1.0

      ! Then, get solver specific data.
      call parlst_getvalue_int (p_rsection, 'icycle', &
                                rsolverNode%p_rsubnodeMultigrid2%icycle,&
                                rsolverNode%p_rsubnodeMultigrid2%icycle)

      ! Coarse grid correction parameters
      call parlst_getvalue_int (p_rsection, 'ccorrectionTypeAlpha', &
           rsolverNode%p_rsubnodeMultigrid2%rcoarseGridCorrection%ccorrectionType,&
           rsolverNode%p_rsubnodeMultigrid2%rcoarseGridCorrection%ccorrectionType)

      call parlst_getvalue_double (p_rsection, 'dalphaMin', &
           rsolverNode%p_rsubnodeMultigrid2%rcoarseGridCorrection%dalphaMin,&
           rsolverNode%p_rsubnodeMultigrid2%rcoarseGridCorrection%dalphaMin)

      call parlst_getvalue_double (p_rsection, 'dalphaMax', &
           rsolverNode%p_rsubnodeMultigrid2%rcoarseGridCorrection%dalphaMax,&
           rsolverNode%p_rsubnodeMultigrid2%rcoarseGridCorrection%dalphaMax)

    case (LINSOL_ALG_DEFLGMRES)

      ! Deflated GMRES method.

      ! Get the relaxation parameter
      call parlst_getvalue_double (p_rsection, 'dkrylowRelax', &
           rsolverNode%p_rsubnodeDeflGMRES%drelax,&
           rsolverNode%p_rsubnodeDeflGMRES%drelax)

      ! Get the Krylow space dimension for the different levels.
      call parlst_getvalue_string (p_rsection, "NkrylowDimension", sstring, "4 2 2")
      read (sstring,*) Ikdim(1), Ikdim(2), Ikdim(3)

      ! Method how to determine the maximum eigenvalue
      call parlst_getvalue_int (p_rsection, 'cmaxEigenvalMethod', &
           rsolverNode%p_rsubnodeDeflGMRES%cmaxEigenvalMethod,&
           rsolverNode%p_rsubnodeDeflGMRES%cmaxEigenvalMethod)

      ! Impose the number of preconditioning steps
      nlevels = linsol_getDeflGMRESLevelCount(rsolverNode)

      do i1 = 1,nlevels
        ! Get the level
        call linsol_getDeflGMRESLevel(rsolverNode,i1,p_rlevelInfoMK)

        ! Dimension of the Krylow space
        p_rlevelInfoMK%ikrylowDim = Ikdim(min(3,i1))
      end do

    case default

      ! Initialise the main solver parameters

      call parlst_getvalue_double (p_rsection, 'domega', &
                                  rsolverNode%domega, rsolverNode%domega)

      call parlst_getvalue_int (p_rsection, 'nminIterations', &
                                rsolverNode%nminIterations, rsolverNode%nminIterations)

      call parlst_getvalue_int (p_rsection, 'nmaxIterations', &
                                rsolverNode%nmaxIterations, rsolverNode%nmaxIterations)

      call parlst_getvalue_string (p_rsection, 'iresCheck', sstring, 'YES')
      call sys_toUpper(sstring)
      if (sstring .eq. 'NO') then
        ! if the string is exactly 'NO' then disable checking of residual norm
        rsolverNode%iresCheck = NO
      else
        ! in all other cases enable checking of residual norm
        rsolverNode%iresCheck = YES
      endif

      call parlst_getvalue_double (p_rsection, 'depsRel', &
                                  rsolverNode%depsRel, rsolverNode%depsRel)

      call parlst_getvalue_double (p_rsection, 'depsAbs', &
                                  rsolverNode%depsAbs, rsolverNode%depsAbs)

      call parlst_getvalue_double (p_rsection, 'depsDiff', &
                                  rsolverNode%depsDiff, rsolverNode%depsDiff)

      call parlst_getvalue_int (p_rsection, 'iresNorm', &
                                rsolverNode%iresNorm, rsolverNode%iresNorm)

      call parlst_getvalue_int (p_rsection, 'niteAsymptoticCVR', &
                                rsolverNode%niteAsymptoticCVR, &
                                rsolverNode%niteAsymptoticCVR)

      call parlst_getvalue_int (p_rsection, 'istoppingCriterion', &
                                rsolverNode%istoppingCriterion, &
                                rsolverNode%istoppingCriterion)

      call parlst_getvalue_int (p_rsection, 'ioutputLevel', &
                                rsolverNode%ioutputLevel, rsolverNode%ioutputLevel)

      call parlst_getvalue_int (p_rsection, 'niteResOutput', &
                                rsolverNode%niteResOutput, rsolverNode%niteResOutput)

      call parlst_getvalue_int (p_rsection, 'isolverSubgroup', &
                                rsolverNode%isolverSubgroup, &
                                rsolverNode%isolverSubgroup)

    end select ! isolverType

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
! To compile

! ifort -check bounds -warn all -o mpas_to_latlon src/mpas_to_latlon.f90
! src/mpas_filter_cells_by_area.f src/flip_to_cf.f90 src/mpas_vort_cell.f
! src/datetime_module.f90 -lnetcdf

! Original code from Michael Duda duda@ucar.edu
! adapted by Dave Ahijevych ahijevyc@ucar.edu

!
! This code interpolates fields on an MPAS mesh to a lat-lon grid.
! For fields defined on cells we find use the nearest cell and smooth by X km.
!
! Reads list of fields from standard input.
!

! Note from Michael:
! For fields defined on vertexes (e.g. vorticity), using the lat-lon edge point,
! we seek the nearest edge point and use that edge and the 4 edges
! on that edge to interpolate the velocities.
! The interpolation uses both the nearest edge and neighboring edges, and uses both
! the normal and reconstructed-tangential velocities

! Note from ahijevych:
! For a while I had problems with ifort. nlat would change
! from 180 to 14 when I inquired the dimension id of a variable.  Strange.
! Solved when I decreased the size of the dimension for var_output_id and
! var_input_id from 100 to 40, or lower. Started changing allocatable arrays to
! static arrays because I thought that was the problem.  Recompiled in pgf90 and avoided the
! problem.

program mpas_to_latlon
    use netcdf
    implicit none

    integer, parameter :: RKIND = selected_real_kind(12)
    integer, parameter :: max_lon_dimsize = 7200
    integer, parameter :: max_lat_dimsize = 3600
    integer, parameter :: max_nCells_dimsize = 36864002

    real (kind=RKIND) :: pii

    include 'netcdf.inc'

    integer, external :: iargc

    real (kind=RKIND), dimension(max_nCells_dimsize) :: lat_cell, lon_cell
    real (kind=RKIND), allocatable, dimension(:) :: lat_vertex, lon_vertex, areaCell
    real (kind=RKIND), allocatable, dimension(:,:) :: input_var, input_var_vertices
    real (kind=RKIND), allocatable, dimension(:,:,:) :: output_var, output_var_cf

    integer, allocatable, dimension(:,:) :: verticesOnCell, cellsOnVertex, cellsOnCell
    integer, allocatable, dimension(:,:) :: nearest_vertex_ll
    integer, dimension(max_nCells_dimsize) :: nedgesOnCell
    ! used to allocate(nedgesoncell(nCells))  ! BUG FIX (used to be nEdges)

    character (len=8)  :: now_date
    character (len=10) :: now_time
    character (len=5)  :: now_zone
    character (len=30) :: junkc
    character (len=30), dimension(100) :: interp_var_input, interp_var_output
    integer, dimension(100) :: interp_var_nz, var_output_id, var_input_id
    real (kind=RKIND) :: grid_spacing ! lat-lon mesh spacing in degrees
    real :: filter_radius_km  ! smoothing radius km (zero for no smoothing)
    integer :: nvars, ivar, ierr, ncid, ncells, ncells_id
    integer :: latcell_id, loncell_id, areaCell_id
    integer :: nvertices, nvertices_id, levels_id, nvertlevels
    integer :: verticesOnCell_id, latvertex_id, lonvertex_id
    integer :: time_id, time_axis_id, ntimes
    integer :: i, j, cellsOnVertex_id, maxedges, maxedges_id
    integer :: cellsOnCell_id, nedgesOnCell_id, ndims, natts, xtype
    ! dimids has to be allocated with enough space for at least *ndimsp integers
    ! to be returned - Dave Feb 17, 2013
    integer, dimension(NF_MAX_VAR_DIMS) :: dimids
    integer, dimension(4) :: ncount
    integer :: ncid_ll, lat_ll_id, lon_ll_id, grid_spacing_id, filter_radius_km_id
    integer :: lat_ll_v_id, lon_ll_v_id

    real (kind=RKIND), dimension(max_lat_dimsize,max_lon_dimsize) :: lat_ll, lon_ll
    real (kind=RKIND), dimension(3,max_lat_dimsize,max_lon_dimsize) :: interp_weights
    integer, dimension(3,max_lat_dimsize,max_lon_dimsize) :: interp_cells
    integer :: nlat, nlon, start_vertex, idimn, iatt
    real (kind=RKIND) :: latitude, longitude, startLon, degreesToRadians, RadiansToDegrees, lat0, lat1
    character (len=256) :: fname, arg, savfile, outfile, tmpdir
    character (len=256) :: meshid, attname
    logical :: file_exists

    pii = 2.0_RKIND * asin(1.0_RKIND)
    degreesToRadians = pii / 180.0_RKIND
    RadiansToDegrees = 180.0_RKIND / pii

    if (iargc() < 3) then
       write(0,*)
       write(0,*) 'Usage:'
       write(0,*) 'mpas_to_latlon fname outfile grid_spacing filter_radius_km meshid [lat0 [lat1 [startLon]]]'
       write(0,*)
       write(0,*) ' fname            : input MPAS filename'
       write(0,*) ' outfile          : output filename'
       write(0,*) ' grid_spacing     : output grid spacing in degrees'
       write(0,*) ' filter_radius_km : radius of circular smoothing filter in km'
       write(0,*) ' meshid           : string to identify MPAS mesh (uni, wp, us, etc.)'
       write(0,*) ' lat0             : minimum latitude in deg (default -90)'
       write(0,*) ' lat1             : maximum latitude in deg (default +90)'
       write(0,*) ' startLon         : starting longitude for output grid (default -180)'
       write(0,*)
       write(0,*) 'Example:'
       write(0,*) 'mpas_to_latlon diagnostics.2013-09-01_00:00:00.nc latlon_0.500deg_025km/diagnostics.2013-09-01_00:00:00.nc 0.5 25 mpas3 -5. 50. -180.'
       write(0,*)
       call exit(1)
    end if
    call getarg(1,fname)
    call getarg(2,outfile)

    if (fname == outfile) then
        write(0,*)' mpas_to_latlon: infile equals outfile. Exiting'
        call exit(1)
    end if
    call getarg(3,arg)
    read(arg,*) grid_spacing
    call getarg(4,arg)
    read(arg,*) filter_radius_km
    call getarg(5,arg)
    read(arg,*) meshid
    lat0 = -90.0
    lat1 =  90.0
    startLon = -180.
    call getarg(6,arg)
    if (len_trim(arg).ne.0) read(arg,*) lat0
    call getarg(7,arg)
    if (len_trim(arg).ne.0) read(arg,*) lat1
    call getarg(8,arg)
    if (len_trim(arg).ne.0) read(arg,*) startLon

    !
    !  create lat lon grid dimensions
    !

    nlat = ceiling((lat1-lat0) / grid_spacing)
    nlon = ceiling(360. / grid_spacing)
    if(nlat>max_lat_dimsize.or.nlon>max_lon_dimsize)then
       write(0,*)
       write(0,*)' mpas_to_latlon: either max_lat_dimsize or max_lon_dimsize is too small. &
           recompile after making it bigger.'
       write(0,*)
       call exit(1)
    endif

    !
    ! Populate list of variables to interpolate
    !

    ! read list of 2D vars from standard input
    write(0,*) ' mpas_to_latlon: reading from standard input (use ctrl-d to end manual entry) '
    read(5,*,iostat=ierr) junkc
    ivar = 1
    do while(ierr == 0)
       interp_var_input(ivar) = junkc
       interp_var_output(ivar) = interp_var_input(ivar)
       ivar = ivar+1
       read(5,*,iostat=ierr) junkc
    end do
    nvars = ivar-1

    !
    !  open MPAS file and read description of MPAS mesh
    !

    ierr = nf_open(trim(fname), NF_NOWRITE, ncid)
    write(0,*)'opened '//trim(fname)//' ncid=',ncid
    if (ierr.ne.NF_NOERR) then
        write(0,*) ' nf_open error opening mpas file '//trim(fname)
        call handle_err(ierr)
    end if

    !
    ! Input must have dimensions 'nCells', 'nVertices', and 'Time'.
    !

    ierr = nf_inq_dimid(ncid, 'nCells', ncells_id)
    if(ierr.ne.NF_NOERR)then
        write(0,*) ' nf_inq_dimid err for ncells'
        call handle_err(ierr)
    end if

    ierr = nf_inq_dim(ncid,ncells_id,junkc,ncells)
    if(ierr.ne.NF_NOERR)then
        write(0,*) ' nf_inq_dim err for ncells'
        write(0,*) ' ncells for mesh ',ncells
        call handle_err(ierr)
    end if
    if(ncells.gt.max_nCells_dimsize)then
        write(0,*) 'ncells ',ncells,' exceeds limit ',max_nCells_dimsize
        write(0,*) 'increase max_nCells_dimsize parameter and recompile'
        stop
    end if

    ierr = nf_inq_dimid(ncid,'nVertices',nVertices_id)
    if(ierr.ne.NF_NOERR)then
        write(0,*) ' nf_inq_dimid err for nVertices id'
        call handle_err(ierr)
    end if

    ierr = nf_inq_dim(ncid,nVertices_id,junkc,nVertices)
    if(ierr.ne.NF_NOERR)then
        write(0,*) ' nf_inq_dim err for nVertices'
        write(0,*) ' nVertices for mesh ',nVertices
        call handle_err(ierr)
    end if

    ierr = nf_inq_dimid(ncid, 'Time', time_id)
    if(ierr.ne.NF_NOERR)then
        write(0,*) ' nf_inq_dimid err for Time id'
        call handle_err(ierr)
    end if

    ierr = nf_inq_dim(ncid, time_id, junkc, ntimes)
    if(ierr.ne.NF_NOERR.or.ntimes.ne.1)then
        write(0,*) ' err for times'
        write(0,*) ' number of times ', ntimes
        call handle_err(ierr)
    end if

    !
    ! create output file
    !
    ierr = nf_create(outfile,NF_64BIT_OFFSET, ncid_ll)
    if (ierr.ne.NF_NOERR) then
        write(0,*) ' err for nf_create '//trim(outfile)
        write(0,*) ' do you need to create a subdirectory first?'
        call handle_err(ierr)
    end if

    ierr = nf_def_dim( ncid_ll, 'time', NF_UNLIMITED, time_id )
    if(ierr.ne.NF_NOERR)write(0,*) ' err on time_id dim def'
    call handle_err(ierr)

    ierr = nf_def_dim( ncid_ll, 'lat', nlat, lat_ll_id )
    if (ierr.ne.NF_NOERR) write(0,*) ' err defining lat dim'
    call handle_err(ierr)

    ierr = nf_def_dim( ncid_ll, 'lon', nlon, lon_ll_id )
    if (ierr.ne.NF_NOERR) write(0,*) ' err defining lon dim'
    call handle_err(ierr)

    ncount(1) = time_id
    ncount(2) = 0
    ncount(3) = 0
    ncount(4) = 0
    ierr = nf_def_var( ncid_ll, 'time', NF_DOUBLE, 1, ncount, time_axis_id )
    call handle_err(ierr)
    ierr = nf_put_att_text(ncid_ll, time_axis_id, 'long_name', len('time'),      'time')
    call handle_err(ierr)
    write(junkc,'(a11,a)') 'days since ', trim('1970-01-01 00:00:00')
    ierr = nf_put_att_text(ncid_ll, time_axis_id, 'units',     len_trim(junkc),  trim(junkc))
    call handle_err(ierr)
    ierr = nf_put_att_text(ncid_ll, time_axis_id, 'calendar',  len('gregorian'), 'gregorian')
    call handle_err(ierr)

    ncount(1) = lat_ll_id
    ncount(2) = 0
    ncount(3) = 0
    ncount(4) = 0
    ierr = nf_def_var( ncid_ll, 'lat', NF_DOUBLE, 1, ncount, lat_ll_v_id )
    if (ierr.ne.NF_NOERR) write(0,*) ' err on lat_ll_id dim def'
    call handle_err(ierr)
    ierr = nf_put_att_text(ncid_ll, lat_ll_v_id, 'long_name', len('latitude'),      'latitude')
    call handle_err(ierr)
    ierr = nf_put_att_text(ncid_ll, lat_ll_v_id, 'units',     len('degrees_north'), 'degrees_north')
    call handle_err(ierr)
    ierr = nf_put_att(ncid_ll, lat_ll_v_id, 'valid_min', NF_DOUBLE, 1, lat0)
    call handle_err(ierr)
    ierr = nf_put_att(ncid_ll, lat_ll_v_id, 'valid_max', NF_DOUBLE, 1, lat1)
    call handle_err(ierr)

    ncount(1) = lon_ll_id
    ncount(2) = 0
    ncount(3) = 0
    ncount(4) = 0
    ierr = nf_def_var( ncid_ll, 'lon', NF_DOUBLE, 1, ncount, lon_ll_v_id )
    if (ierr.ne.NF_NOERR) write(0,*) ' err on lon_ll_id dim def'
    call handle_err(ierr)
    ierr = nf_put_att_text(ncid_ll, lon_ll_v_id, 'long_name', len('longitude'),    'longitude')
    call handle_err(ierr)
    ierr = nf_put_att_text(ncid_ll, lon_ll_v_id, 'units',     len('degrees_east'), 'degrees_east')
    call handle_err(ierr)

    ierr = nf_def_var( ncid_ll, 'mesh_spacing', NF_DOUBLE, 0, 0, grid_spacing_id )
    call handle_err(ierr)
    ierr = nf_put_att_text(ncid_ll, grid_spacing_id, 'long_name', len('lat lon grid spacing (command line argument passed to mpas_to_latlon)'), &
                                                                      'lat lon grid spacing (command line argument passed to mpas_to_latlon)')
    call handle_err(ierr)
    ierr = nf_put_att_text(ncid_ll, grid_spacing_id, 'units',     len('degrees'),    'degrees')
    call handle_err(ierr)
    ierr = nf_def_var( ncid_ll, 'filter_radius_km', NF_DOUBLE,    0, 0, filter_radius_km_id )
    call handle_err(ierr)
    ierr = nf_put_att_text(ncid_ll, filter_radius_km_id, 'long_name', len('radius of influence for mpas_to_latlon smoothing filter'), &
                                                                     'radius of influence for mpas_to_latlon smoothing filter')
    call handle_err(ierr)
    ierr = nf_put_att_text(ncid_ll, filter_radius_km_id, 'units',     len('km'),    'km')
    call handle_err(ierr)

    do i=1,nvars
        ierr = nf_inq_varid(ncid,interp_var_input(i),var_input_id(i))
        call handle_err(ierr)
        write(0,*) 'got var_input_id for '//TRIM(interp_var_input(i)), var_input_id(i)
        if (ierr.ne.NF_NOERR) then
            write(0,*) ' ierr ',ierr,' for input var ',interp_var_input(i),var_input_id(i)
            call handle_err(ierr)
        end if
        ierr = nf_inq_var(ncid, var_input_id(i), junkc, xtype, ndims, dimids, natts)
        if (ierr.ne.NF_NOERR) then
            write(0,*) ' ierr ',ierr,' for nf_inq_var ',interp_var_input(i),var_input_id(i)
            call handle_err(ierr)
        end if

        ! TODO: pass through 1-D coordinate variables without interpolating
        ncount = -1
        ncount(1) = lon_ll_id
        ncount(2) = lat_ll_id
        ! You add 1 dimension by converting nCells or nVertices to (lat, lon).
        ncount(ndims+1) = time_id
        interp_var_nz(i) = 1
        ! See if any dimensions are vertical. If so, define dimension in output file.

        do idimn=1,ndims
           ierr = nf_inq_dim(ncid,dimids(idimn),junkc,nVertLevels)
           if (ierr.ne.NF_NOERR) write(0,*) ' err inq dimn ',idimn,junkc
           call handle_err(ierr)
           if (trim(junkc) .eq. 'nVertLevels' .or. junkc(1:10) .eq. 'nIsoLevels' .or. junkc(2:12) == '_iso_levels') then
              write(0,*) ' found input vertical dimension '//TRIM(junkc)
              if (ndims < 3) then
                 write(0,*) '  but ndims=',ndims
                 stop ! sanity check
              end if
              ierr = nf_inq_dimid(ncid_ll, junkc, levels_id)
              if(ierr.ne.NF_NOERR) then
                ierr = nf_def_dim(ncid_ll, junkc, nVertLevels, levels_id)
                write(0,*) ' define output vertical dimension '//TRIM(junkc)
                if (ierr.ne.NF_NOERR) write(0,*) ' err defining vertical dimension for output',idimn,junkc
                call handle_err(ierr)
              end if
              ncount(3) = levels_id
              interp_var_nz(i) = nVertLevels
           else
              write(0,*) ' input dimension '//TRIM(junkc)//' not vertical'
           end if
        end do
        if (ndims >= 3 .and. ncount(3) == -1) then
           write(0,'(A,I1,A)') ' no vertical dimension found in ',ndims,'-D input. stopping.'
           stop ! sanity check
        end if
        ! Definine output variable. Dimensions are longitude, latitude, and
        ! time. If vertical dimension present, dim 3 is vertical and dim 4
        ! is time.
        !write(0,*)'defining '//interp_var_output(i)
        !write(0,*)'ndims+1,ncount= ',ndims+1,ncount
        ! Used to force NF_DOUBLE - changed to xtype (the type of the input
        ! variable)
        ierr = nf_def_var( ncid_ll, interp_var_output(i), xtype, ndims+1, ncount, var_output_id(i) )
        if(ierr.ne.NF_NOERR) write(0,*) ' output var ',interp_var_output(i),' id ',var_output_id(i)
        call handle_err(ierr)

        ! Copy attributes from input variable to output variable
        do iatt=1,natts
            ierr = nf_inq_attname(ncid, var_input_id(i), iatt, attname)
            if (ierr.ne.NF_NOERR) then
                write(0,*) ' err for nf_inq_attname '//interp_var_input(i),iatt
                call handle_err(ierr)
            end if
            ierr = nf_copy_att(ncid, var_input_id(i), attname, ncid_ll, var_output_id(i))
            if (ierr.ne.NF_NOERR) then
                write(0,*) ' err for nf_copy_att ',attname
                call handle_err(ierr)
            end if
        end do
    end do

    ierr = nf_put_att_text(ncid_ll, NF_GLOBAL, 'Conventions', len('CF-1.0'), 'CF-1.0')
    call handle_err(ierr)
    ierr = nf_put_att_text(ncid_ll, NF_GLOBAL, 'model', len('mpas'), 'mpas')
    call handle_err(ierr)
    ierr = nf_put_att_text(ncid_ll, NF_GLOBAL, 'grid', len('interp_latlon'), 'interp_latlon')
    call handle_err(ierr)
    ierr = nf_put_att_text(ncid_ll, NF_GLOBAL, 'native_grid', len('voronoi'), 'voronoi')
    call handle_err(ierr)
    ierr = nf_put_att_text(ncid_ll, NF_GLOBAL, 'original_filename', len(trim(fname)), TRIM(fname))
    call handle_err(ierr)
    call date_and_time(DATE=now_date, TIME = now_time, ZONE = now_zone)
    ierr = nf_put_att_text(ncid_ll, NF_GLOBAL, 'run_date_and_time', len('mpas_to_latlon: '//now_date//' '//now_time//' '//now_zone), &
                                                                        'mpas_to_latlon: '//now_date//' '//now_time//' '//now_zone)
    call handle_err(ierr)

    ! Copy global attributes to output file
    ierr = nf90_inquire(ncid, nAttributes=natts)
    if (ierr.ne.NF_NOERR) call handle_err(ierr)
    do iatt=1,natts
        ierr = nf_inq_attname(ncid, NF_GLOBAL, iatt, attname)
        if (ierr.ne.NF_NOERR) then
            write(0,*) ' err for nf_inq_attname global att ',iatt
            call handle_err(ierr)
        end if
        ierr = nf_copy_att(ncid, NF_GLOBAL, attname, ncid_ll, NF_GLOBAL)
        if (ierr.ne.NF_NOERR) then
            write(0,*) ' err for nf_copy_att ',attname
            call handle_err(ierr)
        end if
    end do

    ! netcdf - exiting define mode for output file
    ierr = nf_enddef(ncid_ll)
    if(ierr.ne.NF_NOERR)write(0,*) ' err for nf_enddef(ncid_ll)'
    call handle_err(ierr)

    !
    !  create lat lon mesh for interp
    !
    do j=1, nlon
        longitude = startLon + 0.5*grid_spacing + float(j-1)*grid_spacing
        do i=1, nlat
            latitude = (lat0+0.5*grid_spacing) + float(i-1)*grid_spacing
            lat_ll(i,j) = latitude*degreesToRadians
            lon_ll(i,j) = longitude*degreesToRadians
        end do
    end do

    ! get mpas->latlon geometry
    ! use a save file if it exists
    ! Assumed triangulation may be uniquely identified by nCells and lat range of output.
    ! This assumption burned me when I tried the mpas3 mesh, which is simply the
    ! mpas2 mesh shifted around so the fine part is over the W Pacific.
    ! Add a mesh ID (mpas, mpas2, or mpas3)
    ! Finally added a mesh ID and recompiled on Jun 5, 2014.
    call get_environment_variable("TMPDIR", tmpdir)
    ! Assumes $TMPDIR
    ! available to store savefile. Perhaps replace with home directory or /tmp
    ! directory for other institutions to use.
    ! Tried to zero-pad lat and lon so we don't get spaces in the filename, but
    ! fortran only zero-pads integers.
    ! SP turns on leading + or - for the rest of the format string.
    write(fmt='(A,"/",A,"_",I8.8,"_",F5.3,"deg",SP,F7.3,"N",F7.3,"N",F8.3,"E")', unit=savfile) &
      TRIM(tmpdir),trim(meshid), nCells, grid_spacing, lat0, lat1, startLon
    write(0,'(A)')'looking for save file "'//trim(savfile)//'"'
    inquire(file=savfile, exist=file_exists)
    if (file_exists) then
        write(0,*) ' mpas_to_latlon: using save file '//trim(savfile)
        open(unit=2, file=savfile, form='unformatted')
        read(2) maxEdges, interp_weights, interp_cells
        !write(0,*) maxEdges, nCells, nVertices
        !write(0,*) 'interp_weights=',interp_weights, 'interp_cells=',interp_cells
        read(2) nEdgesOnCell
        allocate(verticesOnCell(maxEdges,nCells))
        read(2) verticesOnCell
        allocate(cellsOnCell(maxEdges,nCells))
        read(2) cellsOnCell
        allocate(areaCell(nCells))
        read(2) areaCell
        close(2)
    else

        write(0,'(A)')'did not find save file'
        ierr = nf_inq_dimid(ncid, 'maxEdges', maxedges_id)
        if (ierr.ne.NF_NOERR) then
            write(0,*) ' err for maxedges id ', NF_STRERROR(ierr)
            write(0,*) ' process init.nc first, so maxEdges and nEdges are defined in save file'
            call exit(2)
        end if

        ierr = nf_inq_dim(ncid, maxedges_id, junkc, maxedges)
        if (ierr.ne.NF_NOERR) then
            write(0,*) ' err for nf_inq_dim maxEdges ', NF_STRERROR(ierr)
            call handle_err(ierr)
            call exit(2)
        end if
        write(0,*) ' maxedges for mesh ', maxedges

        allocate(lat_vertex(nvertices))
        allocate(lon_vertex(nvertices))

        !allocate(nedgesoncell(nCells))  ! BUG FIX (used to be nEdges)
        allocate(verticesOnCell(maxedges,ncells))
        allocate(cellsOnCell(maxedges,ncells))
        allocate(cellsOnVertex(3,nVertices))
        allocate(areaCell(nCells))

        ierr = nf_inq_varid(ncid,'latCell',latcell_id)
        call handle_err(ierr)
        ierr = nf_inq_varid(ncid,'lonCell',loncell_id)
        call handle_err(ierr)
        ierr = nf_get_var_double(ncid,latcell_id,lat_cell)
        call handle_err(ierr)
        ierr = nf_get_var_double(ncid,loncell_id,lon_cell)
        call handle_err(ierr)

        ierr = nf_inq_varid(ncid,'latVertex',latVertex_id)
        call handle_err(ierr)
        ierr = nf_inq_varid(ncid,'lonVertex',lonVertex_id)
        call handle_err(ierr)
        ierr = nf_get_var_double(ncid,latVertex_id,lat_Vertex)
        call handle_err(ierr)
        ierr = nf_get_var_double(ncid,lonVertex_id,lon_Vertex)
        call handle_err(ierr)

        ierr = nf_inq_varid(ncid,'verticesOnCell',verticesOnCell_id)
        call handle_err(ierr)
        ierr = nf_get_var_int(ncid,verticesOnCell_id,verticesOnCell)
        call handle_err(ierr)

        ierr = nf_inq_varid(ncid,'cellsOnVertex',cellsOnVertex_id)
        call handle_err(ierr)
        ierr = nf_get_var_int(ncid,cellsOnVertex_id,cellsOnVertex)
        call handle_err(ierr)

        ierr = nf_inq_varid(ncid,'nEdgesOnCell',nedgesOnCell_id)
        call handle_err(ierr)
        ierr = nf_get_var_int(ncid,nedgesoncell_id,nEdgesOnCell)
        call handle_err(ierr)

        ierr = nf_inq_varid(ncid,'cellsOnCell',cellsOnCell_id)
        call handle_err(ierr)
        ierr = nf_get_var_int(ncid,cellsoncell_id,cellsOnCell)
        call handle_err(ierr)

        ierr = nf_inq_varid(ncid,'areaCell',areaCell_id)
        call handle_err(ierr)
        ierr = nf_get_var_double(ncid,areaCell_id,areaCell)
        call handle_err(ierr)

        !
        ! find nearest vertex for interpolation
        !
        allocate(nearest_vertex_ll(nlat,nlon))
        start_vertex = 1
        do j=1,nlon
            do i=1,nlat
                !write(0,*) 'find nearest vertex to',j,i
                nearest_vertex_ll(i,j) = nearest_vertex( lat_ll(i,j), lon_ll(i,j), start_vertex, ncells, nvertices, maxedges, &
                                                         nedgesOnCell, verticesOnCell, cellsOnVertex, &
                                                         lat_Cell, lon_Cell, lat_vertex, lon_vertex )
                !if (nearest_vertex_ll(i,j) .eq. 0) write(0,*), 'after nearest_vertex() the nearest vertex=0 when j=',j,' i=',i
                start_vertex = nearest_vertex_ll(i,j)
            end do
            start_vertex = nearest_vertex_ll(1,j)
        end do

        !
        !  compute interp weights for all the points on the lat-lon mesh
        !

        write(0,*) ' computing interp weights '
        call compute_interp_weights_all ( lat_ll, lon_ll, nearest_vertex_ll, nCells,           &
                                          nVertices, cellsOnVertex, lat_cell, lon_cell,        &
                                          interp_weights, interp_cells, nlat, nlon )

        write(0,*) ' interp weights complete '

        ! write expensive variables to save file, which can be used instead of
        ! rerunning this section of code.
        write(0,'(A)')'writing to save file'
        open(unit=2, file=savfile, form='unformatted')
        write(2) maxEdges, interp_weights, interp_cells
        write(2) nEdgesOnCell
        write(2) verticesOnCell
        write(2) cellsOnCell
        write(2) areaCell
        close(2)
        write(0,'(A)')'wrote expensive variables to save file "'//trim(savfile)//'"'
    end if

    ierr = nf_put_vara_double(ncid_ll, time_axis_id, 1, 1, days_since_1970(trim(fname)) )
    call handle_err(ierr)
    ierr = nf_put_var_double(ncid_ll, lat_ll_v_id, lat_ll(:,1)*RadiansToDegrees )
    call handle_err(ierr)
    ierr = nf_put_var_double(ncid_ll, lon_ll_v_id, lon_ll(1,:)*RadiansToDegrees )
    call handle_err(ierr)
    ierr = nf_put_var_double(ncid_ll, grid_spacing_id, grid_spacing)
    call handle_err(ierr)
    ierr = nf_put_var_double(ncid_ll, filter_radius_km_id, filter_radius_km)
    call handle_err(ierr)

    ! for each variable - read, interp, and write
    do i=1,nvars
        nVertLevels = interp_var_nz(i)
        allocate(input_var(nVertLevels,nCells))
        allocate(output_var(nVertLevels,nlat,nlon))
        allocate(output_var_cf(nlon,nlat,nVertLevels))

        ierr = nf_inq_var(ncid, var_input_id(i), junkc, xtype, ndims, dimids, natts)
        call handle_err(ierr)

        !write(0,*) '1 nlon=',nlon,' nlat=',nlat,' ncid=',ncid,' var_input_id(i)=',var_input_id(i)
        ierr = nf_inq_vardimid(ncid, var_input_id(i), dimids)
        call handle_err(ierr)
        !write(0,*) '2 dimids=',dimids, ' nVertices_id=',nVertices_id,' nlon=',nlon,' nlat=',nlat, ierr
        if (dimids(1) .eq. nVertices_id) then
            write(0,*) trim(interp_var_output(i))//' is a vertex variable with ',nVertices,'vertices'
            allocate(input_var_vertices(nVertLevels,nVertices))
            ierr = nf_get_var_double( ncid,var_input_id(i),input_var_vertices)
            if(ierr.ne.NF_NOERR) then
               write(0,*)'could not read '//trim(interp_var_output(i))//' ierr=',ierr
               write(0,*)'ncid=',ncid,' var_input_id(i)=',var_input_id(i)
               write(0,*)NF_STRERROR(ierr)
               call exit(ierr)
            end if

            ! INPUT: input_var_vertices OUTPUT: input_var (on cells)
            !write(0,*)'calling mpas_vort_cell1'
            call mpas_vort_cell1(nEdgesOnCell, verticesOnCell, maxEdges, nVertLevels, nCells, nVertices, input_var_vertices, input_var)

            deallocate(input_var_vertices)
        else
            ierr = nf_get_var_double( ncid,var_input_id(i),input_var)
            call handle_err(ierr)
        end if

        ! Filter is based on mpas_filter_cells.f, but the number of
        ! smoothing passes is inversely proportional to sqrt(areaCell)
        ! filter_radius_km is in km.
        !write(0,*)'calling mpas_filter_cells_by_area. nEdgesOnCell(1)=',nEdgesOnCell(1), 'cellsOnCell(:,1)=',cellsOnCell(:,1),'maxEdges=',maxEdges,'nCells=',nCells
        !write(0,*)' areaCell(1)=',areaCell(1), ' filter_radius_km=', filter_radius_km
        call mpas_filter_cells_by_area(nEdgesOnCell, cellsOnCell, areaCell, maxEdges, nVertLevels, nCells, input_var, filter_radius_km)
        !write(0,*)'after mpas_filter_cells_by_area: minval(cellsOnCell)=',minval(cellsOnCell)
        !write(0,*)'after mpas_filter_cells_by_area: maxval(cellsOnCell)=',maxval(cellsOnCell)

        call interp_delaunay ( interp_weights, interp_cells,      &
                               input_var, output_var,             &
                               nVertLevels, nVertLevels, &
                               nCells, nlat, nlon                )

        ! lat_ll is still in radians here.
        !call flip_to_cf_and_circle_smooth( output_var, output_var_cf, nlon, nlat, nVertLevels, lat_ll(:,1), filter_radius_km)
        call flip_to_cf( output_var, output_var_cf, nlon, nlat, nVertLevels)

        ierr = nf_put_var_double(ncid_ll, var_output_id(i), output_var_cf )
        call handle_err(ierr)

        write(0,'(A,2E15.6,A,A)') '  mpas_to_latlon: max/min',maxval(output_var_cf), minval(output_var_cf), ' for ',interp_var_output(i)

        deallocate(input_var)
        deallocate(output_var)
        deallocate(output_var_cf)

    end do

    call handle_err( nf_close(ncid_ll) )
    write(0,*) ' mpas_to_latlon: created '//trim(outfile)
    call handle_err( nf_close(ncid) )

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    subroutine interp_delaunay( interp_weights, interp_cells,    &
                                value_cell, value_interp,        &
                                nVertLevels_cells, nVertlevels_ll, nCells, nlat, nlon )

        implicit none

        integer, intent(in) :: nCells, nVertLevels_cells, nVertlevels_ll, nlat, nlon
        real (kind=RKIND), dimension(nVertLevels_cells, nCells), intent(in) :: value_cell
        real (kind=RKIND), dimension(nVertLevels_ll, nlat, nlon), intent(out) :: value_interp
        real (kind=RKIND), dimension(3, nlat, nlon), intent(in) :: interp_weights
        integer, dimension(3, nlat, nlon), intent(in) :: interp_cells

        integer :: i, j, k

        do j=1,nlon
        do i=1,nlat
        do k=1,min(nVertlevels_ll, nVertlevels_cells)
           value_interp(k,i,j) =   interp_weights(1,i,j)*value_cell(k,interp_cells(1,i,j))  &
                                 + interp_weights(2,i,j)*value_cell(k,interp_cells(2,i,j))  &
                                 + interp_weights(3,i,j)*value_cell(k,interp_cells(3,i,j))
        end do
        end do
        end do

    end subroutine interp_delaunay

    ! moved subroutine flip_to_cf_and_circle_smooth to its own file. Mar 12, 2014

    subroutine compute_interp_weights_all( target_lat, target_lon, nearest_vtx,            &
                                        nCells, nVertices, cellsOnVertex,               &
                                        lat_cell, lon_cell,                             &
                                        interp_weights, interp_cells, nlat, nlon )

        implicit none

        integer, intent(in) :: nlat, nlon
        integer, dimension(nlat,nlon), intent(in) :: nearest_vtx
        real (kind=RKIND), dimension(nlat,nlon), intent(in) :: target_lat, target_lon
        integer, intent(in) :: nCells, nVertices
        integer, dimension(3,nVertices), intent(in) :: cellsOnVertex
        real (kind=RKIND), dimension(nCells), intent(in) :: lat_cell, lon_cell
        real (kind=RKIND), intent(out), dimension(3,nlat,nlon) :: interp_weights
        integer, intent(out), dimension(3,nlat,nlon) :: interp_cells

        real (kind=RKIND), dimension(3) :: weights

        integer :: i, j, k
        real (kind=RKIND), parameter :: eps = 1.e-010_RKIND
        real (kind=RKIND), parameter :: invpower = 1.0_RKIND

        real (kind=RKIND) :: lat0, lon0, lat1, lon1, lat2, lon2, lat3, lon3
        real (kind=RKIND) :: area123, area012, area023, area013, sum_weights

        do j=1,nlon
        do i=1,nlat
        lat0 = target_lat(i,j)
        lon0 = target_lon(i,j)
        !if (j.eq.2) write(0,*) j,i,lon0,lat0,nearest_vtx(i,j)
        lat1 = lat_cell(cellsOnVertex(1,nearest_vtx(i,j)))
        lon1 = lon_cell(cellsOnVertex(1,nearest_vtx(i,j)))
        lat2 = lat_cell(cellsOnVertex(2,nearest_vtx(i,j)))
        lon2 = lon_cell(cellsOnVertex(2,nearest_vtx(i,j)))
        lat3 = lat_cell(cellsOnVertex(3,nearest_vtx(i,j)))
        lon3 = lon_cell(cellsOnVertex(3,nearest_vtx(i,j)))

        if (lon0 < 0) lon0 = lon0+2.*pii
        if (lon1 < 0) lon1 = lon1+2.*pii
        if (lon2 < 0) lon2 = lon2+2.*pii
        if (lon3 < 0) lon3 = lon3+2.*pii

        area123 = triangle_area(lat1,lon1,lat2,lon2,lat3,lon3,1.0_RKIND)
        area012 = triangle_area(lat0,lon0,lat1,lon1,lat2,lon2,1.0_RKIND)
        area023 = triangle_area(lat0,lon0,lat2,lon2,lat3,lon3,1.0_RKIND)
        area013 = triangle_area(lat0,lon0,lat1,lon1,lat3,lon3,1.0_RKIND)

        !
        !  check areas
        !
        weights(1) = area023/area123
        weights(2) = area013/area123
        weights(3) = area012/area123

        sum_weights = weights(1)+weights(2)+weights(3)

        weights(1) = weights(1)/sum_weights
        weights(2) = weights(2)/sum_weights
        weights(3) = weights(3)/sum_weights

        do k=1,3
            interp_cells(k,i,j) = cellsOnVertex(k,nearest_vtx(i,j))
            interp_weights(k,i,j) = weights(k)
        end do

        end do
        end do

    end subroutine compute_interp_weights_all

    subroutine convert_lx(x,y,z,radius,lat,lon)

        implicit none

        real (kind=RKIND) x,y,z,lat,lon,radius

        z = radius*sin(lat)
        x = radius*cos(lon)*cos(lat)
        y = radius*sin(lon)*cos(lat)

    end subroutine convert_lx

    integer function nearest_vertex( target_lat, target_lon, &
                                     start_vertex, &
                                     nCells, nVertices, maxEdges, &
                                     nEdgesOnCell, verticesOnCell, &
                                     cellsOnVertex, latCell, lonCell, &
                                     latVertex, lonVertex )

        implicit none

        real (kind=RKIND), intent(in) :: target_lat, target_lon
        integer, intent(in) :: start_vertex
        integer, intent(in) :: nCells, nVertices, maxEdges
        integer, dimension(nCells), intent(in) :: nEdgesOnCell
        integer, dimension(maxEdges,nCells), intent(in) :: verticesOnCell
        integer, dimension(3,nVertices), intent(in) :: cellsOnVertex
        real (kind=RKIND), dimension(nCells), intent(in) :: latCell, lonCell
        real (kind=RKIND), dimension(nVertices), intent(in) :: latVertex, lonVertex

        integer :: i, cell1, cell2, cell3, iCell
        integer :: iVtx
        integer :: current_vertex
        real (kind=RKIND) :: cell1_dist, cell2_dist, cell3_dist
        real (kind=RKIND) :: current_distance, d
        real (kind=RKIND) :: nearest_distance

        nearest_vertex = start_vertex
        current_vertex = -1
        !write(0,*) 'in nearest_vertex(), target_lat=',target_lat,' target_lon=',target_lon
        do while (nearest_vertex /= current_vertex)
            !write(0,*) 'in nearest_vertex(), nearest_vertex=',nearest_vertex,' current_vertex=',current_vertex
            current_vertex = nearest_vertex
            current_distance = sphere_distance(latVertex(current_vertex), lonVertex(current_vertex), &
                                               target_lat,                target_lon,                1.0_RKIND)
            nearest_vertex = current_vertex
            nearest_distance = current_distance
            cell1 = cellsOnVertex(1,current_vertex)
            cell2 = cellsOnVertex(2,current_vertex)
            cell3 = cellsOnVertex(3,current_vertex)
            cell1_dist = sphere_distance(latCell(cell1), lonCell(cell1), target_lat, target_lon, 1.0_RKIND)
            cell2_dist = sphere_distance(latCell(cell2), lonCell(cell2), target_lat, target_lon, 1.0_RKIND)
            cell3_dist = sphere_distance(latCell(cell3), lonCell(cell3), target_lat, target_lon, 1.0_RKIND)
            if (cell1_dist < cell2_dist) then
                if (cell1_dist < cell3_dist) then
                    iCell = cell1
                else
                    iCell = cell3
                end if
            else
                if (cell2_dist < cell3_dist) then
                    iCell = cell2
                else
                    iCell = cell3
                end if
            end if
            !write(0,*) 'in nearest_vertex(), current_distance=',current_distance,' iCell=',iCell
            do i = 1, nEdgesOnCell(iCell)
                iVtx = verticesOnCell(i,iCell)
                d = sphere_distance(latVertex(iVtx), lonVertex(iVtx), target_lat, target_lon, 1.0_RKIND)
                if (d < nearest_distance) then
                    nearest_vertex = iVtx
                    nearest_distance = d
                end if
            end do
        end do
        !write(0,*) 'in nearest_vertex(), nearest_vertex, nearest_distance=',nearest_vertex,nearest_distance

    end function nearest_vertex

    real (kind=RKIND) function sphere_distance(lat1, lon1, lat2, lon2, radius)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
    !   sphere with given radius.
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

       implicit none

        real (kind=RKIND), intent(in) :: lat1, lon1, lat2, lon2, radius

        real (kind=RKIND) :: arg1

        arg1 = sqrt( sin(0.5_RKIND*(lat2-lat1))**2 +  &
                     cos(lat1)*cos(lat2)*sin(0.5_RKIND*(lon2-lon1))**2 )
        sphere_distance = 2.*radius*asin(arg1)

    end function sphere_distance

    real (kind=RKIND) function triangle_area(p1_lat,p1_lon,p2_lat,p2_lon,p3_lat,p3_lon,radius)

        implicit none

        real (kind=RKIND) p1_lat,p1_lon,p2_lat,p2_lon,p3_lat,p3_lon,radius
        real (kind=RKIND) a,b,c,s,e,tanqe

        a = sphere_distance(p1_lat,p1_lon,p2_lat,p2_lon,radius)
        b = sphere_distance(p2_lat,p2_lon,p3_lat,p3_lon,radius)
        c = sphere_distance(p3_lat,p3_lon,p1_lat,p1_lon,radius)
        s = 0.5_RKIND*(a+b+c)

        ! On rare occasions to due numerical precision limitations or compiler
        ! stuff, the product we take the square root of is slightly negative. It happens
        ! along lat=0 and lon=+/-90. The product has always been really small negative,
        ! like smaller than epsilon(x), so setting it to zero isn't changing it much (I guess).
        ! It doesn't seem to do any harm to set the offending parts of the product to zero.
        ! One might think that taking the square root of zero
        ! would still result in a NaN, but it doesn't seem to. That's why I use
        ! the max() function to fix it. Maybe instead of 0.0_RKIND, I should use
        ! TINY(X).   But I don't know what X should be.  I got this intrinsic
        ! function from www.nsc.liu.se/~boein/f77to90/a5.html It may depend on
        ! the compiler.  This is when I used PGI.  20131031
        ! Ahijevych
        if (tan(0.5_RKIND*s)*tan(0.5_RKIND*(s-a))*tan(0.5_RKIND*(s-b))*tan(0.5_RKIND*(s-c)) .le. 0) then
            write(0,*) 'trying to square root a negative in function triangle_area'
            write(0,*) 'p1_lat=',p1_lat,' p1_lon=',p1_lon
            write(0,*) 'p2_lat=',p2_lat,' p2_lon=',p2_lon
            write(0,*) 'p3_lat=',p3_lat,' p3_lon=',p3_lon
            write(0,*) 'radius = ',radius
            write(0,*) ' a=', a, ' b=',b
            write(0,*) ' c=', c, ' s=', s
            write(0,*) ' s-a=', s-a, ' s-b=', s-b
            write(0,*) ' s-c=', s-c
            write(0,*) ' epsilon(s-a)=', epsilon(s-a)
            write(0,*) 'tan(0.5_RKIND*(s))   = ', tan(0.5_RKIND*(s))
            write(0,*) 'tan(0.5_RKIND*(s-a)) = ', tan(0.5_RKIND*(s-a))
            write(0,*) 'tan(0.5_RKIND*(s-b)) = ', tan(0.5_RKIND*(s-b))
            write(0,*) 'tan(0.5_RKIND*(s-c)) = ', tan(0.5_RKIND*(s-c))
            write(0,*) 'attempting to square root', tan(0.5_RKIND*s)*tan(0.5_RKIND*(s-a))*tan(0.5_RKIND*(s-b))*tan(0.5_RKIND*(s-c))
            !tanqe = sqrt(tan(0.5_RKIND*s)*tan(0.5_RKIND*(max(0.0_RKIND,s-a)))*tan(0.5_RKIND*(max(0.0_RKIND,s-b)))*tan(0.5_RKIND*(max(0.0_RKIND,s-c))))
            triangle_area = TINY(triangle_area)
        else
            tanqe = sqrt(tan(0.5_RKIND*s)*tan(0.5_RKIND*(s-a))*tan(0.5_RKIND*(s-b))*tan(0.5_RKIND*(s-c)))
            e = 4.*atan(tanqe)
            triangle_area = radius*radius*e
        endif

    end function triangle_area

    real (kind=RKIND) function triangle_area_1(p1_lat,p1_lon,p2_lat,p2_lon,p3_lat,p3_lon)

        implicit none

        real (kind=RKIND) :: p1_lat,p1_lon,p2_lat,p2_lon,p3_lat,p3_lon

        real (kind=RKIND) :: x1,y1,z1,x2,y2,z2,x3,y3,z3
        real (kind=RKIND) :: angle1, angle2, angle3

        call convert_lx(x1,y1,z1,p1_lat,p1_lon,1.0_RKIND)
        call convert_lx(x2,y2,z2,p2_lat,p2_lon,1.0_RKIND)
        call convert_lx(x3,y3,z3,p3_lat,p3_lon,1.0_RKIND)

        angle1 = sphere_angle(x1,y1,z1,x2,y2,z2,x3,y3,z3)
        angle2 = sphere_angle(x2,y2,z2,x1,y1,z1,x3,y3,z3)
        angle3 = sphere_angle(x3,y3,z3,x2,y2,z2,x1,y1,z1)

        triangle_area_1 = angle1+angle2+angle3 - pii

    end function triangle_area_1

    character(len=256) function basename(input)
       implicit none
       character (len=256), intent(in) :: input
       integer :: i
       logical :: back
       back = .true.
       i = scan(input, '/', back)
       basename = trim(input(i+1:))

    end function basename

    real (kind=RKIND) function days_since_1970(input)

       USE datetime_module, ONLY:datetime,date2num
       implicit none
       TYPE(datetime) :: a, b
       character (len=*), intent(in) :: input
       character (len=19) :: ccyy_mm_dd_hh_mm_ss
       integer :: iyear, imonth, iday, ihr, imin, isec, i, xtime_id

       integer :: ncid_tmp

       ! Open netCDF file
       ierr = nf_open(input, NF_NOWRITE, ncid_tmp)
       if (ierr.ne.NF_NOERR) then
          write(0,*) ' error opening mpas file '//input
          call HANDLE_ERR(ierr)
       end if
       ! Get variable id of 'xtime'
       ierr = nf_inq_varid(ncid_tmp,'xtime',xtime_id)
       if (ierr.eq.NF_NOERR) then
          ! Read 'xtime'
          ierr = nf_get_var_text(ncid_tmp,xtime_id,ccyy_mm_dd_hh_mm_ss)
       else
          ! 'xtime' not available, so try parsing filename for time.
          write(0,*)NF_STRERROR(ierr)
          write(0,*) ' no xtime var. Parsing filename '//input
          if (index(input,'diagnostics').gt.0) then
              i = index(input,'diagnostics')
              ccyy_mm_dd_hh_mm_ss = input(i+12:)
          else if (index(input,'diag').gt.0) then
              i = index(input,'diag')
              ccyy_mm_dd_hh_mm_ss = input(i+5:)
          else if (index(input,'mpas.output.').gt.0) then
              i = index(input,'mpas.output.')
              ccyy_mm_dd_hh_mm_ss = input(i+12:)
          else
              write(0,*)'problem with days_since_1970'
              write(0,*)'input='//input
              call exit(1)
          endif
       endif
       ierr = nf_close(ncid_tmp)

       write(0,*) '"'//ccyy_mm_dd_hh_mm_ss//'"'
       read(ccyy_mm_dd_hh_mm_ss,'(I4,1x,I2,1x,I2,1x, I2,1x,I2,1x,I2)') iyear, imonth, iday, ihr, imin, isec
       write(0,*) 'iyear=',iyear,'imonth=',imonth,'iday=',iday,'ihr=',ihr,'imin=',imin,'isec=',isec

       a = datetime(iyear,imonth,iday,ihr,imin)
       b = datetime(1970,1,1,0,0)
       !write(0,*) 'a=',a.strftime("%Y %B %d %H %M %S")
       !write(0,*) 'b=',b.strftime("%Y %B %d %H %M %S")
       !write(0,*) 'date2num(a)=',date2num(a)
       !write(0,*) 'date2num(b)=',date2num(b)
       days_since_1970 = date2num(a) - date2num(b)
       !write(0,*) 'days_since_1970=',days_since_1970

    end function days_since_1970

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! FUNCTION SPHERE_ANGLE
    !
    ! Computes the angle between arcs AB and AC, given points A, B, and C
    ! Equation numbers w.r.t. http://mathworld.wolfram.com/SphericalTrigonometry.html
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    real (kind=RKIND) function sphere_angle(ax, ay, az, bx, by, bz, cx, cy, cz)

        implicit none

        real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz, cx, cy, cz

        real (kind=RKIND) :: a, b, c          ! Side lengths of spherical triangle ABC

        real (kind=RKIND) :: ABx, ABy, ABz    ! The components of the vector AB
        real (kind=RKIND) :: ACx, ACy, ACz    ! The components of the vector AC

        real (kind=RKIND) :: Dx               ! The i-components of the cross product AB x AC
        real (kind=RKIND) :: Dy               ! The j-components of the cross product AB x AC
        real (kind=RKIND) :: Dz               ! The k-components of the cross product AB x AC

        real (kind=RKIND) :: s                ! Semiperimeter of the triangle
        real (kind=RKIND) :: sin_angle

        a = acos(max(min(bx*cx + by*cy + bz*cz,1.0_RKIND),-1.0_RKIND))   ! Eqn. (3)
        b = acos(max(min(ax*cx + ay*cy + az*cz,1.0_RKIND),-1.0_RKIND))   ! Eqn. (2)
        c = acos(max(min(ax*bx + ay*by + az*bz,1.0_RKIND),-1.0_RKIND))   ! Eqn. (1)

        ABx = bx - ax
        ABy = by - ay
        ABz = bz - az

        ACx = cx - ax
        ACy = cy - ay
        ACz = cz - az

        Dx =   (ABy * ACz) - (ABz * ACy)
        Dy = -((ABx * ACz) - (ABz * ACx))
        Dz =   (ABx * ACy) - (ABy * ACx)

        s = 0.5*(a + b + c)
        sin_angle = sqrt(min(1.0_RKIND,max(0.0_RKIND,(sin(s-b)*sin(s-c))/(sin(b)*sin(c)))))   ! Eqn. (28)

        if ((Dx*ax + Dy*ay + Dz*az) >= 0.0) then
            sphere_angle =  2.0 * asin(max(min(sin_angle,1.0_RKIND),-1.0_RKIND))
        else
            sphere_angle = -2.0 * asin(max(min(sin_angle,1.0_RKIND),-1.0_RKIND))
        end if

    end function sphere_angle

     SUBROUTINE HANDLE_ERR(STATUS)
     INTEGER STATUS
     IF (STATUS .NE. NF_NOERR) THEN
       PRINT *, NF_STRERROR(STATUS)
       STOP 'Stopped'
     ENDIF
     END SUBROUTINE HANDLE_ERR

end program mpas_to_latlon

max_stars_repo_name: NCAR/MPAS-vortex-tracker

max_issues_repo_name: NCAR/MPAS-vortex-tracker

max_forks_repo_name: NCAR/MPAS-vortex-tracker